dolfin team mailing list archive

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

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.

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

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