dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Thursday 05 March 2009 20:36:06 Garth N. Wells wrote:
> 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.

Have you experienced that two matrices that do not have the same nonzero 
pattern are not addable using the EpetraExt::MatrixMatrix::Add? At least 
nothing of that sort is mentioned in the Doxygen generated documentation.

If this is right we only have the "problem" for the PETSc backend. Then we can 
change back the GenericMatrix interface to an axpy without 
same_nonzero_patter. The axpy function in PETScMatrix can use MatAXPY with 
DIFFERENT_NONZERO_PATTERN as default, and then we can add a second axpy 
taking the same_nonzero_pattern as argument.

Johan