dolfin team mailing list archive

[kent-and@xxxxxxxxx]

> 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.
> :-)
>
> --
> Anders


I like it.

Kent