dolfin team mailing list archive

[Anders Logg <logg@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