dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Jun 02, 2006 at 11:47:26AM +0200, Garth N. Wells wrote:
> On Fri, 2006-06-02 at 09:56 +0200, Anders Logg wrote:
> > On Thu, Jun 01, 2006 at 10:23:26PM +0200, Garth N. Wells wrote:
> > > On Thu, 2006-06-01 at 20:37 +0200, Anders Logg wrote:
> > > > On Thu, Jun 01, 2006 at 11:08:40AM +0200, Garth N. Wells wrote:
> > > > > On Thu, 2006-06-01 at 10:39 +0200, DOLFIN wrote:
> > > > > > One or more new changesets pushed to the primary DOLFIN repository.
> > > > > > A short summary of the last three changesets is included below.
> > > > > > 
> > > > > > changeset:   1981:5dd79b41a3d85cde459ff80e85edd968aeb90f78
> > > > > > tag:         tip
> > > > > > user:        "Garth N. Wells <g.n.wells@xxxxxxxxxx>"
> > > > > > date:        Thu Jun  1 09:33:27 2006 +0200
> > > > > > files:       src/kernel/la/uBlasKrylovSolver.cpp
> > > > > > description:
> > > > > > Add uBlasKrylovSolver.cpp
> > > > > > 
> > > > > > 
> > > > > > changeset:   1980:f2770086f902eff6b19f85f6367d45054a998f09
> > > > > > user:        "Garth N. Wells <g.n.wells@xxxxxxxxxx>"
> > > > > > date:        Thu Jun  1 09:31:16 2006 +0200
> > > > > > files:       src/kernel/la/dolfin/uBlasKrylovSolver.h src/kernel/la/dolfin/uBlasLUSolver.h
> > > > > > description:
> > > > > > Add uBlasKrylovSolver and uBlasLUSolver.
> > > > > > 
> > > > > > 
> > > > > > changeset:   1979:6fef8c669a6a21635a3ff05a341d0b4fa9d2ae11
> > > > > > user:        "Garth N. Wells <g.n.wells@xxxxxxxxxx>"
> > > > > > date:        Thu Jun  1 09:15:15 2006 +0200
> > > > > > files:       src/kernel/la/DenseVector.cpp src/kernel/la/GMRES.cpp src/kernel/la/KrylovSolver.cpp src/kernel/la/LU.cpp src/kernel/la/Makefile.am src/kernel/la/Makefile.in src/kernel/la/PETScKrylovSolver.cpp src/kernel/la/PETScLU.cpp src/kernel/la/dolfin/DenseVector.h src/kernel/la/dolfin/GMRES.h src/kernel/la/dolfin/KrylovSolver.h src/kernel/la/dolfin/LU.h src/kernel/la/dolfin/Makefile.am src/kernel/la/dolfin/Makefile.in src/kernel/la/dolfin/PETScKrylovSolver.h src/kernel/la/dolfin/PETScLU.h src/kernel/la/dolfin/Preconditioner.h src/kernel/la/dolfin/dolfin_la.h
> > > > > > description:
> > > > > > Add prefix "PETSc" to PETSc-based solvers and data structures. "KrylovSolver" and "LU" are typedefs to "PETScKrylovSolver" and "PETScLU" if PETSc is enabled, otherwise they are aliases to uBlasKrylovSolver and uBlasLUSolver.
> > > > > > 
> > > > > > This should make things simple throughout the code. If a PETSc-specific feature is required, use PETScFoo and add an "#ifdef", otherwise just use Foo.
> > > > > > 
> > > > > 
> > > > > I'm having some trouble with the base class LinearSolver which has the
> > > > > virtual function
> > > > > 
> > > > >  virtual uint solve(const SomeMatrix& A, SomeVector& x, const
> > > > > SomeVector& b) = 0;
> > > > >   
> > > > > The problem is that I'd like to derive different solvers from it which
> > > > > accept different matrix and vector types. For example
> > > > > 
> > > > >   PETScLU::solve(const PETScMatrix& A, PETScVector& x, const
> > > > > PETScVector& b)
> > > > >   uBlasLUSolver::solve(const uBlasMatrix& A, uBlasVector& x, const
> > > > > uBlasVector& b)
> > > > >  
> > > > > Any ideas how to do this? I think that the Barton-Nackman trick (where
> > > > > the base class is a template) which we tested for matrices might do the
> > > > > job, at the expense of the complication introduced by the templates. 
> > > > > 
> > > > > Garth
> > > > 
> > > > How about the following:
> > > > 
> > > > In many places, we have adopted the practive that we place algorithms
> > > > in separate classes to be able to separate data structures and
> > > > algorithms and to not clutter the implementation. One example of this
> > > > is the mesh library, where we have the class Mesh and a separate class
> > > > MeshRefinement for mesh refinement. One may thus refine a mesh by
> > > > 
> > > >     MeshRefinement::refine(mesh);
> > > > 
> > > > In addition, we have also added a member function to Mesh, which just
> > > > does MeshRefinement::refine(*this) so one may do
> > > > 
> > > >     mesh.refine();
> > > > 
> > > > We could follow the same principle for linear algebra, so that one may
> > > > do either
> > > > 
> > > >     GMRES::solve(A, x, b);
> > > > 
> > > > or
> > > > 
> > > >     A.solveGMRES(x, b);
> > > > 
> > > > There can also be a default function solve() that calls GMRES for some
> > > > matrix types and LU for other matrix types. What we would need to do
> > > > would be to have a class PETScGMRES with a static function solve():
> > > > 
> > > >     static PETScGMRES::solve(const PETScSparseMatrix& A, ...);
> > > > 
> > > > We would also add a function solveGMRES() to PETScSparseMatrix that
> > > > just calls this function with *this. Furthermore, we would have a
> > > > class GMRES that takes a GenericMatrix:
> > > > 
> > > >     static GMRES::solve(GenericMatrix& A, ...);
> > > > 
> > > > This function would just call GenericMatrix::solveGMRES(), which would
> > > > resolve to PETScSparseMatrix::solveGMRES() in the case of a PETSc
> > > > sparse matrix, which in turn calls PETScGMRES::solve(*this).
> > > >
> > > 
> > > Haven't had a chance to go through the above carefully yet, but it
> > > sounds good. Static functions could fix things. The weakness is that the
> > > parameter system isn't as neat as it isn't, for example, possible (in a
> > > simple way) to associate different tolerances for different solver
> > > objects of the same type. 
> > 
> > How do you mean? It's possible to assign tolerances individually to
> > different solver objects:
> > 
> > GMRES gmres1, gmres2;
> > 
> > gmres1.set("Krylov tolerance", 1e-5);
> > gmres2.set("Krylov tolerance", 1e-7);
> 
> Sure, this how it's done now. But if we make "solve" a static function
> solve and call
> 
>   KrylovSolver::solve()
> 
> we can't set parameters which a specific to that particular solve.

