dolfin team mailing list archive

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

On Mon, 2006-04-24 at 11:01 -0500, Anders Logg wrote:
> On Mon, Apr 24, 2006 at 04:06:42PM +0200, Garth N. Wells wrote:
> > I found a nice compromise solution for dense and sparse matrices. It's
> > often called the Barton-Nackman trick. 
> > 
> > I've added a base class GenericMatrix which is a template. The classes
> > DenseMatrix and Matrix (which is sparse) are then derived from
> > GenericMatrix (this has only been done for DenseMatrix at this stage).
> > The user can create matrices as usual,
> > 
> > DenseMatrix A;
> > Matrix B;
> > 
> > and use any of the public member functions of the respective class
> > (which don't have to be the same). Functions which accept either a dense
> > or sparse matrix then accept GenericMatrix<class T> as an argument, so
> > they have to be templated. GenericMatrix only needs to have member
> > functions that are used by the templated functions. At this stage,
> > GenericMatrix contains only the functions required for assembly.
> 
> I haven't seen the GenericMatrix template yet (doesn't seem to be
> checked in) but I think I get the picture.
>  
> > For example, it is possible to assemble a dense or sparse matrix by
> > 
> > Mesh mesh;
> > BilinearForm a;
> > FEM::assemble(a, A, mesh)
> > FEM::assemble(a, B, mesh)
> > 
> > Take a look in src/test/main.cpp where a sparse and a dense matrix are
> > assembled. I've only templated one function in FEM.h for testing. 
> > 
> > The advantages of this approach are:
> > 
> > 1) Avoids virtual functions which are called from within loops.
> > 2) No changes at the users end when declaring or assembling matrices.
> > 3) GenericMatrix is simple as it needs to contain only member functions
> > which are called by functions which accept both dense and sparse
> > matrices.
> 
> Point 1 and 2 can also be handled by the envelope-letter design
> (without using any templates). In particular, we only call
> Matrix::add() inside the assembly and then the overhead of the extra
> function call is small.
> 

I don't see how the simple use of the envelope-letter design can avoid
virtual function calls. If a user was to create a Matrix of type dense,
then virtual functions would be called for all operations on the matrix.
The only way I see to avoid that is to create a DenseMatrix to operate
on, and then create a Matrix from it which could be passed to the
assembler, which is not very elegant if not hidden from the user and
simple to implement. 

For point 2, sure. This is an advantage over the more standard template
design advocated in NewMatrix.

> But point 3 is very valid. It makes life so much easier for us if we
> don't have to declare every new function n times.
>
> > The disadvantage is that functions which accept dense and sparse
> > matrices need to be templated. I don't think that this is such a big
> > deal, especially since FEM (which is where the function template will be
> > needed) has been cleaned up now.
> 
> It seems this will only affect the implementation of FEM, so it should
> be ok.
> 
> If it turns out we need something similar in another place, we might
> want to have several templates (FEMMatrix, SomethingElseMatrix etc)
> with interfaces that contain only the functionality we need for a
> particular algorithm. Or we throw everything into GenericMatrix.
> 
 
I guess we wait and see how large GenericMatric becomes. Also, the
template is needed only when an algorithm operates with more than one
type of matrix. Most will probably use a single matrix type.

> > Let me know what you think. If everyone thinks that it's OK, I'll go
> > ahead an complete the implementation. Also, any preference for "Matrix"
> > or "SparseMatrix" as names for a sparse matrix? Most problems involve
> > sparse matrices only, so I have a very slight preference for Matrix.  
> 
> How about naming it SparseMatrix (to be consistent) and then have a
> typedef
> 
>     typedef SparseMatrix Matrix;
> 

Sounds good.

Garth

