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;
>        }
>    }
> };
> 
> 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. Perhaps the simplest solution is to remove LinearSolver, and do
something like

  KrylovSolver(Preconditioner::LU);

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

>