SymGEigsSolver.h

// Copyright (C) 2016-2025 Yixuan Qiu <yixuan.qiu@cos.name>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
#ifndef SPECTRA_SYM_GEIGS_SOLVER_H
#define SPECTRA_SYM_GEIGS_SOLVER_H
 
#include "HermEigsBase.h"
#include "Util/GEigsMode.h"
#include "MatOp/internal/SymGEigsCholeskyOp.h"
#include "MatOp/internal/SymGEigsRegInvOp.h"
 
namespace Spectra {


 
// Empty class template
template <typename OpType, typename BOpType, GEigsMode Mode>
class SymGEigsSolver
{};

 
// Partial specialization for mode = GEigsMode::Cholesky
template <typename OpType, typename BOpType>
class SymGEigsSolver<OpType, BOpType, GEigsMode::Cholesky> :
    public HermEigsBase<SymGEigsCholeskyOp<OpType, BOpType>, IdentityBOp>
{
private:
    using Scalar = typename OpType::Scalar;
    using Index = Eigen::Index;
    using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
    using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
 
    using ModeMatOp = SymGEigsCholeskyOp<OpType, BOpType>;
    using Base = HermEigsBase<ModeMatOp, IdentityBOp>;
 
    const BOpType& m_Bop;
 
public:
    SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv) :
        Base(ModeMatOp(op, Bop), IdentityBOp(), nev, ncv),
        m_Bop(Bop)
    {}

 
    Matrix eigenvectors(Index nvec) const override
    {
        Matrix res = Base::eigenvectors(nvec);
        Vector tmp(res.rows());
        const Index nconv = res.cols();
        for (Index i = 0; i < nconv; i++)
        {
            m_Bop.upper_triangular_solve(&res(0, i), tmp.data());
            res.col(i).noalias() = tmp;
        }
 
        return res;
    }
 
    Matrix eigenvectors() const override
    {
        return SymGEigsSolver<OpType, BOpType, GEigsMode::Cholesky>::eigenvectors(this->m_nev);
    }

};

 
// Partial specialization for mode = GEigsMode::RegularInverse
template <typename OpType, typename BOpType>
class SymGEigsSolver<OpType, BOpType, GEigsMode::RegularInverse> :
    public HermEigsBase<SymGEigsRegInvOp<OpType, BOpType>, BOpType>
{
private:
    using Index = Eigen::Index;
 
    using ModeMatOp = SymGEigsRegInvOp<OpType, BOpType>;
    using Base = HermEigsBase<ModeMatOp, BOpType>;
 
public:
    SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv) :
        Base(ModeMatOp(op, Bop), Bop, nev, ncv)
    {}
};
 
}  // namespace Spectra
 
#endif  // SPECTRA_SYM_GEIGS_SOLVER_H