dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Apr 06, 2006 at 06:06:38PM +0200, Garth N. Wells wrote:
> On Thu, 2006-04-06 at 10:34 -0500, Anders Logg wrote:
> > On Thu, Apr 06, 2006 at 05:20:00PM +0200, Garth N. Wells wrote:
> > > I've created a NewMatrix template to allow a matrix to be created which
> > > is either dense or sparse. Some benchmarks are in src/test/main.cpp for 
> > > 
> > >   NewMatrix<DenseMatrix> A;
> > > 
> > > There is no overhead when operating on A with respect to directly
> > > creating a uBlas matrix and operating on it.
> > 
> > Very good!
> > 
> > The only drawback as I see it compared to the envelope-letter design
> > is that we can't change the storage dynamically, but we can probably
> > live without that.
> > 
> 
> What do you mean?

I mean you can declare a matrix as sparse (or dense) and then at some
later point, it will (for some reason) change it's storage from one to
the other.

    Matrix A;    // sparse matrix
    A.toDense(); // A is now dense

For the Function class, this is very convenient. It allows to do
things like

    Function u0 = 0.0; // zero initial condition
    Function u1 = u0;  // zero initial condition
    u1 = ...           // compute new end-point value in time-stepping
    u0 = u1;           // update at left end-point

Here both u1 and u0 change their storage dynamically as needed.

But I think this is less important for matrices.

> > It should be enough if we provide functionality for converting from
> > one type to the other.
> > 
> > > I suggest putting the necessary functions which are common to dense and
> > > sparse matrices into NewMatrix, and changing the name of Matrix to
> > > SparseMatrix, and creating a similar template for vectors. Any
> > > objections or suggestions?
> > > 
> > > Garth
> > 
> > No objections. This sounds like a very good solution.
> > 
> > One of the functions we need to have is Matrix::add() for assembly.
> > 
> 
> Sure. I was thinking that we should just implement a limited number of
> essential functions in NewMatrix to cut down on code and maintenance
> (functions for assembly, display and basic algebra). I'll add a function
> which returns the pointer to the underlying matrix, allowing access to
> the full range of SparseMatrix or DenseMatrix member functions.

A suggestion is to just call the Function data() and make it a
reference:

    template <class T> T& NewMatrix<T>::data()
    {
        return *A;
    }

Then you would get a compile-time error if you don't put the correct
type on the left-hand side:

   Matrix<DenseMatrix> A;
   DenseMatrix B(A.data);

Just a suggestion, maybe you already have a better solution.

> > Is it possible to set a default type? So one would do just
> > 
> >     Matrix A;
> > 
> > to get a sparse matrix and would need to do something different to get
> > a dense one. I'd prefer something else than Matrix<DenseMatrix>,
> > perhaps
> > 
> 
> This is "half" possible. Using
> 
>   template < class T = SparseMatrix >
>   class Matrix . . . 
> 
> a matrix is created by NewMatrix<> A. Not very pretty.

Agree, not very pretty.

> >     Matrix<dense> A;
> > 
> > if possible.
> > 
> 
> I prefer DenseMatrix as a class name, as a user can create a
> DenseMatrix, rather than a Matrix. For my applications, I often need
> dense matrices, but it is rare that I need to use them in the assembly.
> I would create a DenseMatrix so that I could access all the uBlas
> functions.

ok.

> > Or we could rename the DenseMatrix and SparseMatrix classes to
> > something like
> > 
> >     DenseMatrixImplementation / SparseMatrixImplementation 
> >     DenseMatrixStorage / SparseMatrixStorage
> >     DenseMatrixReal / SparseMatrixReal
> >     other options?
> > 
> > and then do
> > 
> >    typedef MatrixInterface<DenseMatrixImplementation> DenseMatrix;
> >    typedef MatrixInterface<SparseMatrixImplementation> SparseMatrix;
> >    typedef MatrixInterface<SparseMatrixImplementation> Matrix
> > 
> > where MatrixInterface is the name of the template NewMatrix after we
> > rename it.
> > 
> 
> Problem with this is that all matrix functions need to be included in
> MatrixInterface, DenseMatrixImplementation and
> SparseMatrixImplementation, even functions that are relevant only for a
> dense or a sparse matrix. This leads to a lot of duplication. On top of
> that, it prevents direct access to uBlas functions for dense matrices as
> the user has to access everything via MatrixInterface.
> 
> This seems complicated and I don't see the need for the extra
> complication.

ok, I agree.

What I don't like is that it would seem that a user needs to write

    Matrix<SparseMatrix> A;

to work with a matrix (assemble, solve linear system).

It would be better if they would only need to write

    SparseMatrix A;

or

    DenseMatrix A;

and then work only with A. Is it enough (and possible) to create copy
constructors in the template Matrix (NewMatrix) that accept
DenseMatrix and SparseMatrix?

  template <class T> Matrix<T>::Matrix(const T& A) : A(&A)
  {
    // Do nothing
  }

Then we would template common functions like assembly so the assembly
functions accept input arguments in the form Matrix<T>. A user could
then give a DenseMatrix or a SparseMatrix as the argument and the
template Matrix would be automatically created.

/Anders


> Garth
> 
> > This would make it easy for the user to create matrices, but it would
> > make life harder for us when we write functions that take a general
> > matrix (sparse or dense) as an argument, but we can probably find a
> > nice solution for that too...
> > 
> > /Anders
> > 
> > 
> > > 
> > > On Thu, 2006-04-06 at 16:50 +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:   1845:4c75b74c95a56c23f1bb6d073c7e8261515a2ffd
> > > > tag:         tip
> > > > user:        "Garth N. Wells <g.n.wells@xxxxxxxxxx>"
> > > > date:        Thu Apr  6 16:49:50 2006 +0200
> > > > files:       configure configure.ac 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/kernel/la/dolfin/dolfin_la.h src/test/main.cpp
> > > > description:
> > > > Create template class NewMatrix.
> > > > 
> > > > 
> > > > changeset:   1844:6b24c4aa352caa578d97cf8eee398793f9a65499
> > > > user:        "Anders Logg <logg@xxxxxxxxx>"
> > > > date:        Wed Apr  5 16:37:57 2006 -0500
> > > > files:       configure configure.ac
> > > > description:
> > > > Remove -Werror when compiling with --enable-optimization
> > > > 
> > > > 
> > > > changeset:   1843:a665e3cd78c7ec75ef59cf6f61f26c24630a1f22
> > > > parent:      1842:05da1868689e5ec15c64172ce8f8d56ab9a65ad0
> > > > parent:      1841:80c6859e328196bae99a9f15d475a3fae5b60f83
> > > > user:        "Garth N. Wells <g.n.wells@xxxxxxxxxx>"
> > > > date:        Wed Apr  5 14:48:52 2006 +0200
> > > > files:       
> > > > description:
> > > > merge
> > > > 
> > > > 
> > > > -------------------------------------------------------
> > > > 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/