dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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

When it comes to LinearSolver, we might just remove it. I grepped
through the code and it seems we don't use it.

/Anders