dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Jun 07, 2006 at 10:07:39AM +0200, Garth N. Wells wrote:
> 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;
> >        }
> >    }
> > };
> > 
> > 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
> > 
> 
> I'm starting to think in circles on this point (and not making much
> sense). The original problem was in creating a solver which could be
> either LU or Krylov. I don't need to create a solver at which I can
> throw any type of matrix or vector, just the "default" Matrix and
> Vector.

I think the ideal would be to have something that works for any matrix
type. I'd like to be able to choose matrix type and solver separately:
use this matrix type and use this solver. Not all solvers will support
all matrix types of course.

> Perhaps the simplest solution is to remove LinearSolver, and do
> something like
> 
>   KrylovSolver(Preconditioner::LU);

I don't think this is the best option, since I'd like to be able to
say "LU" without needing to think about it as a preconditioner for
GRMES.

> or alternatively 
> 
>   class NewtonSolver
>   {
>     public:
>       NewtonSolver();  
>       uint solve(...);
>       ...	
>     private:
>        LU* lu;
>        KrylovSolver* krylov_solver;
>   };
> 
>   uint NewtonSolver::NewtonSolver(...) : krylov_solver(0)
>   {
>     lu = new LU;
>   }
>   uint NewtonSolver::solve(...)
>   {
>     if( lu )
>       lu.solve(...);
>     if( gmres )
>       krylov_solver.solve(...); 
>   }
> 
> It's useful if a solver is not re-created at every solve as PETSc has
> some optimisation features for if the matrix elements do not change or
> if the matrix has the same non-zero structure. 
> 
> Garth

I don't think this is a good solution, since one would like to make
the choice once and not on every solve.

I think a simple solution would be to overload on matrix types in
LinearSolver, with virtual functions that return an error message if
the matrix type is not supported. See the attached sketch. It
compiles, it's very simple and I think it has the generality we need.

/Anders

// Test how to design the linear algebra solvers in DOLFIN.
// Anders Logg, 2006-06-07.

#include <iostream>

using namespace std;

class SparseMatrix {};

class DenseMatrix {};

class LinearSolver
{
public:

  virtual void solve(SparseMatrix& A)
  {
    cout << "Error, not implemented for sparse matrices." << endl;
  }

  virtual void solve(DenseMatrix& A)
  {
    cout << "Error, not implemented for dense matrices." << endl;
  }

};

class KrylovSolver : public LinearSolver
{
public:

  void solve(SparseMatrix& A)
  {
    cout << "Solving sparse system with Krylov solver." << endl;
  }

};

class LUSolver : public LinearSolver
{
public:

  void solve(DenseMatrix& A)
  {
    cout << "Solving dense system with LU solver." << endl;
  }

};

void foo(LinearSolver& solver)
{
  SparseMatrix A;
  DenseMatrix B;

  solver.solve(A);
  solver.solve(B);
}

int main()
{
  LinearSolver* lu = new LUSolver();
  LinearSolver* krylov = new KrylovSolver();

  foo(*lu);
  foo(*krylov);

  delete lu;
  delete krylov;

  return 0;
}