SymEigsShiftSolver.h

// Copyright (C) 2016-2022 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_EIGS_SHIFT_SOLVER_H
#define SPECTRA_SYM_EIGS_SHIFT_SOLVER_H
 
#include <Eigen/Core>
 
#include "SymEigsBase.h"
#include "Util/SelectionRule.h"
#include "MatOp/DenseSymShiftSolve.h"
 
namespace Spectra {
 
template <typename OpType = DenseSymShiftSolve<double>>
class SymEigsShiftSolver : public SymEigsBase<OpType, IdentityBOp>
{
private:
    using Scalar = typename OpType::Scalar;
    using Index = Eigen::Index;
    using Array = Eigen::Array<Scalar, Eigen::Dynamic, 1>;
 
    using Base = SymEigsBase<OpType, IdentityBOp>;
    using Base::m_nev;
    using Base::m_ritz_val;
 
    const Scalar m_sigma;
 
    // First transform back the Ritz values, and then sort
    void sort_ritzpair(SortRule sort_rule) override
    {
        // The eigenvalues we get from the iteration is nu = 1 / (lambda - sigma)
        // So the eigenvalues of the original problem is lambda = 1 / nu + sigma
        m_ritz_val.head(m_nev).array() = Scalar(1) / m_ritz_val.head(m_nev).array() + m_sigma;
        Base::sort_ritzpair(sort_rule);
    }
 
public:
    SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma) :
        Base(op, IdentityBOp(), nev, ncv),
        m_sigma(sigma)
    {
        op.set_shift(m_sigma);
    }
};
 
}  // namespace Spectra
 
#endif  // SPECTRA_SYM_EIGS_SHIFT_SOLVER_H