dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Jun 12, 2006 at 11:49:26AM +0200, Garth N. Wells wrote:
> 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'll take care of the ones in src/kernel/ode at least.

> 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).

Very nice.

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

Don't we need normal derivatives in FFC for this to be useful?

/Anders


> Garth
> 
> > _______________________________________________
> > 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