I'm trying to use the Spectra 3.5 Library on my Linux machine, and the SparseGenMatProd wrapper for the Matrix-Vector multiplication only seems to work when the sparse matrix is in ColMajor format. Is this normal behavior and if so, how can I fix it to take RowMajor format? I've included a basic example where the output is "Segmentation fault (core dumped)". I've gone through several other posts and the documentation, but can't seem to find an answer.
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <GenEigsSolver.h>
#include <MatOp/SparseGenMatProd.h>
#include <iostream>
using namespace Spectra;
int main()
{
// A band matrix with 1 on the main diagonal, 2 on the below-main subdiagonal,
// and 3 on the above-main subdiagonal
const int n = 10;
Eigen::SparseMatrix<double, Eigen::RowMajor> M(n, n);
M.reserve(Eigen::VectorXi::Constant(n, 3));
for(int i = 0; i < n; i++)
{
M.insert(i, i) = 1.0;
if(i > 0)
M.insert(i - 1, i) = 3.0;
if(i < n - 1)
M.insert(i + 1, i) = 2.0;
}
// Construct matrix operation object using the wrapper class SparseGenMatProd
SparseGenMatProd<double> op(M);
// Construct eigen solver object, requesting the largest three eigenvalues
GenEigsSolver< double, LARGEST_MAGN, SparseGenMatProd<double> > eigs(&op, 3, 6);
// Initialize and compute
eigs.init();
int nconv = eigs.compute();
// Retrieve results
Eigen::VectorXcd evalues;
if(eigs.info() == SUCCESSFUL)
evalues = eigs.eigenvalues();
std::cout << *emphasized text*"Eigenvalues found:\n" << evalues << std::endl;
return 0;
}
If you change line 15 to:
Eigen::SparseMatrix<double, Eigen::ColMajor> M(n, n);
it will work as expected.
Currently I'm working around this and converting my matrices to ColMajor, but I'd like to understand what's going on. Any help is much appreciated.
The API of SparseGenMatProd
seems to be quite misleading. Looks like you have to specify you are dealing with row-major matrices through the second template parameter:
SparseGenMatProd<double,RowMajor> op(M);
otherwise M
is implicitly converted to a temporary column-major matrix which is then stored by const reference by op
but this temporary is dead right after that line of code.