dolfin team mailing list archive

[Johan Jansson <johanjan@xxxxxxxxxxxxxxxx>]

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

where uBlasMatrixImplemntation can be templated and be either sparse
or dense.

What do you think?

  Johan