dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

2008/4/4, 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

We discussed this on monday after Kent implemented vec() that way, and
I'm guessing you won't be convinced by my counterargument: "it feels
strange" :-)

u = v->instance()->instance()->instance()->instance()->instance();

At least I like the name instance() better than vec() for this.

The function vec() should be defined by a convention, not in the interface,
and return a pointer to the underlying backend-specific type.
I.e. a petsc Vec, Epetra_FEVector, ublas_vector.

One could then write "pretty" code like:

uBlasVector * tmp = dynamic_cast<uBlasVector*>(some_generic_vector->instance());
if(!tmp) dolfin_error("...");
ublas_vector & v = *(tmp->vec());

// use v as an ublas_vector with all its functionality

without using multiple inheritance in uBlasVector.


We could make a helper function for this:

ublas_vector * as_ublas_vector(GenericVector & gv)
{
  uBlasVector * tmp = dynamic_cast<uBlasVector*>(gv.instance());
  if(!tmp) dolfin_error("...");
  return tmp->vec();
}

ublas_vector & v = *as_ublas_vector(some_generic_vector);


And we'll need const versions of vec(), instance(), and as_ublas_vector(...).


But the "fun" part begins when considering cross-language memory management...

--
Martin