dolfin team mailing list archive

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

On Fri, 2006-04-07 at 10:36 -0500, Anders Logg wrote:
> 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... :-)
> 

We can just add an #include at the bottom of FEM.h. Seems to be the
common approach for large templates.

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

What about in the assembly? If DenseMartix is converted to Matrix, won't
there be the overhead problem when inserting values due to virtual
function calls?

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

I didn't know that this works (still learning . . . ). Does the compiler
look in Matrix for a copy constructor from 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 about inserting values in the assembly functions when DenseMatrix
has been been converted to Matrix?

Garth

> 
> > > 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
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev