dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Jan 04, 2008 at 06:36:41PM +0000, Garth N. Wells wrote:
> 
> 
> Martin Sandve Alnæs wrote:
> > Would it be ok now to add arithmetic operations in the Generic*
> > interfaces? They can be limited to objects of the same type, using a
> > type assertion like init already must do for the sparsity pattern. A
> > fairly small set of operations will allow us to write many linear
> > algebra algorithms in PyDolfin for any backend, great for
> > application-specialized algorithms and rapid prototyping.
> > 
> > Vectors:
> > v+u
> > a*v
> > 
> > Matrices:
> > A+B
> > a*A
> > 
> > Matrix-vector:
> > A*v
> > 
> 
> I can't comment on integration pyDOLFIN, but we went to uBlas to get 
> this type of functionality for free. The problem we had with adding a 
> lot of operators to Generic* is that the code for each back-end really 
> grows, especially if the back-end is not well suited to the operation 
> (this was often the case for PETSc).
> 
> Garth

I think being able to add vectors and a few other basic operations is
something we need to support.

I propose we add a few basic operations to GenericMatrix and
GenericVector but implement them as *functions*, not operators. For
example, we can have

  void GenericMatrix::mult(const GenericVector& x, GenericVector& Ax) const;

and a few others. Implementing these as functions rather than
operators makes it possible to later add some syntactic sugar, which
we need to do differently in C++ and Python. SWIG will wrap the
functions for us and then we can do some tricks on the Python side.

The syntactic sugar can be added to the classes Matrix and Vector and
not be part of the GenericFoo interfaces.

> > If there are any good reasons to avoid this (f.ex. virtual operators
> > in C++ doesn't exist), I suggest either defining separate interfaces
> > for arithmetic operations (f.ex. a GenericOperator could be something
> > with a matrix-vector product) or at the very least defining
> > conventions for these operations, which would allow PyDolfin apps to
> > use this by "duck-typing" (since Python is dynamically typed).

The reason to avoid it is as Garth says that it's difficult to
implement some of these for PETSc and for uBLAS we already have them
for free. But if we limit ourselves to a few operators that PETSc
supports and don't add the operators to GenericFoo, then we don't have
any problems with PETSc and the expression templates of uBLAS can
still do their job.

-- 
Anders