dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Aug 15, 2006 at 03:33:21PM +0200, Johan Jansson wrote:
> On Mon, Aug 14, 2006 at 02:42:15PM +0200, Anders Logg wrote:
> > 3. Are all the functions in GenericMatrix/Vector necessary? Could some
> > be removed? Also, there are plenty of functions in PETScMatrix/Vector
> > and uBlasMatrix/Vector that are not in GenericMatrix/Vector. Should we
> > make sure that the interface stays the same for all linear algebra
> > types so for example only GenericMatrix functions are present in the
> > public interface of PETScMatrix?
> 
> I've been thinking a bit about this, but I haven't found a clear
> solution. The problem with using the Generic* interface is that it can
> only support unary operations (i.e. indexing, assembling, ...). If we
> define a binary operation, then we allow cases like adding an
> uBlasVector to a PETScVector, and we don't have that semantic.
> 
> What we should do is define a simple interface for binary operations,
> and then manually enforce that in uBlasVector/Matrix and
> PETScVector/Matrix. And strip away essentially everything else.

Agree, but how do we implement it?

One way would be to make Matrix/Vector envelope classes like we had
before. Then we could do A += B etc with all possible combinations of
matrix types (PETSc, uBlas, sparse or dense). The current design would
remain basically the same, with the addition of a simple envelope
class on top that holds the pointer to GenericMatrix.

Is it possible to add the binary operators without adding the envelope
class? If we just work through GenericMatrix, we would either need to
define

 virtual const GenericMatrix& operator+= (const GenericMatrix& A) = 0;

in GenericMatrix. This could be implemented for example by working
through functions like getRow() and add(block). Is this what you mean?
The disadvantage would then be that we wouldn't use the optimized
versions of += for uBlas and PETSc.

> I guess not much more than the standard vector space operations are
> needed, i.e. multiply by scalar, addition, dot product. But it's
> important that we standardize that. Then if you need something more
> exotic, you use the PETSc/uBlas interfaces.

Agree. Here's a list of suggestions:

    +=
    -=
    *=
    dot()
    norm()

> > 4. Indexing in a PETSc matrix does not work as expected:
> > 
> >   >>> from dolfin import *
> >   >>> A = Matrix(10, 10)
> >   >>> print A[(5,5)]
> >   [0]PETSC ERROR: MatGetValues() line 1324 in src/mat/interface/matrix.c
> >   [0]PETSC ERROR: Object is in wrong state!
> >   [0]PETSC ERROR: Not for unassembled matrix!
> >   1.27319747458e-313
> > 
> > The natural would be for this to return zero.
> > 
> > /Anders
> 
> This is just because a PETSc matrix is not assembled by default. If you do:
> 
> >>> from dolfin import *
> Importing DOLFIN
> >>> A = Matrix(10, 10)
> Initializing PETSc (ignoring command-line arguments).
> >>> A.apply()
> >>> A[5, 5]
> 0.0
> 
> then it works as expected. We should just add apply() in the
> constructor.
> 
>   Johan

ok, let's add it.

/Anders