dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Sunday 28 September 2008 23:26:53 kent-and@xxxxxxxxx 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.
> >
> >> >   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? I do spot an
> > assignment operator for scalars in GenericVector. Should this be
> > implemented
> > in the GenericMatrix too?
> >
> >> 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.
>
> The problem we have had is that the object is no longer the same after +=
> You can check this with
> a = ..
> print id(a)
> a += ...
> print id(a)
>
> If the object is new then something has to be deleted.
>
> This somethimes causes troubles.

Yes, now I remember. In PyCC it caused a segfault long time ago. It seems that 
this is not the case in pyDOLFIN at least, but the case with changed id is 
present. 

To handle this problem we need to %ignore them in la_pre and add them as 
python operators in la_post.

Johan