dolfin team mailing list archive

["Garth N. Wells" <g.n.wells@xxxxxxxxxx>]

On Thu, 2006-04-06 at 17:49 -0500, Anders Logg wrote:
> On Thu, Apr 06, 2006 at 04:33:14PM -0500, Anders Logg wrote:
> > I thought about the linear algebra some more...
> > 
> > If we make the main interface a template, doesn't that mean we need to
> > template all the assembly routines?

Yes, but it only requires 

  template<class NewMatrix>

to be inserted before the function declaration. The annoyance is the
need to include all code in the in the header file. There are some neat
tricks to get around this (see http://www.parashift.com/c
++-faq-lite/templates.html#faq-35.14),
and the keyword "export" is part of the C++ standard which circumvents
this problem, but isn't yet implemented by most compilers, including g
++.

> > 
> > My suggestion would therefore be to do the usual envelope-letter design
> > that we have in other places with a class Matrix as the common
> > interface, but at the same time allow users to work with the
> > DenseMatrix class. If someone needs to do something which is more
> > efficient with the DenseMatrix class or do something which is not
> > supported by the common interface, then use a DenseMatrix.
> > 

I have a few objections to the envelope-letter design for matrices. One
is the overhead for small functions which are called frequently (like
inserting or retrieving entries) due to the virtual functions. This is
not an issue for sparse matrices since entries are added in blocks, but
is important for dense matrices.
I'd like to try using FFC to some of the work I'm currently doing by
hand when constructing small dense matrices. This will require fast
assembly of dense matrices. The other objection is the amount of code
required and the maintenance for the envelope-letter design. It requires
the classes Matrix, GenericMatrix, SparseMatrix and DenseMatrix. Using
the template design we eliminate the need for GenericMatrix.


> > If we then supply constructors in Matrix that accept as arguments
> > either SparseMatrix or DenseMatrix, one can call the assembly routines
> > (which take as input a Matrix) with a DenseMatrix as well as with a
> > Matrix.
> > 
> > This way, there would be no overhead in accessing individual elements
> > (if you work through the DenseMatrix interface) and there would be a
> > common interface Matrix which can be used whenever someone does not
> > have special need for individual accessing.
> > 
> > /Anders
> 
> More thoughts:
> 
> In Matrix we add copy constructors (and assignment) for the types
> SparseMatrix and DenseMatrix, which enables the following:
> 
>     DenseMatrix A;
>     Matrix B(A);
> 
>     SparseMatrix A;
>     Matrix B(A);
> 
> These copy constructors take a (non-const) reference and just use the
> data provided. No copying. This is the same as for the Function class
> where one can create a Function from a Vector (and using the storage
> provided by the given Vector).
> 

Sure. Sounds good.

> We can also create copy constructors (and assignment) in the other
> direction:
> 
>     Matrix A;
>     DenseMatrix B(A);
> 
>     Matrix A;
>     SparseMatrix B(A);
> 
> These would copy the values (not just pick the data).
>
> One should also be able to convert between SparseMatrix and
> DenseMatrix:
> 
>     DenseMatrix A;
>     SparseMatrix B(A);
> 
>     SparseMatrix A;
>     DenseMatrix B(A);
> 
> These also copy the values.
> 
> The thing that might be missing from the picture is how to get the
> DenseMatrix storage from a given Matrix, something like
> 
>     Matrix A;
>     DenseMatrix& B = A.denseMatrixStorage();
> 
> I'd like to say that we don't need this. If someone wants to work with
> dense matrices, then just use DenseMatrix all the time. A DenseMatrix
> can be given as argument to the assembly routines, solvers etc, so
> there is no need to create a Matrix and then get the DenseMatrix
> behind it.
> 

This could be useful. Say a user assembles a Matrix (see next paragraph)
which is dense, and then wants to do some dense matrix-specific
operations which are provided by uBlas. 

Won't passing DenseMatrix to the assembly functions require duplication
of the assembly functions?

  FEM::assemble(BilinearForm&, Matrix&, Mesh&);
  FEM::assemble(BilinearForm&, DenseMatrix&, Mesh&);


> What do think?
> 
> It would be more work to maintain than the template solution, but it
> looks simpler to the user and it doesn't propagate like the template
> solution does (forcing everything to be a template).
> 

If we can find a neat solution to including the definitions of all
templated functions in the header file, I favour the template design for
efficiency and simplicity.

Garth 

> We can also pick a minimum of functions to put in the common interface
> Matrix.
> 
> /Anders
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev