dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Apr 24, 2006 at 06:41:25PM +0200, Garth N. Wells wrote:
> 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. 

I just mean there will be very little overhead since the assembly only
interacts with the matrix through the function Matrix::add() and the
work that add() does should completely dominate the overhead of the
function call.

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

Sounds good. I have no more objections.

/Anders


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