The core dumped can happen when using row major Eigen::SparseMatrix

[qzxuhui]

The core dumped can happen when using row major EigenSparseMatrix.

I am not sure it is a my code bug or spectra's bug.

The follow code can used to reproduce problem.

#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Eigen/Eigenvalues>
#include <Spectra/SymGEigsSolver.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
#include <Spectra/MatOp/SparseCholesky.h>
#include <iostream>
int main()
{
    // We are going to solve the generalized eigenvalue problem A * x = lambda * B * x
    const int n = 10;

    // Define the A matrix, a band matrix with 3 on the diagonal and 1 on the subdiagonals
    // Define the B matrix, a band matrix with 2 on the diagonal and 1 on the subdiagonals
    Eigen::SparseMatrix<double, Eigen::RowMajor> A(n, n);
    Eigen::SparseMatrix<double, Eigen::RowMajor> B(n, n);
    std::vector<Eigen::Triplet<double>> coefficients_A;
    std::vector<Eigen::Triplet<double>> coefficients_B;
    for (int i = 0; i < n; i++)
    {
        coefficients_A.push_back(Eigen::Triplet<double>(i, i, 3));
        coefficients_B.push_back(Eigen::Triplet<double>(i, i, 2));
        if (i > 0)
        {
            coefficients_A.push_back(Eigen::Triplet<double>(i - 1, i, 1));
            coefficients_B.push_back(Eigen::Triplet<double>(i - 1, i, 1));
        }

        if (i < n - 1)
        {
            coefficients_A.push_back(Eigen::Triplet<double>(i + 1, i, 1));
            coefficients_B.push_back(Eigen::Triplet<double>(i + 1, i, 1));
        }
    }
    A.setFromTriplets(coefficients_A.begin(), coefficients_A.end());
    B.setFromTriplets(coefficients_B.begin(), coefficients_B.end());

    // Construct matrix operation object using the wrapper classes
    Spectra::SparseSymMatProd<double> Aop(A);
    Spectra::SparseCholesky<double> Bop(B);

    // Construct generalized eigen solver object, requesting the smallest three generalized eigenvalues
    const int num_solution = 5;
    const int convergence_parameter = 8;

    Spectra::SymGEigsSolver<double, Spectra::SMALLEST_MAGN,
                            Spectra::SparseSymMatProd<double>,
                            Spectra::SparseCholesky<double>,
                            Spectra::GEIGS_CHOLESKY>
        eigen_solver(&Aop, &Bop, num_solution, convergence_parameter);

    // Initialize and compute
    eigen_solver.init();
    int num_convergence = eigen_solver.compute();
    std::cout << "num convergence = " << num_convergence << std::endl;

    // Retrieve results
    Eigen::VectorXd eigen_value;
    Eigen::MatrixXd eigen_vector;
    if (eigen_solver.info() == Spectra::SUCCESSFUL)
    {
        eigen_value = eigen_solver.eigenvalues();
        eigen_vector = eigen_solver.eigenvectors();
    }
    std::cout << "Generalized eigenvalues found:\n"
              << eigen_value << std::endl;
    std::cout << "Generalized eigenvectors found:\n"
              << eigen_vector << std::endl;

    // Verify results using the generalized eigen solver in Eigen
    Eigen::MatrixXd Adense = A;
    Eigen::MatrixXd Bdense = B;
    Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> standard_eigen_solver(Adense, Bdense);
    std::cout << "Generalized eigenvalues:\n"
              << standard_eigen_solver.eigenvalues() << std::endl;
    std::cout << "Generalized eigenvectors:\n"
              << standard_eigen_solver.eigenvectors() << std::endl;

    return 0;
}

The follow script can used to run the project on my computer

eigen_path=/opt/Eigen/eigen-3.3.7
spectra_path=/opt/spectra/spectra-0.8.1/include/

rm EigenSolverExample.out

g++ -I ${eigen_path} -I ${spectra_path} EigenSolverExample.cpp -o EigenSolverExample.out

./EigenSolverExample.out

One can see in code's top I defined two matrix

Eigen::SparseMatrix<double, Eigen::RowMajor> A(n, n);
Eigen::SparseMatrix<double, Eigen::RowMajor> B(n, n);

If I modified to

Eigen::SparseMatrix<double, Eigen::ColMajor> A(n, n);
Eigen::SparseMatrix<double, Eigen::ColMajor> B(n, n);

there will be no core dumped happen.

Could any one help me? Thanks for your time.

[yixuan]

Hi, try to use the following two lines:

