dolfin team mailing list archive

[kent-and@xxxxxxxxx]

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

Kent