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
signature.asc
Description: Digital signature
_______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