Spectra::SparseSymMatProd<double, Eigen::Lower, Eigen::RowMajor> Aop(A);
Spectra::SparseCholesky<double, Eigen::Lower, Eigen::RowMajor> Bop(B);

[qzxuhui]

Sir, I try to use your suggestion, that is

#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Eigen/Eigenvalues>
#include <Spectra/SymGEigsSolver.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
#include <Spectra/MatOp/SparseCholesky.h>
#include <iostream>
int main()
{
    // We are going to solve the generalized eigenvalue problem A * x = lambda * B * x
    const int n = 10;

    // Define the A matrix, a band matrix with 3 on the diagonal and 1 on the subdiagonals
    // Define the B matrix, a band matrix with 2 on the diagonal and 1 on the subdiagonals
    Eigen::SparseMatrix<double, Eigen::RowMajor> A(n, n);
    Eigen::SparseMatrix<double, Eigen::RowMajor> B(n, n);
    std::vector<Eigen::Triplet<double>> coefficients_A;
    std::vector<Eigen::Triplet<double>> coefficients_B;
    for (int i = 0; i < n; i++)
    {
        coefficients_A.push_back(Eigen::Triplet<double>(i, i, 3));
        coefficients_B.push_back(Eigen::Triplet<double>(i, i, 2));
        if (i > 0)
        {
            coefficients_A.push_back(Eigen::Triplet<double>(i - 1, i, 1));
            coefficients_B.push_back(Eigen::Triplet<double>(i - 1, i, 1));
        }

        if (i < n - 1)
        {
            coefficients_A.push_back(Eigen::Triplet<double>(i + 1, i, 1));
            coefficients_B.push_back(Eigen::Triplet<double>(i + 1, i, 1));
        }
    }
    A.setFromTriplets(coefficients_A.begin(), coefficients_A.end());
    B.setFromTriplets(coefficients_B.begin(), coefficients_B.end());

    // Construct matrix operation object using the wrapper classes
    Spectra::SparseSymMatProd<double, Eigen::Lower, Eigen::RowMajor> Aop(A);
    Spectra::SparseCholesky<double, Eigen::Lower, Eigen::RowMajor> Bop(B);

    // Construct generalized eigen solver object, requesting the smallest three generalized eigenvalues
    const int num_solution = 5;
    const int convergence_parameter = 8;

    Spectra::SymGEigsSolver<double, Spectra::SMALLEST_MAGN,
                            Spectra::SparseSymMatProd<double>,
                            Spectra::SparseCholesky<double>,
                            Spectra::GEIGS_CHOLESKY>
        eigen_solver(&Aop, &Bop, num_solution, convergence_parameter);

    // Initialize and compute
    eigen_solver.init();
    int num_convergence = eigen_solver.compute();
    std::cout << "num convergence = " << num_convergence << std::endl;

    // Retrieve results
    Eigen::VectorXd eigen_value;
    Eigen::MatrixXd eigen_vector;
    if (eigen_solver.info() == Spectra::SUCCESSFUL)
    {
        eigen_value = eigen_solver.eigenvalues();
        eigen_vector = eigen_solver.eigenvectors();
    }
    std::cout << "Generalized eigenvalues found:\n"
              << eigen_value << std::endl;
    std::cout << "Generalized eigenvectors found:\n"
              << eigen_vector << std::endl;

    // Verify results using the generalized eigen solver in Eigen
    Eigen::MatrixXd Adense = A;
    Eigen::MatrixXd Bdense = B;
    Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> standard_eigen_solver(Adense, Bdense);
    std::cout << "Generalized eigenvalues:\n"
              << standard_eigen_solver.eigenvalues() << std::endl;
    std::cout << "Generalized eigenvectors:\n"
              << standard_eigen_solver.eigenvectors() << std::endl;

    return 0;
}

when I compile the code using script above, the compiler compained that

EigenSolverExampleRowMajor.cpp:50:69: error: no matching function for call to ‘Spectra::SymGEigsSolver<double, 4, Spectra::SparseSymMatProd<double>, Spectra::SparseCholesky<double>, 0>::SymGEigsSolver(Spectra::SparseSymMatProd<double, 1, 1>*, Spectra::SparseCholesky<double, 1, 1>*, const int&, const int&)’
         eigen_solver(&Aop, &Bop, num_solution, convergence_parameter);
                                                                     ^
