dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Aug 11, 2006 at 05:08:38PM +0200, Garth N. Wells wrote:
> 
> 
> Johan Jansson wrote:
> > On Fri, Aug 11, 2006 at 02:42:49PM +0200, Anders Logg wrote:
> >> Could we define a GenericMatrixElement in the same way as we have
> >> defined a PETScMatrixElement to allow access according to
> >>
> >>     a = A(i, j)
> >>     A(i, j) = a
> >>
> >> ?
> >>
> >> This will be *slow* for PETSc, but fast for uBlas. The question is if
> >> we can do this without introducing any overhead when indexing a
> >> matrix that is known to be a uBlasMatrix at compile-time.
> >>
> >> /Anders
> > 
> > I don't like this kind of workaround (i.e. MatrixElement,
> > VectorElement) for just defining an operator. Here's my suggestion:
> > 
> 
> Could we get rid of these classes (I find them confusing)?

Don't think it's possible. They are indeed confusing but I don't know
of another way to differentiate between access and assignment. If
someone does

  a = A(i, j)

then we need to call MatGetValue and if someone does

  A(i, j) = a

then we need to call MatSetValue. The Element classes are a trick to
solve this.

> > Let's define functions getelement()/setelement() for GenericVector and
> > GenericMatrix. The implementation for PETSc will be slow but for dense
> > uBlas very fast.
> >
> 
> The problem is that any functions that are defined in 
> GenericMatrix/GenericVector have some overhead due to the functions 
> being virtual. When accessing elements of a uBlas vector/matrix within a 
> loop this becomes noticeable.

There won't be any overhead as long as you work in terms of a
DenseMatrix/uBlasDenseMatrix/uBlasSparseMatrix. Then the type will be
known at compile-time and the virtual functions don't have to be
resolved.

/Anders

> > We can then define nice operators in Python for these functions with
> > __getitem__()/__setitem__():
> > 
> >     def __setitem__(self, i, a):
> >         if(isinstance(i, int)):
> >             self.setelement(i, a)
> >         elif(isinstance(i, tuple)):
> >             self.setelement(i[0], i[1], a)
> > 
> > In C++, the uBlas indexing interface already exists, so we can use that.
> > 
> 
> Is the issue with SWIG being unable to wrap uBlas a problem here? I can 
> add indexing functions to uBlasMatrix and uBlasVector with no 
> performance penalty if that helps.
> 
> Garth
> 
> > What do you think?
> > 
> >   Johan
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/mailman/listinfo/dolfin-dev
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev