dolfin team mailing list archive

[Kent-Andre Mardal <kent-and@xxxxxxxxx>]

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. 

Kent