Reputation: 63
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
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