dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jun 06, 2006 at 08:22:21PM +0200, Garth N. Wells wrote:
> On Tue, 2006-06-06 at 20:10 +0200, Anders Logg wrote:
> > On Tue, Jun 06, 2006 at 04:41:31PM +0200, Garth N. Wells wrote:
> > > On Tue, 2006-06-06 at 13:50 +0200, Anders Logg wrote:
> > > > On Tue, Jun 06, 2006 at 12:51:58PM +0200, Garth N. Wells wrote:
> > > > > On Sat, 2006-06-03 at 10:26 +0200, Anders Logg wrote:
> > > > > 
> > > > > > > > > 
> > > > > > > > > Just done some test and this works fine, but the function cannot be a
> > > > > > > > > pure virtual function. 
> > > > > > > > 
> > > > > > > > Strange. The following compiles for me:
> > > > > > > > 
> > > > > > > > class Abstract
> > > > > > > > {
> > > > > > > > public:
> > > > > > > > 
> > > > > > > >   virtual void foo() = 0;
> > > > > > > > 
> > > > > > > > };
> > > > > > > > 
> > > > > > > > class Foo
> > > > > > > > {
> > > > > > > > public:
> > > > > > > > 
> > > > > > > >   virtual void bar(Abstract& abstract) = 0;
> > > > > > > > 
> > > > > > > > };
> > > > > > > > 
> > > > > > > > class Bar : public Foo
> > > > > > > > {
> > > > > > > > public:
> > > > > > > > 
> > > > > > > >   void bar(Abstract& abstract) {}
> > > > > > > > 
> > > > > > > > };
> > > > > > > > 
> > > > > > > > int main()
> > > > > > > > {
> > > > > > > >   Bar bar;
> > > > > > > > 
> > > > > > > >   return 0;
> > > > > > > > }
> > > > > > > > 
> > > > > > > > Is this not what you mean?
> > > > > > > > 
> > > > > > > 
> > > > > > > This doesn't work for me. Can you send the actual code? I'll try it out.
> > > > > > > 
> > > > > > > Garth
> > > > > > 
> > > > > > The code above is the actual code. I've attached it. Just compile with
> > > > > > with
> > > > > > 
> > > > > >     g++ main.cpp
> > > > > > 
> > > > > > Works for me with g++ 4.0.3.
> > > > > > 
> > > > > > /Anders
> > > > > > 
> > > > > > 
> > > > > 
> > > > > This works for me, but the problem I have is slightly different. The
> > > > > problem I have is that the equivalent of Bar::bar (a solver) does't want
> > > > > "Abstract" as an argument, rather "SomeAbstract" type. I've pasted some
> > > > > code below. The solution seems be not to make the functions in
> > > > > LinearSolver pure virtual functons.
> > > > > 
> > > > > Garth
> > > > > 
> > > > > 
> > > > > class GenericMatrix
> > > > > {
> > > > > public:
> > > > >   
> > > > >   virtual void foo() = 0;
> > > > >   
> > > > > };
> > > > > 
> > > > > 
> > > > > class SomeMatrix : public GenericMatrix
> > > > > {
> > > > > public:
> > > > >   
> > > > >   
> > > > > };
> > > > > 
> > > > > class LinearSolver
> > > > > {
> > > > > public: 
> > > > >   
> > > > >   virtual void bar(GenericMatrix& matrix) = 0;
> > > > >   
> > > > > };
> > > > > 
> > > > > class LU : public LinearSolver
> > > > > {
> > > > > public:
> > > > >   
> > > > >   // This doesn't work
> > > > >   void bar(SomeMatrix& matrix) {}
> > > > > 
> > > > >   // This works
> > > > > //  void bar(GenericMatrix& matrix) {}
> > > > >   
> > > > > };
> > > > > 
> > > > > int main()
> > > > > {
> > > > >   LinearSolver* solver;
> > > > > 
> > > > >   solver = new LU;
> > > > >   
> > > > >   return 0;
> > > > > } 
> > > > 
> > > > ok, I see the problem.
> > > > 
> > > > We will need to do something like the following:
> > > > 
> > > >   class LinearSolver
> > > >   {
> > > >   public: 
> > > >     virtual void solve(GenericMatrix& A, ...) = 0;
> > > >   };
> > > > 
> > > >   class GenericMatrix
> > > >   {
> > > >   public:
> > > >     virtual void solveLU(Vector& x, const Vector& b)
> > > >     {
> > > >       dolfin_info("Informative error message, not implemented.");
> > > >     }
> > > >   }
> > > > 
> > > >   class LU : public LinearSolver
> > > >   {
> > > >   public: 
> > > >     void solve(GenericMatrix& A, ...)
> > > >     {
> > > >       A->solveLU(x, b);
> > > >     }
> > > >   };
> > > > 
> > > >   class PETScSparseMatrix : public GenericMatrix
> > > >   {
> > > >   public:
> > > >     void solveLU(Vector& x, const Vector& b)
> > > >     {
> > > >       PETScLUSolver solver;
> > > >       solver.solve(this, x, b);
> > > >     }
> > > >   };
> > > > 
> > > > PETScLUSolver does not inherit from LinearSolver. It could inherit
> > > > from PETScLinearSolver if we want.
> > > > 
> > > 
> > > A simpler solution is to not make the virtual function in LinearSolver
> > > pure. This generates a compile-time error if the declared solver doesn't
> > > accept the matrix/vector types given to it.
> > 
> > But that seems to make the class LinearSolver unusable? Then you can't
> > do
> > 
> > LinearSolver* solver = new LU();
> > solver.solve(A, x, b);
> > 
> > since LU does not need to implement solve() with the matrix type of A.
> > 
> 
> It does work, so long the function "solve" is not a pure virtual
> function. I've just checked it in so try it out.

