← Back to team overview

dolfin team mailing list archive

Re: uBlas Krylov solver and preconditioner

 

On Tue, Jul 04, 2006 at 08:26:47AM +0200, Garth N. Wells wrote:
> On Mon, 2006-07-03 at 17:31 +0200, Anders Logg wrote:
> > I've extended the uBlas Krylov solver so that it now works with
> > matrices of type uBlasKrylovMatrix, only defined by having two
> > functions size() and mult(). I had to make some minor changes to avoid
> > using axpy_prod() but the tests I made showed now reduction in
> > speed. (There is one extra copy of values into a vector compared to
> > the original code.)
> > 
> 
> Make sure that you use noalias() whenever possible when working with
> uBlas data types. In the code, r += b should be noalias(r) += b or
> r.plus_assign(b). Otherwise, a copy of r is make internally which gets
> expensive for large vectors. 

ok, I'll fix.

> > I have temporarily disabled the preconditioning. See all comments
> > 
> >     // FIXME: Preconditioner temporarily disabled
> > 
> > This was necessary since uBlasPreconditioner assumes uBlasSparseMatrix.
> > We need to make it so that it only assumes uBlasKrylovMatrix. Then we
> > can have a special subclass, like
> > 
> >     class uBlasILUPreconditioner : public uBlasPreconditioner
> > 
> > or similarly. This class should then take in the constructor as
> > argument a uBlasSparseMatrix. Then we can overload solve()
> > uBlasKrylovSolver so we have two versions: one with and one without an
> > extra preconditioner argument.
> > 
> > I don't think it's possible to let the Krylov solver create a
> > preconditioner on its own since it only knows how to do mult(), that
> > is, it only knows that the matrix is uBlasKrylovMatrix. So it can't
> > create a preconditioner that assumes it can do M.assign(A).
> > 
> > Does that make sense?
> > 
> 
> Sounds ok. The preconditioner interface does need to be sorted out. My
> first implementation was very simple just to get something working.
> 
> Garth  

ok, I'll make an attempt to get the preconditioning going again and
then you can go ahead with the templating. :-)

/Anders


References