In file included from EigenSolverExampleRowMajor.cpp:4:0:
/opt/spectra/spectra-0.8.1/include/Spectra/SymGEigsSolver.h:203:5: note: candidate: Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::SymGEigsSolver(OpType*, BOpType*, Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Index, Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Index) [with Scalar = double; int SelectionRule = 4; OpType = Spectra::SparseSymMatProd<double>; BOpType = Spectra::SparseCholesky<double>; Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Index = long int]
     SymGEigsSolver(OpType* op, BOpType* Bop, Index nev, Index ncv) :
     ^~~~~~~~~~~~~~
/opt/spectra/spectra-0.8.1/include/Spectra/SymGEigsSolver.h:203:5: note:   no known conversion for argument 1 from ‘Spectra::SparseSymMatProd<double, 1, 1>*’ to ‘Spectra::SparseSymMatProd<double>*’
/opt/spectra/spectra-0.8.1/include/Spectra/SymGEigsSolver.h:167:7: note: candidate: Spectra::SymGEigsSolver<double, 4, Spectra::SparseSymMatProd<double>, Spectra::SparseCholesky<double>, 0>::SymGEigsSolver(const Spectra::SymGEigsSolver<double, 4, Spectra::SparseSymMatProd<double>, Spectra::SparseCholesky<double>, 0>&)
 class SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, GEIGS_CHOLESKY>:
       ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/opt/spectra/spectra-0.8.1/include/Spectra/SymGEigsSolver.h:167:7: note:   candidate expects 1 argument, 4 provided

Thanks for your time.

[qzxuhui]

If I using

Spectra::SparseSymMatProd<double, Eigen::Lower, Eigen::RowMajor> Aop(A);
Spectra::SparseCholesky<double, Eigen::Lower, Eigen::RowMajor> Bop(B);
const int num_solution = 5;
const int convergence_parameter = 8;
Spectra::SymGEigsSolver<double, Spectra::SMALLEST_MAGN,
                            Spectra::SparseSymMatProd<double, Eigen::Lower, Eigen::RowMajor>,
                            Spectra::SparseCholesky<double, Eigen::Lower, Eigen::RowMajor>,
                            Spectra::GEIGS_CHOLESKY>
        eigen_solver(&Aop, &Bop, num_solution, convergence_parameter);

the compiler will complain

In file included from /opt/Eigen/eigen-3.3.7/Eigen/SparseCholesky:37:0,
                 from /opt/spectra/spectra-0.8.1/include/Spectra/MatOp/SparseCholesky.h:12,
                 from EigenSolverExampleRowMajor.cpp:6:
/opt/Eigen/eigen-3.3.7/Eigen/src/SparseCholesky/SimplicialCholesky.h: In instantiation of ‘static Eigen::internal::traits<Eigen::SimplicialLLT<_MatrixType, _UpLo, _Ordering> >::MatrixU Eigen::internal::traits<Eigen::SimplicialLLT<_MatrixType, _UpLo, _Ordering> >::getU(const MatrixType&) [with _MatrixType = Eigen::SparseMatrix<double, 1>; int _UpLo = 1; _Ordering = Eigen::AMDOrdering<int>; Eigen::internal::traits<Eigen::SimplicialLLT<_MatrixType, _UpLo, _Ordering> >::MatrixU = Eigen::TriangularView<const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >, 2>; typename Eigen::SparseMatrix<typename Lhs::Scalar, 0, typename Lhs::StorageIndex>::AdjointReturnType = Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >; Eigen::internal::traits<Eigen::SimplicialLLT<_MatrixType, _UpLo, _Ordering> >::MatrixType = Eigen::SparseMatrix<double, 1>]’:
/opt/Eigen/eigen-3.3.7/Eigen/src/SparseCholesky/SimplicialCholesky.h:360:28:   required from ‘const MatrixU Eigen::SimplicialLLT<_MatrixType, _UpLo, _Ordering>::matrixU() const [with _MatrixType = Eigen::SparseMatrix<double, 1>; int _UpLo = 1; _Ordering = Eigen::AMDOrdering<int>; Eigen::SimplicialLLT<_MatrixType, _UpLo, _Ordering>::MatrixU = Eigen::TriangularView<const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >, 2>]’
/opt/spectra/spectra-0.8.1/include/Spectra/MatOp/SparseCholesky.h:103:42:   required from ‘void Spectra::SparseCholesky<Scalar, Uplo, Flags, StorageIndex>::upper_triangular_solve(const Scalar*, Scalar*) const [with Scalar = double; int Uplo = 1; int Flags = 1; StorageIndex = int]’
/opt/spectra/spectra-0.8.1/include/Spectra/SymGEigsSolver.h:225:13:   required from ‘Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Matrix Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::eigenvectors(Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Index) const [with Scalar = double; int SelectionRule = 4; OpType = Spectra::SparseSymMatProd<double, 1, 1>; BOpType = Spectra::SparseCholesky<double, 1, 1>; Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Matrix = Eigen::Matrix<double, -1, -1>; Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Index = long int]’
/opt/spectra/spectra-0.8.1/include/Spectra/SymGEigsSolver.h:234:100:   required from ‘Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Matrix Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::eigenvectors() const [with Scalar = double; int SelectionRule = 4; OpType = Spectra::SparseSymMatProd<double, 1, 1>; BOpType = Spectra::SparseCholesky<double, 1, 1>; Spectra::SymGEigsSolver<Scalar, SelectionRule, OpType, BOpType, 0>::Matrix = Eigen::Matrix<double, -1, -1>]’
EigenSolverExampleRowMajor.cpp:59:50:   required from here
/opt/Eigen/eigen-3.3.7/Eigen/src/SparseCholesky/SimplicialCholesky.h:283:60: error: no matching function for call to ‘Eigen::TriangularView<const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >, 2>::TriangularView(const AdjointReturnType)’
   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
                                                            ^~~~~~~~~~~~~~~~~~~~
