dolfin team mailing list archive

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

On Fri, 2006-06-02 at 12:39 +0200, Anders Logg wrote:
> 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.
> 

Just done some test and this works fine, but the function cannot be a
pure virtual function. 

Garth

> /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
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev