dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Mar 05, 2009 at 07:36:06PM +0000, 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.

ok.

> 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.

I think it's ok to provide a simple but slow alternative if it is not
possible to make the fast option simple, as long as it is well
documented. Simple linear algebra operations like indexing, addition
etc are very useful for debugging and writing simple test programs. I
don't agree that everything that is potentially slow must be
prohibited.

-- 
Anders

Attachment: signature.asc
Description: Digital signature