dolfin team mailing list archive

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

On Thu, 2006-04-06 at 11:24 -0500, Anders Logg wrote:
> 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.
> 

We can add functions to convert from different types. 

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

Sounds good.

> > > 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 think that it's a good to specify 

  Matrix<dense> A; 
  Matrix<sparse> B; 

as it forces the user to consider what type of matrix they want to use
and hopefully avoid an inappropriate selection.

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

This would be a good function to have, and I think that it is possible.

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

This would work, but couldn't we then get rid of Matrix completely and
just have SparseMatrix and DenseMatrix since the assembly functions will
be templated? We could add
  tyepdef SparseMatrix Matrix;

Garth

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