dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Apr 06, 2006 at 07:42:03PM +0200, Garth N. Wells wrote:
> On Thu, 2006-04-06 at 12:11 -0500, Anders Logg wrote:
> > > 
> > >   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.
> > 
> > This looks pretty neat. On the other hand, we don't have much <>
> > anywhere else in the code (at least not in the user interface) so
> > I don't have a clear opinion.
> > 
> >     SparseMatrix A;
> > 
> >     Matrix<sparse> A;
> > 
> > Which looks best?
> > 
> 
> I think it comes down to whether we want a templated matrix class, or
> templated functions in other classes that work with matrices. The
> advantage of a templated matrix class is that we localise the work -
> just one class - and don't have to touch the code in FEM.cpp. The only
> disadvantage I see is the introduction of <> into the user part of the
> code. It could however be useful in the future for PDE's among other
> things. The problem I see with the letter-envelope design is when
> different PDE's are not easily distinguishable on the basis of the
> arguments passed to the constructor. We could have a neater PDE
> implementation with less code if we could do
> 
>   PDE<linear> pde_linear(BilinearForm&, LinearForm&,
> BoundaryConditions&, Mesh&);
>   PDE<nonlinear>
> pde_nonlin(BilinearForm&,LinearForm&,BoundaryConditions&,Mesh&);

Looks good to me. Just because we haven't used much <> in the user
interface before does not mean we must avoid it.

> 
> > > > 
> > > > 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.
> > 
> > ok, we'll have to try...
> > 
> > > > 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;
> > 
> > Yes, you're right... maybe that's the simplest solution? We would then
> > have
> > 
> >     Matrix       (typedef of SparseMatrix)
> >     SparseMatrix (current Matrix)
> >     DenseMatrix  (current DenseMatrix)
> > 
> > Sounds very simple. We wouldn't need the template matrix class and no
> > class hierarchy... Maybe since it is the simplest solution it is also
> > the best?
> > 
> > But maybe the drawback is that we would need to put all the assembly
> > code in FEM.h. It's a lot of code...
> 
> and anywhere else both dense and sparse matrices might be or are used.

I'd like to avoid that if possible...

So should we go with

    Matrix<dense> A:
    Matrix<sparse> A;

?

In addition, we do the following:

1. sparse is the default type, so Matrix<> A works.

2. a suitable copy constructor in Matrix so DenseMatrix and
   SparseMatrix can be passed as arguments to functions that want a
   Matrix<T> as argument

3. copy constructors in DenseMatrix and SparseMatrix so we can create
   one from the other

If 2 works, would that remove the need for Matrix<dense> and
Matrix<sparse>? One would then always work with DenseMatrix or
SparseMatrix and pass these as arguments to assembly functions.

/Anders