dolfin team mailing list archive

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

On Wed, 2006-07-12 at 18:50 +0200, Anders Logg wrote:
> On Mon, Jul 10, 2006 at 08:50:09PM +0200, Garth N. Wells wrote:
> > On Mon, 2006-07-10 at 20:22 +0200, Anders Logg wrote:
> > > On Mon, Jul 10, 2006 at 08:06:16PM +0200, Garth N. Wells wrote:
> > > > On Mon, 2006-07-10 at 19:09 +0200, Johan Jansson wrote: 
> > > > > On Mon, Jul 10, 2006 at 02:26:19PM +0200, Garth N. Wells wrote:
> > > > > 
> > > > > > My first implementation of DenseMatrix worked something like this. I
> > > > > > changed it later because it didn't seem very elegant given that uBlas is
> > > > > > written in C++. Also, it's not possible to access member functions of
> > > > > > uBlasMatrixImplementation without a wrapper function or directly working
> > > > > > with the pointer to uBlasMatrixImplementation. 
> > > > > 
> > > > > Why is it a problem to work with the pointer? Isn't it a good idea to
> > > > > work with uBlasMatrix if you want to use the GenericMatrix interface,
> > > > > and to work with uBlasMatrixImplementation (or
> > > > > ublas_sparse_matrix/ublas_dense_matrix) if you want to use the uBlas
> > > > > interface?
> > > > > 
> > > > 
> > > > It might be  possible to first get the pointer to
> > > > uBlasMatrixImplementation before doing anything, but I prefer not to
> > > > (could be because I'm lazy). Besides my laziness, it might not be so
> > > > simple. If we use the same design as PETScSparseMatrix, I'm pretty sure
> > > > that it's not possible to have a pointer
> > > > 
> > > >     uBlasMatrixImplementation* mat;
> > > > 
> > > > if uBlasMatrixImplementation is templated. It needs to be
> > > > 
> > > >     uBlasMatrixImplementation<matrix_type>* mat;
> > > > 
> > > > How do we deal then with the different matrix_type.
> > > 
> > > If we should follow the same approach as for the PETSc wrappers, then
> > > the pointer should be to a uBlas native type, not to another DOLFIN
> > > class.
> > > 
> > > And it should probably a reference (two versions, one const) instead
> > > of a pointer. The pointer in the PETSc case is just because PETSc Mat
> > > is a pointer and the PETSc functions assume to get a pointer.
> > > 
> > > If we decide to do it like this, it shouldn't increase the amount of
> > > code we need to maintain since uBlasMatrix already needs to implement
> > > the interface of GenericMatrix.
> > > 
> > 
> > This looks like a project which is using swig together with uBlas. I
> > don't know enough (actually anything) about swig to know if it's useful.
> > 
> > http://www2.laas.fr/~tlemaire/viewsvn_jafar/?do=log&project=jafarModules&path=/trunk/jmath/include/jmath/ublasMatrix.i&rev=HEAD
> 
> I think it's a good principle to keep things as simple as possible
> when possible (as Johan J has pointed out to me many times before when
> I try to do something too complicated). If we can generate the swig
> interface for the current implementation without too much hassle, then
> fine. 

Johan, any news on wrapping the uBlas functions? It would be good to
sort this out before proceeding too far with the linear algebra clean
up. We could change uBlasMatrix to look like

  template<class Mat>
  class uBlasMatrix : public Variable, public GenericMatrix, 
                      public uBlasKrylovMatrix
  {
    public:
      .
      .

    private:
      Mat* mat;
  }

Would this work with Swig?

If we make this change, we should consider to what extent the
limitations of Swig should affect the future development of DOLFIN. 

Another possibility is to consider Boost.python. This came up on the
DOLFIN mailing list last year when PyDOLFIN started. It seems to offer
better compatibility with C++. After a quick look on the web, I saw that
the guys at Simula have experience with both Swig and Boost.python, so
they could be a very useful resource. There is a oldish presentation
comparing different tools at
http://people.web.psi.ch/geus/talks/europython2004_geus.pdf .


Garth



> Otherwise, maybe it's not so inconvenient to do
> 
>     ( uBlasSparseMatrix A; )
>     ublas_sparse_matrix& AA = A.mat();
> 
> whenever one needs to work with special uBlas operations on a matrix?
> 
> Typically, this should only have to be done in the uBlas
> solvers. Otherwise, the functions in the GenericMatrix interface
> should be enough?
> 
> /Anders
> 
> 
> > 
> > Garth
> > 
> > > /Anders
> > > 
> > > 
> > > > Garth 
> > > > 
> > > > > > 
> > > > > > Would it help if rather than having 
> > > > > > 
> > > > > >    typedef uBlasMatrix<ublas_sparse_matrix> uBlasSparseMatrix;
> > > > > > 
> > > > > > we have 
> > > > > > 
> > > > > >    class  uBlasSparseMatrix : uBlasMatrix<ublas_sparse_matrix> {};?
> > > > > > 
> > > > > > Garth
> > > > > > 
> > > > > 
> > > > > This doesn't help unfortunately. uBlasMatrix<ublas_sparse_matrix> is
> > > > > not instantiable without knowing what ublas_sparse_matrix is. The
> > > > > result is that SWIG ignores it and uBlasSparseMatrix becomes empty.
> > > > > 
> > > > >   Johan
> > > > 
> > > > _______________________________________________
> > > > DOLFIN-dev mailing list
> > > > DOLFIN-dev@xxxxxxxxxx
> > > > http://www.fenics.org/mailman/listinfo/dolfin-dev
> > > _______________________________________________
> > > DOLFIN-dev mailing list
> > > DOLFIN-dev@xxxxxxxxxx
> > > http://www.fenics.org/mailman/listinfo/dolfin-dev
> >