Exporting a matrix from FEniCS 1.6.0 / Problems with EigenMatrix

Question

My problem is with exporting sparse matrices created in FEniCS, (e.g., to SciPy or Matlab).
I know that there were Q&As on the issue, but it seems that the solutions I have found:
http://fenicsproject.org/qa/8400/csr-representation-of-ublas-matrix
http://fenicsproject.org/qa/2257/how-to-read-data-from-class-matrix
http://fenicsproject.org/qa/3221/how-to-import-matrix-in-matlab-to-fenics
no longer apply in FEniCS 1.6.0.

I also know that the uBLAS linear algebra backend is no longer supported and it is advised to use the Eigen backend instead:
http://fenicsproject.org/qa/8166/query-on-parameters-linear_algebra_backend-ublas
However, it seems that the Eigen backend is not working properly. In the below code, I'm trying to assemble a matrix using the Eigen backend and then access its data:

from dolfin import *

parameters['linear_algebra_backend'] = 'Eigen'

N    = 4
mesh = UnitSquareMesh(N,N)
V    = FunctionSpace(mesh,"Lagrange",1)
u, v = TrialFunction(V), TestFunction(V)

MassMat = assemble(u*v*dx)

print MassMat.data()

getting an error:

Traceback (most recent call last):
  File "error.py", line 15, in <module>
    print MassMat.data()
AttributeError: 'Matrix' object has no attribute 'data'

It seems that MassMat is not recognized as an EigenMatrix class instance.

However, the class EigenMatrix itself works correctly. Running the below gives no errors:

from dolfin import *

M = EigenMatrix(3,3)

print M.data()

But what I want is to play with the MassMat created with the assemble function. Is there anything I'm doing wrong or is it a bug? If a bug, is there any other method to export a sparse matrix created in FEniCS 1.6.0? If not, the only option seems to be downgrading to older FEniCS releases, with uBLAS support, which is not a solution I dream of.

=====
P.S. There is also an issue concerning the sparray() method which I don't understand. The below code:

from dolfin import *

M  = EigenMatrix(3,3)

print M.sparray()

gives an error:

Traceback (most recent call last):
  File "error.py", line 8, in <module>
    print M.sparray()
  File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 2245, in sparray
    C = csr_matrix((data[2], data[1], data[0]))
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 77, in __init__
    raise ValueError('unable to infer matrix dimensions')
ValueError: unable to infer matrix dimensions

is there something more I need to do with the M.sparray() data to handle it properly?

Answer 1

Hi,

To get an EigenMatrix from a Matrix try this (not tested)

M2 = M.as_backend_type()

Another option is this code which I got from MiroK:

def mat_to_csr(dolfin_matrix):
    """
    Convert any dolfin.Matrix to csr matrix in scipy.
    Based on code by Miroslav Kuchta
    """
    assert df.MPI.size(df.mpi_comm_world()) == 1, 'mat_to_csr assumes single process'

    rows = [0]
    cols = []
    values = []
    for irow in range(dolfin_matrix.size(0)):
        indices, values_ = dolfin_matrix.getrow(irow)
        rows.append(len(indices)+rows[-1])
        cols.extend(indices)
        values.extend(values_)

    shape = dolfin_matrix.size(0), dolfin_matrix.size(1)

    return scipy.sparse.csr_matrix((numpy.array(values, dtype='float'),
                                    numpy.array(cols, dtype='int'),
                                    numpy.array(rows, dtype='int')),
                                    shape)

Comments

Sorry for a year too late response, I had decided to downgrade to 1.5.0 before you posted your answer. Now that I have finally decided to upgrade and installed the latest release (2017.1.0), I can test your solution.

Yes, as_backend_type solved my problem. However, it seems that it should be called in a function form, not as a method (as your answer suggests). At least in 2017.1.0, the following does not work

MassMat.as_backend_type()

Instead, it is legal to do

as_backend_type(MassMat)

To sum up, assuming the use of the Eigen backend and 2017.1.0 release, the CSR data of the matrix MassMat can be extracted with

rows,cols,vals = as_backend_type(MassMat).data()

Answer 2

You can't get sparrray() from an uniniitialised EigenMatrix (maybe there should be better error checks here).

What you need is

as_backend_type(MassMat).sparray()

Comments

I didn't response a year ago because prior to your answer I decided to downgraded to 1.5.0 and stick to the uBLAS backend. Please excuse me. Now I have upgraded to 2017.1.0 and become forced to abandon uBLAS, so here I write again.

Yes, as_backend_type is what I missed. Both of these two lines works as expected:

rows,cols,vals = as_backend_type(MassMat).data() # Get CSR data
M = as_backend_type(MassMat).sparray() # Get SciPy sparse matrix

The only thing I wonder about is as follows. One can directly get the SciPy sparse matrix from a dolfin EigenMatrix using the sparray() method. So why it is so frequently advised to get the SciPy sparse matrix by using the data() method which subsequently involves an additional operation of manually creating a SciPy sparse matrix from the data? Using data() is the most frequent solution in the links I provided in my original post and appears even in the FEniCS book (Chap. 10.4, p. 223).