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