dolfin team mailing list archive

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

On Mon, 2006-06-12 at 11:04 +0200, Anders Logg wrote:
> On Mon, Jun 12, 2006 at 10:21:09AM +0200, Garth N. Wells wrote:
> > On Mon, 2006-06-12 at 10:10 +0200, Anders Logg wrote:
> > > On Mon, Jun 12, 2006 at 09:37:59AM +0200, Garth N. Wells wrote:
> > > > On Mon, 2006-06-12 at 08:39 +0200, Anders Logg wrote:
> > > > > On Wed, Jun 07, 2006 at 10:56:18AM +0200, Anders Logg wrote:
> > > > > > 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;
> > > > > > }
> > > > > 
> > > > > Looks good now. The ODE solver works again now after the latest fix.
> > > > > 
> > > > > Should we also move the solve() functions from the uBlas matrix types
> > > > > (DenseMatrix and uBlasSparseMatrix) to the solver classes, so there
> > > > > are no solve() functions in the matrix classes?
> > > > > 
> > > > 
> > > > I think that it should stay in DenseMatrix. It is intended for solving
> > > > very small dense matrices rapidly, and is not suitable for large
> > > > systems. From this point of view, it would be better if it returned a
> > > > vector,
> > > > 
> > > >   x = A.solve(b)
> > > > 
> > > > so it could be used as part of matrix-vector algebra, 
> > > > 
> > > >   real a = inner_prod(A.solve(b), d);
> > > > 
> > > > A solver class should perform more checks, have more options, and for
> > > > the LU solver to be any good, we need to do some renumbering.
> > > > 
> > > > Garth
> > > 
> > > Maybe, but since it can also be used to repeatedly solve (small or
> > > large) dense systems in for example the ODE solver, it would be useful
> > > to have a solver class that can be reused where one supplies the x.
> > > 
> > 
> > Sure. Just as long as there remains a slimmed-down solver which is
> > intended for solving very small problems rapidly. I have a lot of
> > problems which solve regularly 4x4 and 6x6 matrices. If a 6x6 matrix is
> > solved using the "full scale" linear solver, the solve time could be
> > dominated by setting parameters, doing checks, renumbering, etc.
> 
> This sounds like something for a specialized solver class. If one
> needs a solver that is dirty and quick, then use a specialized solver.
> But right now, I just suggest moving the very small amount of code
> from DenseMatrix::solvefoo() to a separate solver class.
> 
> > > I think we should separate data structures from algorithms and
> > > put all solve() functions in solver classes. It's also more consistent
> > > with the rest of the solvers/matrix classes. But that doesn't stop us
> > > from having an x = A.solve(b) in the DenseMatrix class. It would just
> > > need to call the solver class.
> > > 
> > > Here's another suggestion for the linear algebra: let's have the base
> > > class LinearSolver like we have now, and from that inherit two
> > > subclasses KrylovSolver and LUSolver that work just like they do
> > > now. Then we can use the two names GMRES and LU for static versions
> > > that just call the solvers:
> > > 
> > >     static GMRES::solve(A, x, b)
> > >     {
> > >         KrylovSolver solver(KrylovSolver::gmres);
> > >         solver.solve(A, x, b);
> > >     } 
> > > 
> > >     static LU::solve(A, x, b)
> > >     {
> > >         LUSolver solver;
> > >         solver.solve(A, x, b);
> > >     } 
> > > 
> > > Then one can do
> > > 
> > >     GMRES::solve(A, x, b);
> > >     LU::solve(A, x, b);
> > > 
> > > for quick access. For more advanced usage (like reuse of a solver,
> > > setting parameters etc) one creates a solver object like now.
> > > 
> > 
> > Sounds good.
> 
> ok.
> 
> I'm currently keeping pretty busy with the new mesh library, but will
> take a closer look at the linear algebra in a while.
> 
> How about the following plan:
> 
>  - I merge the new mesh library in a day or two (NewMesh) to
>    get some feedback
> 
>  - We work out any remaining issues with the linear algebra library
>    and ODE solvers
>   
>  - We release a very good 0.6.2
> 
>  - We concentrate on porting everything to the new mesh library
> 
> /Anders
> 

We should also remove all possible direct PETSc calls. There a quite a
few which I haven't touched in src/kernel/ode and in src/modules.

I also have a very simple KrylovSolver for the uBlas data types. The
high-level syntax of uBlas makes these very compact. I'll add it
shortly. Most of the demos won't run with the uBlas LU solver (way too
slow).

What about using Neumann boundary conditions in more of the demos?

Garth

> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev