dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

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.

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