dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Apr 04, 2008 at 12:42:47PM +0200, Kent-Andre Mardal wrote:
> 
> fre, 04.04.2008 kl. 12.03 +0200, skrev Anders Logg:
> > Ola and I have been discussing some of the issues with the linear
> > algebra design.
> > 
> > There are two classes of functions: those that work through the
> > GenericFoo interface, like assemble(), and those that need to resolve
> > the backend, like FooMatrix::mult(). The reason we need to resolve the
> > backend in some functions is that we want to call the backend for
> > efficiency. For example, in PETScMatrix::mult(), we want to call
> > MatMult in PETSc and not reimplement A*x through the GenericFoo
> > interface.
> > 
> > So, if we want to have mult() in the GenericMatrix interface, we need
> > to implement it in the PETScMatrix subclass and try to resolve the
> > PETSc backend. The first iteration would be something like this:
> > 
> >   const PETScMatrix& PETScMatrix::mult(const GenericVector& x,
> >                                        GenericVector& Ax)
> >   {
> >     const PETScMatrix* xx  = dynamic_cast<PETScMatrix*>(&x);
> >     const PETScMatrix* Axx = dynamic_cast<PETScMatrix*>(&Ax);
> > 
> >     if (!xx || !Axx)
> >       error("Incompatible types for matrix-vector multiplication.");
> > 
> >     MatMult(A, xx->vec(), Axx->vec());
> > 
> >     return *this;
> >   }
> > 
> > This works if the arguments are PETScVectors.
> > 
> > It also works if the arguments are uBlasVectors (an exception will be
> > thrown which is the correct behaviour).
> > 
> > Now, the problem is that it won't work if the arguments are Vectors
> > that are wrappers for PETScVectors. Thus, even if you have two Vectors
> > which in principle are PETScVectors, it won't work, but it should.
> > 
> > A simple solution is to add one function to the GenericVector (and
> > other GenericFoo classes):
> > 
> >   // Return instance (implementation)
> >   GenericVector* instance();
> > 
> > In PETScVector, the implementation will be
> > 
> >   PETScVector* instance() { return this; }
> > 
> > while in Vector, the implementation will be
> > 
> >   GenericVector* instance() { return this->vector; }
> > 
> > The second (and final?) iteration of mult() would then be
> > 
> >   const PETScMatrix& PETScMatrix::mult(const GenericVector& x,
> >                                        GenericVector& Ax)
> >   {
> >     const PETScMatrix* xx  = dynamic_cast<PETScMatrix*>(x.instance());
> >     const PETScMatrix* Axx = dynamic_cast<PETScMatrix*>(Ax.instance());
> > 
> >     if (!xx || !Axx)
> >       error("Incompatible types for matrix-vector multiplication.");
> > 
> >     MatMult(A, xx->vec(), Axx->vec());
> > 
> >     return *this;
> >   }
> > 
> > Will this solve all/some problems, or are are there some issues we
> > haven't thought of. I guess Martin will have some nice counter-examples. :-)
> > 
> 
> The consts in the above code is placed a bit random. 

Yes. I'm not really good at coding when gcc is not there to help me.

-- 
Anders