dolfin team mailing list archive

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

On Mon, 2006-07-10 at 14:08 +0200, Johan Jansson wrote: 
> Hi,
> 
> I've looked a bit at SWIG and uBlas. I started looking at the linear
> algebra interface from a few days back, and had some success in
> wrapping it. However, the change to make uBlasMatrix a template class
> makes wrapping difficult.
> 
> Previously an uBlas matrix in DOLFIN was defined like so:
> 
> class uBlasSparseMatrix : public Variable,
>                           public GenericMatrix,
>                           public uBlasKrylovMatrix,
>                           public ublas_sparse_matrix
> 
> If SWIG is not aware of uBlas, it simply ignores the part of the
> interface coming from ublas_sparse_matrix, and can still generate a
> wrapper for the rest of uBlasSparseMatrix.
> 
> However, an uBlas matrix in DOLFIN is now defined like so:
> 
> template< class Mat >
> class uBlasMatrix : public Variable,
>                     public GenericMatrix,
>                     public uBlasKrylovMatrix,
>                     public Mat
> 
> with the instantiation in the sparse case:
> 
> typedef uBlasMatrix<ublas_sparse_matrix> uBlasSparseMatrix;
> 
> If SWIG does not know what ublas_sparse_matrix is, it cannot
> instantiate the template class at all.
> 
> For SWIG to be aware of uBlas, SWIG needs to wrap uBlas as well. I've
> made an attempt, but the SWIG parser cannot parse the uBlas
> headers. They use quite advanced C++ constructs and are several
> thousand lines long. One example SWIG gets stuck at:
> 
> template<class T>
> const typename zero_matrix<T>::value_type zero_matrix<T>::zero_ (0);
> 
> I cannot even manually parse this, which makes understanding what SWIG
> gets stuck at very difficult.
> 
> My suggestion is to enforce a stricter separation between interface
> and implementation for the linear algebra classes. It would be good if
> the essential public interface (i.e. the class definition) does not
> depend on the underlying implementation. I think the PETSc wrappers
> achieve this quite well:
> 
> class PETScSparseMatrix : public GenericMatrix, public Variable
> {
> ...
>     /// Return PETSc Mat pointer
>     Mat mat() const;
> ...
> }
> 
> It's possible to use PETScSparseMatrix without knowing anything about
> PETSc. Only when I want to access the underlying implementation
> (through mat()) do I need to know what a Mat actually is.
> 
> The uBlas wrappers could be defined in a similar way:
> 
> class uBlasSparseMatrix : public Variable,
>                           public GenericMatrix,
>                           public uBlasKrylovMatrix
> {
> ...
>     /// Return uBlas matrix pointer
>     ublas_sparse_matrix* mat() const;
> ...
> }
> 
> Or alternatively:
> 
> class uBlasMatrix : public Variable,
>                     public GenericMatrix,
>                     public uBlasKrylovMatrix
> {
> ...
>     /// Return uBlas matrix pointer
>     uBlasMatrixImplementation* mat() const;
> ...
> }
> 

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. 

Would it help if rather than having 

   typedef uBlasMatrix<ublas_sparse_matrix> uBlasSparseMatrix;

we have 

   class  uBlasSparseMatrix : uBlasMatrix<ublas_sparse_matrix> {};?

Garth

> where uBlasMatrixImplemntation can be templated and be either sparse
> or dense.
> 
> What do you think?
> 
>   Johan
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev