dolfin team mailing list archive

Thread
Date

Re: [HG] 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.

To: Discussion of DOLFIN development <dolfin-dev@xxxxxxxxxx>
From: Anders Logg <logg@xxxxxxxxx>
Date: Mon, 24 Apr 2006 11:01:04 -0500
In-reply-to: <1145887603.8995.39.camel@localhost.localdomain>
Mail-followup-to: Discussion of DOLFIN development <dolfin-dev@xxxxxxxxxx>
Reply-to: Discussion of DOLFIN development <dolfin-dev@xxxxxxxxxx>
User-agent: Mutt/1.5.9i

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.

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.

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

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

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

Follow ups

Re: [HG] 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.
From: Garth N. Wells, 2006-04-24

References

[HG] 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.
From: DOLFIN, 2006-04-24
Re: [HG] 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.
From: Garth N. Wells, 2006-04-24