dolfin team mailing list archive

["Garth N. Wells" <g.n.wells@xxxxxxxxxx>]

On Tue, 2006-06-06 at 21:18 +0200, Anders Logg wrote:
> 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;
>        }
>    }
> };
> 

Couldn't me make it straightforward by using template specialisation?
Would something like the below work?

class LU : public LinearSolver
{
public:

   template<class MAT>
   void solve(MAT& matrix, . . . . .)
     { dolfin_error("No solver.."); }

   template<>
   void solve<PETScSparseMatrix>(PETScSparseMatrix& matrix, . . . )
     { PETScLUSolver solver;
       solver.solve(matrix, . . . );
     }

   template<>
   void solve<uBlasSparseMatrix>(uBlasSparsematrix& matrix, . . . )
     { uBlasLUSolver solver;
       solver.solve(matrix, . . . );
     }
};

Garth

> 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
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
-- 
Dr. Garth N. Wells
Faculty of Civil Engineering and Geosciences
Delft University of Technology
Stevinweg 1
2628 CN Delft
Netherlands

tel.   +31 15 278 7922
fax.   +31 15 278 6383
e-mail g.n.wells@xxxxxxxxxxx
url    http://www.mechanics.citg.tudelft.nl/~garth