> /Anders
> 
> 
> > Garth
> > 
> > On Mon, 2006-04-24 at 16:02 +0200, DOLFIN wrote:
> > > One or more new changesets pushed to the primary DOLFIN repository.
> > > A short summary of the last three changesets is included below.
> > > 
> > > changeset:   1893:545d61886439a362a658b0a48e5ef2808441320c
> > > tag:         tip
> > > user:        "Garth N. Wells <g.n.wells@xxxxxxxxxx>"
> > > date:        Mon Apr 24 15:54:19 2006 +0200
> > > files:       src/kernel/fem/dolfin/FEM.h src/kernel/la/DenseMatrix.cpp src/kernel/la/Makefile.am src/kernel/la/Makefile.in src/kernel/la/dolfin/DenseMatrix.h src/kernel/la/dolfin/Makefile.am src/kernel/la/dolfin/Makefile.in src/test/main.cpp
> > > description:
> > > Add GenericMatrix class. This templated class acts as a base class for dense and sparse matrices. It provides polymorphism without having to use virtual functions.
> > > 
> > > 
> > > changeset:   1892:779ec524d4b0c72ec15c7f77f09015645d696f62
> > > user:        "Anders Logg <logg@xxxxxxxxx>"
> > > date:        Sun Apr 23 07:44:47 2006 -0500
> > > files:       src/bench/ode/reaction/main.cpp src/bench/ode/wave/bench.log src/bench/ode/wave/bench.py src/bench/ode/wave/main.cpp
> > > description:
> > > Fixes for second multi-adaptive benchmark (wave equation).
> > > Seems to work fine now:
> > > 
> > > CPU time cG(1):  188.093 s
> > > CPU time mcG(1): 45.4958 s
> > > 
> > > 
> > > changeset:   1891:abcc54c4a3f446f20c29cf8eb032bbf5f9fbaa19
> > > user:        "Anders Logg <logg@xxxxxxxxx>"
> > > date:        Sat Apr 22 15:41:44 2006 -0500
> > > files:       src/bench/ode/reaction/bench-domain.log src/bench/ode/reaction/bench-domain.py src/bench/ode/reaction/bench-tol.log src/bench/ode/reaction/bench-tol.py src/bench/ode/reaction/benchutil.py src/bench/ode/reaction/main.cpp src/bench/ode/reaction/parameters-bench.xml
> > > description:
> > > Recompute benchmarks with improved adaptivity.
> > > 
> > > Results look very good:
> > > 
> > > bench-tol: cg/fixed-point
> > > ---------------------------------
> > > 
> > > tol     N       Error           CPU time        Steps   Iterations      Index
> > > --------------------------------------------------------------------------------------
> > > 1.0e-06 1000    2.345e-05       28.050          117089 (1)      4.000 (0.000)   1.0
> > > 5.0e-07 1000    1.173e-05       39.520          165586 (1)      4.000 (0.000)   1.0
> > > 1.0e-07 1000    2.344e-06       71.940          370254 (1)      3.000 (0.000)   1.0
> > > 5.0e-08 1000    1.172e-06       101.730         523615 (1)      3.000 (0.000)   1.0
> > > 
> > > bench-tol: mcg/fixed-point
> > > ---------------------------------
> > > 
> > > tol     N       Error           CPU time        Steps   Iterations      Index
> > > --------------------------------------------------------------------------------------
> > > 1.0e-06 1000    1.771e-05       14.240          1922 (5)        3.990 (1.498)   95.266
> > > 5.0e-07 1000    1.091e-05       23.300          1912 (9)        4.822 (1.544)   138.240
> > > 1.0e-07 1000    1.921e-06       48.050          1929 (7)        4.905 (1.594)   142.619
> > > 5.0e-08 1000    8.979e-07       49.770          1917 (7)        4.131 (1.680)   172.388
> > > 
> > > 
> > > bench-domain: cg/fixed-point
> > > ---------------------------------
> > > 
> > > tol     N       Error           CPU time        Steps   Iterations      Index
> > > --------------------------------------------------------------------------------------
> > > 1.0e-06 1000    2.345e-05       28.100          117089 (1)      4.000 (0.000)   1.0
> > > 1.0e-06 2000    2.233e-05       64.840          117091 (1)      4.000 (0.000)   1.0
> > > 1.0e-06 4000    2.166e-05       101.290         117090 (1)      4.000 (0.000)   1.0
> > > 1.0e-06 8000    2.219e-05       175.110         117089 (1)      4.000 (0.000)   1.0
> > > 1.0e-06 16000   2.242e-05       327.680         117089 (1)      4.000 (0.000)   1.0
> > > 
> > > bench-domain: mcg/fixed-point
> > > ---------------------------------
> > > 
> > > tol     N       Error           CPU time        Steps   Iterations      Index
> > > --------------------------------------------------------------------------------------
> > > 1.0e-06 1000    1.771e-05       13.640          1922 (5)        3.990 (1.498)   95.266
> > > 1.0e-06 2000    1.680e-05       17.310          1923 (5)        3.991 (1.154)   140.461
> > > 1.0e-06 4000    1.643e-05       24.000          1920 (6)        3.991 (1.004)   184.969
> > > 1.0e-06 8000    1.691e-05       33.670          1918 (5)        3.991 (1.004)   218.782
> > > 1.0e-06 16000   1.697e-05       57.850          1919 (5)        3.990 (1.004)   239.992
> > > 
> > > 
> > > Newton benchmarks are temporarily disabled: takes longer
> > > time anyway and there seems to be a bug in the multi-adaptive
> > > Newton solver (takes more iterations than fixed point). Also
> > > changed kmax to 1e-3 in bench-domain (works better with new
> > > choice of threshold). kmax is now same for both benchmars and
> > > set in parameters-bench.xml.
> > > 
> > > 
> > > -------------------------------------------------------
> > > For more details, visit http://www.fenics.org/hg/dolfin
> > > 
> > > _______________________________________________
> > > 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
> >