dolfin team mailing list archive

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

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.

OK.

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

Unfortunately we can't be sure that the backend doesn't modify the 
non-zero pattern. For example, the matrix might be initialised with a 
sparsity pattern, but if an entry is added to a 'zero' entry, the 
pattern changes.

Garth

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