I just checked. It looks to me like the meaning with LinearSolver is
lost if solve is not pure virtual. If the solvers that inherit from
LinearSolver don't implement a single function specified in the
interface, then in practice they are not subclasses of LinearSolver.

With the current solution, one can't do the following:

// Choose type of solver in some way
LinearSolver* solver;
if ( ... )
  solver = new LU();
else
  solver = new GMRES();

solver->solve(A, x, b);

The latest update also breaks the ODE solvers where this is used:

logg@gwaihir:~/[snip]/dolfin-dev/src/demo/ode/courtemanche$ ./dolfin-courtemanche 
Creating ODE of size 21.
Solving ODE on time interval (0, 300).
  Initializing continous Galerkin method cG(1).
  Initializing PETSc (ignoring command-line arguments).
  Using mono-adaptive Newton solver.
  Using discrete tolerance tol = 0.01.
  Using iterative linear solver: GMRES (default).
  *** Error: Solver not available for matrix/vector type.
  [./dolfin/LinearSolver.h:51: solve()]

/Anders


> > > > Thus, all the work takes place in the specialized classes. We have
> > > > PETScSparseMatrix and we have PETScLUSolver (or PETScSparseLUSolver?).
> > > > 
> > > > The classes LinearSolver and LU, GMRES etc are just very small classes
> > > > that we add for convenience.
> > > 
> > > The matrix classes will become very large. I think that we should keep
> > > the solvers separate from the matrix classes because the solvers are
> > > independent of whether a matrix is dense or sparse in the case of PETSc
> > > and uBlas.  
> > 
> > Yes, I think they should be separate. But is it possible without
> > working putting in a small function in the matrix interface that
> > chooses the correct solver (or by templating)?
> > 
> 
> Sorry, I wasn't thinking too straight on this point before. It would be
> nice and clean to keep the solvers completely separate from the matrix
> classes though. 
> 
> Garth
> 
> > Another option is to overload on all possible matrix types in
> > LinearSolver and LU, and choose the specific solver there.
> > 
> > /Anders
> > 
> > 
> > > 
> > > Garth
> > > 
> > > > 
> > > > When it comes to parameters, we could maybe add a ParameterList as an
> > > > extra parameter to the solve function in GenericMatrix. Thus,
> > > > LinearSolver would inherit from Parametrized so one can do
> > > > 
> > > >     Matrix A; // defaults to PETScSparseMatrix
> > > >     LU solver;
> > > >     solver.set("...", ...);
> > > >     solver.solve(A, x, b);
> > > > 
> > > > The last call would call LU::solve() which calls
> > > > PETScSparseMatrix::solve() which calls PETScLUSolver::solve().
> > > > 
> > > > /Anders
> > > > 
> > > > _______________________________________________
> > > > DOLFIN-dev mailing list
> > > > DOLFIN-dev@xxxxxxxxxx
> > > > http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> > > 
> > > 
> > > _______________________________________________
> > > DOLFIN-dev mailing list
> > > DOLFIN-dev@xxxxxxxxxx
> > > http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> > 
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> 
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev