Support of -DEIGEN_DEFAULT_TO_ROW_MAJOR

[armingeiser]

I have compared the result of SymGEigsShiftSolver if compiled with and without -DEIGEN_DEFAULT_TO_ROW_MAJOR.
It seems that internally some types are hard-coded to default to col-major ordering and ignore this flag.
This leads to wrong results if the Eigen storage order is defaulting to row-major.

Are there plans to support this flag?

without -DEIGEN_DEFAULT_TO_ROW_MAJOR

Generalized eigenvalues found:
  0.185391
-0.0214049
 -0.128334

with -DEIGEN_DEFAULT_TO_ROW_MAJOR

Generalized eigenvalues found:
 5.88503e-55
-3.12854e-55
-1.31134e-36

Example code

The test example i have used is from the documentation of the SymGEigsShiftSolver class:

#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Spectra/SymGEigsShiftSolver.h>
#include <Spectra/MatOp/SymShiftInvert.h>
#include <Spectra/MatOp/SparseSymMatProd.h>
#include <iostream>

using namespace Spectra;

int main()
{
    // We are going to solve the generalized eigenvalue problem
    //     A * x = lambda * B * x,
    // where A is symmetric and B is positive definite
    const int n = 100;

    // Define the A matrix
    Eigen::MatrixXd M = Eigen::MatrixXd::Random(n, n);
    Eigen::MatrixXd A = M + M.transpose();

    // Define the B matrix, a tridiagonal matrix with 2 on the diagonal
    // and 1 on the subdiagonals
    Eigen::SparseMatrix<double> B(n, n);
    B.reserve(Eigen::VectorXi::Constant(n, 3));
    for (int i = 0; i < n; i++)
    {
        B.insert(i, i) = 2.0;
        if (i > 0)
            B.insert(i - 1, i) = 1.0;
        if (i < n - 1)
            B.insert(i + 1, i) = 1.0;
    }

    // Construct matrix operation objects using the wrapper classes
    // A is dense, B is sparse
    using OpType = SymShiftInvert<double, Eigen::Dense, Eigen::Sparse>;
    using BOpType = SparseSymMatProd<double>;
    OpType op(A, B);
    BOpType Bop(B);

    // Construct generalized eigen solver object, seeking three generalized
    // eigenvalues that are closest to zero. This is equivalent to specifying
    // a shift sigma = 0.0 combined with the SortRule::LargestMagn selection rule
    SymGEigsShiftSolver<OpType, BOpType, GEigsMode::ShiftInvert>
        geigs(op, Bop, 3, 6, 0.0);

    // Initialize and compute
    geigs.init();
    int nconv = geigs.compute(SortRule::LargestMagn);

    // Retrieve results
    Eigen::VectorXd evalues;
    Eigen::MatrixXd evecs;
    if (geigs.info() == CompInfo::Successful)
    {
        evalues = geigs.eigenvalues();
        evecs = geigs.eigenvectors();
    }

    std::cout << "Number of converged generalized eigenvalues: " << nconv << std::endl;
    std::cout << "Generalized eigenvalues found:\n" << evalues << std::endl;
    std::cout << "Generalized eigenvectors found:\n" << evecs.topRows(10) << std::endl;

    return 0;
}

[armingeiser]

I tried to replace all occurences of Eigen::ColMajor in Spectra with EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION (which is Eigen::ColMajor/Eigen::RowMajor depending on the -DEIGEN_DEFAULT_TO_ROW_MAJOR flag).

Then i get a segfault however.

[yixuan]

It seems that DEIGEN_DEFAULT_TO_ROW_MAJOR is only for testing and debugging (https://eigen.tuxfamily.org/dox/TopicPreprocessorDirectives.html), and the Eigen documentation does not recommend using it in production code.

If the input matrix is row-majored, then there are flags in SymShiftInvert and SparseSymMatProd to support the layout.

[armingeiser]

If the input matrix is row-majored, then there are flags in SymShiftInvert and SparseSymMatProd to support the layout.

Thanks! This works, for me and i have used it also in #36. (compiled without the mentioned flag)

Along the way i have noted that i got wrong results with -DEIGEN_DEFAULT_TO_ROW_MAJOR and was just curious if this flag was/is intended do be supported.

Maybe a warning in the docs of spectra that it should not be used in combination with this Eigen flag (since it is giving wrong results, something that i did not observe yet for Eigen alone) could prevent future users from this pitfall.

[yixuan]

Yeah, internally Spectra relies on the column-majored layout for efficient operations (e.g. copying columns), and I was even unaware of the existence of the EIGEN_DEFAULT_TO_ROW_MAJOR macro, so at least in the current stage Spectra does not support this flag.

I once hoped that when the matrix layout is different from the wrapper class (e.g. SparseSymMatProd), a compiler error can occur. But it turns out that Eigen is so "robust" that such code still compiles. I'll add more documentation to emphasize this.

By the way, I noticed that you had a PR in KratosMultiphysics/Kratos#8440. Do you need to support the EIGEN_DEFAULT_TO_ROW_MAJOR flag globally? It may be easier if only the input matrix is declared row-majored, and no internal code is affected.

[armingeiser]

The KratosMultiphysics project uses row-majored ublas matrices internally. We are interfacing Eigen for access to some linear algebra functionality and for convenience we have used the -DEIGEN_DEFAULT_TO_ROW_MAJOR flag globally until now, such that everything has the same ordering by default.

In PR 8440 i have interfaced spectra on top of that, and in this process i have detected the issues with the -DEIGEN_DEFAULT_TO_ROW_MAJOR flag. Thats why i have switched to explicitly define only the input matrices to be row-major and have removed the -DEIGEN_DEFAULT_TO_ROW_MAJOR flag.
The PR is still in review, so lets see if it goes in. At the moment it seems that this solution will work.