dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Apr 07, 2006 at 10:29:47AM +0200, Garth N. Wells wrote:
> 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
> ++.

I don't mind the template<class NewMatrix> but moving all the code in
FEM.cpp to FEM.h does not seem to be practical.

And export should probably be avoided if we want to be able to compile
the code... :-)

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

There is no overhead as long as you work with the class DenseMatrix.

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

Agree, it requires more maintenance, but it seems we can't avoid it
(if we can't do the template solution).

> > > 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&);

No, we just need to write the functions for the class Matrix.

Whenever you do

    DenseMatrix A;
    FEM::assemble(a, A, mesh);

you would automatically call the copy constructor for Matrix:

    Matrix::Matrix(DenseMatrix& A) : A(&)
    {
      // Do nothing (but remember not to delete A)
    }

This is very cheap since the constructor will just take the pointer of
the given DenseMatrix.

Then the assembly will assemble into the Matrix calling Matrix::add()
and the overhead should be very small for that function, even with the
extra function call and the resolution of the virtual function.

And you can yourself index all you want in your DenseMatrix A with no
overhead.

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

I don't see how, but if we can figure it out I also prefer the
template solution.

/Anders