In file included from /opt/Eigen/eigen-3.3.7/Eigen/Core:489:0,
                 from EigenSolverExampleRowMajor.cpp:1:
/opt/Eigen/eigen-3.3.7/Eigen/src/Core/TriangularMatrix.h:217:21: note: candidate: Eigen::TriangularView<MatrixType, Mode>::TriangularView(Eigen::TriangularView<MatrixType, Mode>::MatrixType&) [with _MatrixType = const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >; unsigned int _Mode = 2; Eigen::TriangularView<MatrixType, Mode>::MatrixType = const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >]
     explicit inline TriangularView(MatrixType& matrix) : m_matrix(matrix)
                     ^~~~~~~~~~~~~~
/opt/Eigen/eigen-3.3.7/Eigen/src/Core/TriangularMatrix.h:217:21: note:   no known conversion for argument 1 from ‘const AdjointReturnType {aka const Eigen::Transpose<const Eigen::SparseMatrix<double, 1> >}’ to ‘Eigen::TriangularView<const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >, 2>::MatrixType& {aka const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >&}’
/opt/Eigen/eigen-3.3.7/Eigen/src/Core/TriangularMatrix.h:186:58: note: candidate: constexpr Eigen::TriangularView<const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >, 2>::TriangularView(const Eigen::TriangularView<const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >, 2>&)
 template<typename _MatrixType, unsigned int _Mode> class TriangularView
                                                          ^~~~~~~~~~~~~~
/opt/Eigen/eigen-3.3.7/Eigen/src/Core/TriangularMatrix.h:186:58: note:   no known conversion for argument 1 from ‘const AdjointReturnType {aka const Eigen::Transpose<const Eigen::SparseMatrix<double, 1> >}’ to ‘const Eigen::TriangularView<const Eigen::Transpose<const Eigen::SparseMatrix<double, 0, int> >, 2>&’

[yixuan]

I see. It should be an issue of the SimplicialLLT class in the Eigen library, which only supports column-majored sparse matrices. I'll submit a bug report there. Right now you have to convert B to the column-majored form.

[qzxuhui]

Thank you for your time, sir.

[dburov190]

I'd like to bump this up. Please include the information about Spectra::SparseSymMatProd<double, Eigen::Lower, Eigen::RowMajor> in your documentation. I spent three days debugging my code only to find out that it's a bug (a known one for that matter!) in your library.

Also, as a side note, it's considered bad practice to let code crash with a seg fault, no matter whether it's expected to execute properly in that branch or not. A better way is to throw exceptions.

[yixuan]

@qzxuhui The issue regarding SimplicialLLT has been fixed by the Eigen group (https://gitlab.com/libeigen/eigen/-/issues/1910), and it will appear in the next release of Eigen. You can also manually fix the Eigen source code based on that commit.

[yixuan]

@dburov190 The segment fault you encountered is not related to this one, but I will add more documentation to the SparseSymMatProd class per your suggestion. Basically I can hardly control the happening of segment fault or convert it to exceptions since all the low-level operations are done inside the Matrix class. The crash is caused by a wrong assumption on the layout of the matrix, but the compiler or runtime execution is unaware of that.

[yixuan]

The template parameters of the wrapper classes have been documented, and a static assert is added to test the consistency of the input matrix and the Flags template parameter. A compiler error will occur if a row-major matrix is passed to a SparseSymMatProd with Flags = Eigen::ColMajor.

See https://github.com/yixuan/spectra/commit/dba88dc6471cead3cf07e822fcbfded68d90a772.