dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jun 06, 2006 at 09:03:12PM +0200, Garth N. Wells wrote:
> On Tue, 2006-06-06 at 20:22 +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 take this back! Just tested it some more and it doesn't work, so
> another solution is required.
> 
> Garth 

There seem to be two options. In either case we need to resolve which
specific implementation of solve() to call and we need to get from
GenericMatrix to a specific matrix type.

The first option is to call GenericMatrix::solve(x, b) in
LinearSolver::solve(). Then GenericMatrix::solve() will resolve to the
solve() function of the specific matrix type, which in turn can call a
specific solver with *this (as in my sketch above).

The other option is to avoid having solve() functions in the matrix
types, but then each matrix type needs to know its type (something
like an enum Type { petsc_sparse, ublas_sparse, ...} etc in
GenericMatrix). Each type of solver, like LU inheriting from
GenericMatrix then needs to check the value of the type, do a cast and
call the specific solver:

class LU : public LinearSolver
{
public:

   void solve(GenericMatrix& A, ...)
   {
       switch ( A.type() )
       {
       case GenericMatrix::petsc_sparse:
           PETScLUSolver solver;
           solver.solve(static_cast<PETScMatrix>(A), ...);
           break;
       }
   }
};

The first option has the advantage that we don't need to introduce the
extra type variable and everything is resolved automatically. The
second option has the advantage that we don't need to put solve
functions into the matrices, but one could argue that only the matrix
types know what solvers are best for them.

/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
> 
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev