dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, Apr 07, 2006 at 05:52:50PM +0200, Garth N. Wells wrote:
> 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.

Yes, that's a possibility. But it's a fairly large amount of code to
include in every file that depends on FEM.h.

> > > > > 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 would think the overhead is small since the assembly would only
interact with the Matrix through Matrix::add() which is just one
virtual function call. Then the actual addition of the values will be
done in DenseMatrix::add() which can access all the values very
quickly.

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

Yes and no... It seems I was slightly off track. The following works:

  class Foo
  {
  public:
    Foo() : a(3) {}
    int a;
  };

  class Bar
  {
  public:
    Bar(const Foo& foo) : b(2*foo.a) {}
    int b;
  };

  void bar_function(Bar bar)
  {
    cout << "Value of bar.b = " << bar.b << endl;
  }

  Foo foo;
  bar_function(foo);

Here you call bar_function() with a Foo as an argument when
bar_function() expects a Bar. But it works fine since Bar has a
constructor that creates a Bar from a Foo.

But if you modify the constructor in Bar so it takes a Bar&, then gcc
complains:

  error: invalid initialization of reference of type 'Bar&' from
         expression of type 'Foo'

And since we do want to pass by reference (to be able to pick the
pointer) when creating a Matrix from a DenseMatrix, this won't work.

On the other hand, we wouldn't need to duplicate the code but would
have to add something like this:

    void assemble(BilinearForm& a, DenseMatrix& A, Mesh& mesh)
    {
        Matrix AA(A);
        assemble(a, AA, mesh);
    }

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

I would think the overhead is small since we work through the add()
function which does a lot of work (a doubly nested loop) inside
DenseMatrix. It should dominate the extra function call.

/Anders


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

-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/