dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Sunday 28 September 2008 23:08:34 Anders Logg wrote:
> On Sun, Sep 28, 2008 at 10:56:41PM +0200, Johan Hake wrote:
> > On Sunday 28 September 2008 09:51:41 Anders Logg wrote:
> > > On Sat, Sep 27, 2008 at 11:46:42PM +0200, Johan Hake wrote:
> > > > Hello!
> > > >
> > > > Occasionally we have had some discussion about implementing a basic
> > > > aritmetic operators for matrices in DOLFIN. I have been able to add
> > > > -=, +=, +, -, operators, for the PETScMatrix interface, using the
> > > > petsc function MatAXPY. I also implemented both assignment operators.
> > > >
> > > > I added a public function add(A,a), which add a GenericMatrix, A,
> > > > scaled by a, to the present matrix. This is a usefull funtion when
> > > > combining matrices in a linear algebra setting. The function is also
> > > > used to implement the +=, and -=, which then are used to implement +,
> > > > - together with the assignment operator.
> > > >
> > > > Is this something we want? I would at least like it to happen :)
> > >
> > > I think this sounds excellent.
> > >
> > > > Do I miss any important aspects? What about two distributed matrices?
> > > > Will PETSc AXPY take care of this, (supposing of course that the
> > > > matrices have the same number of nonzeros)?
> > >
> > > I think so (if that's what "Collective on Mat" means).
> > >
> > > > If we add the "add" function in the GenericMatrix interface we could
> > > > implement the +=, -= operators directly in the GenericMatrix
> > > > interface, together with both + and - using +=, -=. This is a good
> > > > solution for at least PETSc and Epetra Matrix that do not implement
> > > > its own += and -=, which uBLAS and MTL4 do.
> > > >
> > > > I have it going for PETScMatrix. Should I implement this interface to
> > > > GenericMatrix, and update the other dolfin Matrix classes too? The
> > > > changes in GenericMatrix would be:
> > > >
> > > >   1) remove explicit from the copy constructor, to allow "return by
> > > >   value"
> > >
> > > Why is the copy constructor explicit? I might have forgotten something
> > > but wasn't the conclusion that we should make all constructors
> > > explicit except copy constructors? Martin?
> >
> > In email from 8.7.2008 you come with this conclusion. The email thread
> > ended with your post, and it seems that nothing more happened. Not having
> > too much c++ expearience, it took me two hours to figure out that it was
> > the explicit keyword that made my implementation not work...
> >
> > Anyway, in my upcomming patch I can remove the explicit in the copy
> > constructors.
>
> Fine with me. The only problem I see might be that Martin has some
> obscure (but valid) reason for keeping it. But he'll let us know. :-)

Ok, maybe I wait then. This is only needed for the + and - operators, which 
need to return by value. And this solution is probably not preferable as we 
then create one extra tmp matrix. 

Should we wait with the + and - operators in c++. They can be added in the 
python interface as allready is done for Vectors, see below.

> > > >   2) add virtual add(const GenericMatrix, real a) = 0
> > > >   3) implement +=, -= using the add function.
> > > >   4) implement +, - using the +=, -=
> > > >
> > > > The changes in the other Matrix classes would be:
> > > >
> > > > For PETSc, and Epetra Matrices.
> > > >
> > > >   1) Implement add(A,a)
> > > >   2) Implement operator=(A)
> > > >
> > > > I can make the changes in GenericMatrix, and implement the interface
> > > > for PETSc, Epetra, and I could update GenericMatrix too. If we want
> > > > to use uBLAS's += -=, and I suppose we want, we need to overload
> > > > these functions in the uBLAS interface, with the proposed
> > > > implementation. Are there any similare functionality to PETSc's AXPY
> > > > in uBLAS, for implementing an "add" function?
> > > >
> > > > I suppose uBLAS and MTL4 are similare with respect to implementation.
> > >
> > > Some further comments:
> > >
> > > 1. Name the function axpy(). We already have axpy() in the vector
> > > classes.
> >
> > Will do!
> >
> > > 2. Make sure the same operators are supported for both vectors and
> > > matrices.
> >
> > This means that we just have to add + and - operator, right?
>
> Yes.
>
> > I do spot an assignment operator for scalars in
> > GenericVector. Should this be implemented in the GenericMatrix too?
>
> Probably not. It doesn't really make sense to assign all values of a
> matrix to something (and perhaps not for vectors?).
>
> > > 3. Place all operators in the GenericClasses (as you suggest).
> > >
> > > 4. We also need to make sure that the same operators are supported in
> > > Python, by first ignoring the operators from C++ and then adding the
> > > operators back in Python using axpy(). Some of this is already there
> > > but last I tried using += for vectors I got a segmentation fault and
> > > had to us axpy().
> >
> > In my implementation, I added the operators direct in the PETScMatrix
> > class, and all of them were nicely wrapped to python, and I do not get a
> > segfault using the implemented operators in Vector. Here Vector has a
> > PETScVector, doing the work.
> >
> >    >>> u = Vector(10)
> >    >>> v = Vector(10)
> >    >>> u.assign(10)
> >
> >    <dolfin.dolfin.Vector; proxy of <Swig Object of type 'dolfin::Vector
> > *' at 0x87193b0> >
> >
> >    >>> v.assign(5)
> >
> >    <dolfin.dolfin.Vector; proxy of <Swig Object of type 'dolfin::Vector
> > *' at 0x8191dd0> >
> >
> >    >>> u-=v
> >    >>> u[0]
> >
> >    5.0
> >
> > same with uBLAS, PETSc and Eptra Vector.
>
> Very nice. In that case, just remove the stuff in
>
>   dolfin/swig/dolfin_la_post.i

I looked at it and it seems that what you are refering too is the + and - 
operators. These are implemented as python code in dolfin_la_post. I also got 
a segfault for these operators, while trying to implement them in c++, but 
that was due to the explicit in the copy constructor. When I removed this it 
worked.

Johan