dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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

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

-- 
Anders

Attachment: signature.asc
Description: Digital signature