dolfin team mailing list archive

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

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

Probably not... It's very very difficult (or impossible) to make a
failsafe system involving a mixture of different libraries, so I'll
probably be satisfied with keeping references in Python manually
myself where strange situations occur.

--
Martin