dolfin team mailing list archive

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

On Wed, 2006-07-05 at 09:14 +0200, Anders Logg wrote:
> On Tue, Jul 04, 2006 at 08:39:48PM +0200, Garth N. Wells wrote:
> > On Tue, 2006-07-04 at 11:32 +0200, Garth N. Wells wrote:
> > > On Tue, 2006-07-04 at 11:24 +0200, Anders Logg wrote:
> > > > On Tue, Jul 04, 2006 at 09:12:48AM +0200, Anders Logg wrote:
> > > > > I'm wondering about how to specify the interface for uBlas
> > > > > preconditioners. I think the natural would be
> > > > > 
> > > > >     void solve(DenseVector& x, const DenseVector& b);
> > > > > 
> > > > > but it seems that the Krylov solver prefers an in-place
> > > > > 
> > > > >     void solve(DenseVector& x);
> > > > > 
> > > > > where x = b before the call.
> > > > > 
> > > > > I think we need to have just one interface. There are two
> > > > > options. Either we force the solve(x, b) interface, which I think
> > > > > we can do without sacrificing performance, or we force the solve(x)
> > > > > interface, which might not be natural for all types of
> > > > > preconditioners.
> > > > > 
> > > > > A Jacobi preconditioner would need to store a separate vector and make
> > > > > an additional copy which seems unnessary.
> > > > > 
> > > > > /Anders
> > > > > _______________________________________________
> > > > > DOLFIN-dev mailing list
> > > > > DOLFIN-dev@xxxxxxxxxx
> > > > > http://www.fenics.org/mailman/listinfo/dolfin-dev
> > > > 
> > > > I have modified the uBlas preconditioning interface and moved the ILU
> > > > preconditioner to a separate class. Here are the timings for the
> > > > benchmark, before and after the change:
> > > > 
> > > > Before change
> > > > -------------
> > > > 
> > > > Solving linear system of size 38809 x 38809 (Krylov solver).
> > > > Krylov solver (gmres, ilu) converged in 2851 iterations.
> > > > norm(x) 12.2
> > > > PETSc Krylov solve time = 33.2
> > > > uBlas assembly time = 0.29
> > > > Solving linear system of size 38809 x 38809 (uBlas Krylov solver).
> > > > Krylov solver converged in 2771 iterations.
> > > > norm(x) 12.2
> > > > uBlas Krylov solve time = 26.7
> > > > Factor ( < 1 mean uBlas is faster ) 0.805
> > > > 
> > > > After change
> > > > ------------
> > > > 
> > > > Solving linear system of size 38809 x 38809 (Krylov solver).
> > > > Krylov solver (gmres, ilu) converged in 2851 iterations.
> > > > norm(x) 12.2
> > > > PETSc Krylov solve time = 33.2
> > > > uBlas assembly time = 0.3
> > > > Solving linear system of size 38809 x 38809 (uBlas Krylov solver).
> > > > Krylov solver converged in 2744 iterations.
> > > > norm(x) 12.2
> > > > uBlas Krylov solve time = 27
> > > > Factor ( < 1 mean uBlas is faster ) 0.814
> > > > 
> > > > Some more work (a few extra copies) is made so it might be marginally
> > > > slower than yesterday (26.7 --> 27.0, about 1%).
> > > > 
> > > > Here's the new interface for solving with GMRES and ILU
> > > > preconditioning:
> > > > 
> > > >   uBlasKrylovSolver solver(uBlasKrylovSolver::gmres);
> > > >   uBlasILUPreconditioner pc(A);
> > > >   solver.solve(A, x, b, pc);
> > > > 
> > > > Garth, take a look and see if you like it. I'll keep my hands of the
> > > > linear algebra now for a while so if you want to do the templating and
> > > > renaming of the matrix/vector classes, please go ahead.
> > > > 
> > > 
> > > Looks fine.
> > > 
> > > Garth.
> > > 
> > 
> > One thing - wouldn't it be better if the Krylov solver is preconditioned
> > by default?
> > 
> > Garth
> 
> Maybe, but I don't think it's possible. The Krylov solver can only
> assume that the matrix it gets can do mult() and size(). To create
> a preconditioner, it normally needs to know more than that. The ILU
> preconditioner needs to know a whole lot more.
> 

Why do we have just one uBlasKrylovSolver::solve() function for both
sparse and virtual matrices, whereas we have two
PETScKrylovSolver:solve() functions (one or sparse matrices, one for
virtual matrices)? What if we had something like

    uint solve(const uBlasSparseMatrix& A, DenseVector& x, 
               const DenseVector& b);

    uint solve(const uBlasKrylovMatrix& A, DenseVector& x, 
               const DenseVector& b, const uBlasPreconditioner& pc);

where the first function would apply a default preconditioner. If we
don't precondition by default, I would prefer to remove the function

    uint solve(const uBlasKrylovMatrix& A, DenseVector& x, 
               const DenseVector& b);

because it will be useless for 99% of problems and encourages the
unknowing the use of non-preconditioned solvers.

> One could argue that the Krylov solver should not precondition by
> default. A user needs to choose the Krylov method and the
> preconditioner. The combination of a Krylov solver and preconditioner
> could be something on top of that. For example, the modules will
> choose a suitable combination, and the automated solvers (LinearPDE
> etc) will also choose a suitable combination.
> 

This isn't consistent with the approach we've taken with the PETSc
Krylov solver. I prefer to have "good" default options which can be
changed if necessary for a particular problem. 

Garth 

> We could add an abstraction layer in between, something like
> LinearSolver that takes a matrix and depending on the matrix type (by
> overloading) it chooses a Krylov solver and preconditioner or a direct
> solver if it gets something dense.
> 
> One could also argue that maybe the matrix itself could be allowed to
> suggest a suitable preconditioner. We could add a member function pc()
> to uBlasKrylovMatrix that returns a suitable preconditioner. This
> function could be virtual but not pure virtual and return an
> uBlasDummyPreconditioner by default.
> 
> /Anders
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev