dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

Anders Logg wrote:
> On Thu, Mar 05, 2009 at 08:12:29PM +0100, Johan Hake wrote:
> > On Thursday 05 March 2009 17:13:02 Garth N. Wells wrote:
> > > We have a high-level interface which provides operations like
> > >
> > >    A += B
> > >
> > > for matrices. 
> > For the record: This is not only for the python interface. The c++ operator+= 
> > also use axpy.
> >
> > > This causes some problems as it hides the complexity which 
> > > is inherent in sparse matrices, in particularly with respect to whether
> > > or not the two matrices have the same non-zero pattern. Unless someone
> > > has a good idea as to get around this robustly, I suggest that we remove
> > > these functions for matrices. A user can still do
> > >
> > >      A.axpy(1.0, B)
> > >
> > > or
> > >
> > >      A.axpy(1.0, B, True)
> > >
> > > where the last argument indicates whether or not the sparsity patterns
> > > of A and B are the same (default is False).
> > A naive suggestion: 
> > Is it possible to compute a unique number based on our SparsityPattern? This 
> > number could be stored as a private member of a Matrix? Then only matrices 
> > from the same SparsityPattern will be addable. When ever a Matrix is changed, 
> > by other means than M.init(sp), this number will be set to some default 
> > incompatible number.
> >
> > If we do not find a robust way to do this I am fine with removing the 
> > operators that use axpy, but I really think it is a neat feature which I 
> > would like to keep.
>
> How about having both axpy() and operators like +, += but implemented
> differently.
>
> The axpy() call can make assumptions about the sparsity pattern (or
> take a flag) and use library calls to PETSc and Epetra to compute the
> result.

I would hope that we can make axpy() do the work for +=, =+, etc,
with the assumption that the sparsity patterns do differ. The tricky bit 
is getting this to work with some of the backends, namely Epetera which 
appears to me to be less well documented than PETSc so I couldn't figure 
it out easily.

> The other calls can use getrow() / setrow() and do everything through
> the GenericMatrix interface (so no separate implementation is needed
> for each backend). This will be slower than axpy() but it is very
> convenient to write A + B and sometimes speed is not important. If it
> is, then use axpy().

This is tricky with sparse matrices because we don't want to encourage 
users to unknowingly use operations which are slow (if they know that 
it's slow, that's ok). This is why we removed simple element-wise access 
to sparse matrices and added the uBLAS backend.

Garth

> ------------------------------------------------------------------------
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev