dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Apr 04, 2008 at 12:57:02PM +0200, Martin Sandve Alnæs wrote:
> 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. :-)
> >
> >
> 
> 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" :-)

Yes, it's strange but that's life.

> 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);

I like that.

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

Yes. Anyway, it sounds like you don't have any major objections?
 
> But the "fun" part begins when considering cross-language memory management...

I trust you will fix this. ;-)

-- 
Anders