dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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

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

/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