ok, so maybe it shouldn't be static so we can set parameters, but the
above design should work anyway, we just put the following in
PETScSparseMatrix::solve():

    PETScGMRES gmres;
    gmres.solve(*this, x, b);

If one would still like to have a static function, it must look
different from the non-static (they can't be overloaded). So we could
have something like

    static GMRES::solve(A, x, b, ParameterList& parameters);

that is, a static function that takes as additional arguments a list
of parameters. Then one can have different parameter sets instead of
different solver objects. We can add this later if we think it's
necessary.

> > > > When it comes to LinearSolver, we might just remove it. I grepped
> > > > through the code and it seems we don't use it.
> > > > 
> > > 
> > > I tried first not to derive the solver class from LinearSolver, but ran
> > > into problems with NewtonSolver which has a pointer to a LinearSolver.
> > > Hopefully the above will solve this.
> > 
> > Aha, ok. So there might indeed be use for LinearSolver. We also use it
> > in the ODE solvers, typically to pick a suitable linear solver at
> > startup (based on some parameter) and then work through the pointer.
> > 
> > But this should be easy to fix by just letting LinearSolver work with
> > GenericMatrix:
> > 
> >     virtual LinearSolver::solve(const GenericMatrix& A,
> >                                 GenericVector& x,
> >                                 const GenericVector& b);
> > 
> 
> I don't think that this is allowed (?) since GenericMatrix is an
> abstract base class. It didn't work when I tested yesterday, but this
> could be have been due to another problem I had.
> 
> Garth

This is allowed. There is nothing that stops you from having a virtual
function (or any other function) that takes an abstract (interface)
class reference as argument.

/Anders


> 
> > /Anders
> > 
> > > Garth
> > > 
> > > 
> > > > /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