Trouble Creating Custom SymShiftInvert operation that uses Eigen::UmfpackLU instead of Eigen::SparseLU

[otmanon]

Hi! I'm trying to solve a generalized eigenvalue problem on huge 100kx100k sparse matrices (that are not generally positive definite). I'm using the default shift invert solver. The SparseLU solver that SymShiftInvert uses seems to take forever.

To this end, it seems Matlab uses Umfpack and seeing as Eigen has a nice Wrapper for it, I wanted to see if it'd be as simple as replacing the SparseLU solver in my custom SymShiftInvert class.

Unfortunately it seems the GEigsSymShiftInvert solver uses move constructors that seem to destroy/overwrite crucial memory in the UmfpackLU solver. Specifically the mp_matrix member variable of the UmfpackLU solver becomes out of scope once we pass it to our generalized eigen solver, causing my programme to crash.

Is there a method you recommend using UmfpackLU instead of the default SparseLU for my SymShiftInvertOp? Could you shed some light as to why this overwriting is occurring?
Thank you

[otmanon]

I've also confirmed separately that the UmfpackLU solver works as expected on the same system matrix I'm trying to solve. So it really is when I give it to my SymGEigsShiftSolver that it gets overwritten and doesn't behave as expected

[otmanon]

Now that I have more time I'll share some of my code:
Here I initialize my custom UMFPACKSymShiftInvert operation, as well as the default SparseMatProd for B defined by Spectra.

    using OpType = UMFPACKSymShiftInvert;
    using BopType = Spectra::SparseSymMatProd<double>;
    Eigen::UmfPackLU<Eigen::SparseMatrix<double>> solver = Eigen::UmfPackLU<Eigen::SparseMatrix<double>>(A);
    OpType op(A, B, solver);
    BopType Bop(B);

Looking into the data pointed to by solver, everything is accessible, and I've even checked the I can use the solver to solve a system like Ax=0 . So it seems all good for now.

Next, we pass our ops into the SymGeigsShiftSolver

SymGEigsShiftSolver<OpType, BopType, GEigsMode::ShiftInvert> geigs(op, Bop, r, r*2, 0);

as soon as this line executes, all the data pointed to by our solver is destroyed! I can no longer solved the same system Ax=0 as before with my solver!

Specifically, doing so leads me to a Memory out of Range error in UmfpackLU.solve(), as that method seems to make use of member variables stored in the UmfpackLU object that have now been overwritten. Specifically the mp_matrix member variable defined in line 556 of UMFPACK LU been destroyed.

For completeness, here is my custom UMFPACKSymShiftInvert class

class UMFPACKSymShiftInvert
{
public: 
    using Scalar = double;
private:
    SparseMatrix<double> m_matA;
    SparseMatrix<double> m_matB;
    const int m_n;
    Eigen::UmfPackLU<Eigen::SparseMatrix<double>>* m_solver;
public:
    UMFPACKSymShiftInvert(const SparseMatrix<double>& A, const SparseMatrix<double>& B, Eigen::UmfPackLU<Eigen::SparseMatrix<double>>* solver) :
        m_matA(A), m_matB(B), m_n(A.rows()), m_solver(solver)
    {
        if (m_n != A.cols() || m_n != B.rows() || m_n != B.cols())
            throw std::invalid_argument("UMFPACKSymShiftInvert: A and B must be square matrices of the same size");
    }
    int rows() const { return m_n; }
    int cols() const { return m_n; }
    void set_shift(const double& sigma)
    {
        Eigen::SparseMatrix<double> C = m_matA + sigma * m_matB;
        m_solver->compute(C);
        const bool success = (m_solver->info() == Eigen::Success);
        if (!success)
            throw std::invalid_argument("UMFPACKSymShiftInvert: factorization failed with the given shift");
    }
    void perform_op(const double* x_in,double* y_out) const
    {
        Eigen::VectorXd x = Eigen::Map<const Eigen::VectorXd>(x_in, m_n);
        Eigen::VectorXd y = Eigen::Map<Eigen::VectorXd>(y_out, m_n);
        y = m_solver->solve(x);
    }
};

I've played around with the format a bunch, such as storing my solver pointer on a heap, vs not. Making my solver a reference rather than a pointer to the solver defined outside... even creating and storing the solver entirely in the UMFPACKSymShiftInvert class (rather than a pointer). This same problem happens with anything I try.

Is there any bandaid I can put on this that will let me use UMFPACKSLU?

[otmanon]

I did a little more investigating and came up with an inelegant bandaid fix.

 OpType op(A, B, solver);
 BopType Bop(B);

It seems that the data stored in either operator's memory gets cleared after being passed through

SymGEigsShiftSolver<OpType, BopType, GEigsMode::ShiftInvert> geigs(op, Bop, r, r*2, 0);

I'll admit I'm not a C++ expert but I've noticed things like std::move that might clear some of the source data.

My temporary fix is to simply re-initalize the relevant member variables in both my custom op and Bop by calling new setters like so:

    //initialize operators
    using OpType = UMFPACKSymShiftInvert;
    using BopType = MatProd;
    Eigen::UmfPackLU<Eigen::SparseMatrix<double>>* solver = new Eigen::UmfPackLU<Eigen::SparseMatrix<double>>(A);
    OpType op(A, B, solver);
    BopType Bop(B);

     //this line corrupts op and Bop
    SymGEigsShiftSolver<OpType, BopType, GEigsMode::ShiftInvert> geigs(op, Bop, r, r*2, 0);
    //reinitialize invert operator
    solver = new Eigen::UmfPackLU<Eigen::SparseMatrix<double>>(A);
    op.set_mat(A, B);
    op.set_solver(solver);
    op.set_shift(0);
    //reinitialize matrix product operator
    Bop.set_mat(B);
    geigs.init();

And the code for my custom OpType UMFPACKSymShiftInvert

class UMFPACKSymShiftInvert
{
public: 
    using Scalar = double;
private:
    SparseMatrix<double> m_matA;
    SparseMatrix<double> m_matB;
    Eigen::SparseMatrix<double> C;
    
    const int m_n;
    Eigen::UmfPackLU<Eigen::SparseMatrix<double>>* m_solver;

public:
    ///
    /// Constructor to create the matrix operation object.
    ///
    UMFPACKSymShiftInvert(const SparseMatrix<double>& A, const SparseMatrix<double>& B, Eigen::UmfPackLU<Eigen::SparseMatrix<double>>* solver) :
        m_matA(A), m_matB(B), m_n(A.rows()), m_solver(solver)
    {
        if (m_n != A.cols() || m_n != B.rows() || m_n != B.cols())
            throw std::invalid_argument("UMFPACKSymShiftInvert: A and B must be square matrices of the same size");
    }
    int rows() const { return m_n; }
    int cols() const { return m_n; }
    void set_mat(Eigen::SparseMatrix<double>& A, Eigen::SparseMatrix<double> & B)
    {
        m_matA = A;
        m_matB = B;
    }
    void set_solver(Eigen::UmfPackLU<Eigen::SparseMatrix<double>>* solver)
    {
        m_solver = solver;
    }
    void set_shift(const double& sigma)
    {
        C = m_matA + sigma * m_matB;
        m_solver->compute(C);
        const bool success = (m_solver->info() == Eigen::Success);
        if (!success)
            throw std::invalid_argument("UMFPACKSymShiftInvert: factorization failed with the given shift");
    }
    void perform_op(const double* x_in,double* y_out) const
    {
        Eigen::Map<const Eigen::VectorXd>x(x_in, m_n);
        Eigen::Map<Eigen::VectorXd>y(y_out, m_n);
        y = m_solver->solve(x);
    }
};

And finally my new custom MatProd (Which I'm surprised was also being overwritten..?)

class MatProd
{
public:
    using Scalar = double;
private:

    SparseMatrix<double> m_matB;
    const int m_n;
public:
    ///
    /// Constructor to create the matrix operation object.
    ///
    MatProd(const SparseMatrix<double>& B) :
        m_matB(B), m_n(B.rows())
    {
        if (m_n != B.rows())
            throw std::invalid_argument("UMFPACKSymShiftInvert: A and B must be square matrices of the same size");
    }
    int rows() const { return m_n; }
    int cols() const { return m_n; }
    void set_mat(Eigen::SparseMatrix<double>& B)
    {
        m_matB = B;
    }
    void perform_op(const double* x_in, double* y_out) const
    {
        Eigen::Map<const Eigen::VectorXd>x(x_in, m_n);
        Eigen::Map<Eigen::VectorXd> y(y_out, m_n);
        y = m_matB*(x);
    }
};

So yeah! This was my bandaid fix to making UmfPackLU work for me with Spectra! I'm wondering if you know of a more elegant /permanent solution, or if you can shed any light as to what design choices are incompatible with each other! Maybe I shouldn't have spent so much time on it but I'm glad I did because my Eigen decompositions go wayy faster now !