L.Ferron
L.Ferron

Reputation: 63

Mapping complex sparse matrix in Eigen from MATLAB workspace

I am working on solving the linear algebraic equation Ax = b by using Eigen solvers through mex function of Matlab. Given a complex sparse matrix A and a sparse vector b from Matlab workspace, I want to map matrix A and vector b in Eigen sparse matrix format. After that, I need to use Eigen's linear equation solvers to solve it. At the end I need to transfer the results x to Matlab workspace.

However, since I am not good at C++ and not familiar with Eigen either. I am stuck at the first step, namely constructing the complex sparse matrix in Eigen accepted format.

I have found there is the following function in Eigen,

Eigen::MappedSparseMatrix<double,RowMajor>  mat(rows, cols, nnz, row_ptr, col_index, values);

And I can use mxGetPr, mxGetPi, mxGetIr, mxGetJc, etc, these mex functions to get the info for the above "rows, cols, nnz, row_ptr, col_index, values". However, since in my case, matrix A is a complex sparse matrix, I am not sure whether "MappedSparseMatrix" can do that.

If it can, how the format of "MappedSparseMatrix" should be ? Is the following correct ?

Eigen::MappedSparseMatrix<std::complex<double>> mat(rows, cols, nnz, row_ptr, col_index, values_complex);

If so, how should I construct that values_complex ? I have found about a relevant topic before. I can use the following codes to get a complex dense matrix.

MatrixXcd mat(m,n);
mat.real() = Map<MatrixXd>(realData,m,n);
mat.imag() = Map<MatrixXd>(imagData,m,n);

However, since my matrix A is a sparse matrix, it seems that it will produce errors if I define mat as a complex sparse matrix like the following:

SparseMatrix<std::complex<double> > mat;
mat.real()  =  Map<SparseMatrix>(rows, cols, nnz, row_ptr, col_index, realData);
mat.imag() = Map<SparseMatrix>(rows, cols, nnz, row_ptr, col_index, imagData);

So can anyone provide some advice for that?

Upvotes: 0

Views: 552

Answers (1)

ggael
ggael

Reputation: 29205

MatlLab stores complex entries in two separate buffers: one for the real components and one for the imaginary components, whereas Eigen needs them to be interleaved:

value_ptr = [r0,i0,r1,i1,r2,i2,...]

so that it is compatible with std::complex<>. So in your case, you will have to create yourself a temporary buffer holding the values in that interleaved format to be passed to MappedSparseMatrix, or, if using Eigen 3.3, to Map<SparseMatrix<double,RowMajor> >.

Moreover, you will have to adjust the buffer of indices so that they are zero-based. To this end, decrement by one all entries of col_ptr and row_ptr before passing them to Eigen, and increment them by one afterward.

Upvotes: 0

Related Questions