SymGEigsShiftSolver.h

// Copyright (C) 2020-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_GEIGS_SHIFT_SOLVER_H
#define SPECTRA_SYM_GEIGS_SHIFT_SOLVER_H
 
#include <utility>  // std::move
#include "SymEigsBase.h"
#include "Util/GEigsMode.h"
#include "MatOp/internal/SymGEigsShiftInvertOp.h"
#include "MatOp/internal/SymGEigsBucklingOp.h"
#include "MatOp/internal/SymGEigsCayleyOp.h"
 
namespace Spectra {
 
 
// Empty class template
template <typename OpType, typename BOpType, GEigsMode Mode>
class SymGEigsShiftSolver
{};
 
 
// Partial specialization for mode = GEigsMode::ShiftInvert
template <typename OpType, typename BOpType>
class SymGEigsShiftSolver<OpType, BOpType, GEigsMode::ShiftInvert> :
    public SymEigsBase<SymGEigsShiftInvertOp<OpType, BOpType>, BOpType>
{
private:
    using Scalar = typename OpType::Scalar;
    using Index = Eigen::Index;
    using Array = Eigen::Array<Scalar, Eigen::Dynamic, 1>;
 
    using ModeMatOp = SymGEigsShiftInvertOp<OpType, BOpType>;
    using Base = SymEigsBase<ModeMatOp, BOpType>;
    using Base::m_nev;
    using Base::m_ritz_val;
 
    const Scalar m_sigma;
 
    // Set shift and forward
    static ModeMatOp set_shift_and_move(ModeMatOp&& op, const Scalar& sigma)
    {
        op.set_shift(sigma);
        return std::move(op);
    }
 
    // 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:
    SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) :
        Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv),
        m_sigma(sigma)
    {}
};
 
 
// Partial specialization for mode = GEigsMode::Buckling
template <typename OpType, typename BOpType>
class SymGEigsShiftSolver<OpType, BOpType, GEigsMode::Buckling> :
    public SymEigsBase<SymGEigsBucklingOp<OpType, BOpType>, BOpType>
{
private:
    using Scalar = typename OpType::Scalar;
    using Index = Eigen::Index;
    using Array = Eigen::Array<Scalar, Eigen::Dynamic, 1>;
 
    using ModeMatOp = SymGEigsBucklingOp<OpType, BOpType>;
    using Base = SymEigsBase<ModeMatOp, BOpType>;
    using Base::m_nev;
    using Base::m_ritz_val;
 
    const Scalar m_sigma;
 
    // Set shift and forward
    static ModeMatOp set_shift_and_move(ModeMatOp&& op, const Scalar& sigma)
    {
        if (sigma == Scalar(0))
            throw std::invalid_argument("SymGEigsShiftSolver: sigma cannot be zero in the buckling mode");
        op.set_shift(sigma);
        return std::move(op);
    }
 
    // 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 = lambda / (lambda - sigma)
        // So the eigenvalues of the original problem is lambda = sigma * nu / (nu - 1)
        m_ritz_val.head(m_nev).array() = m_sigma * m_ritz_val.head(m_nev).array() /
            (m_ritz_val.head(m_nev).array() - Scalar(1));
        Base::sort_ritzpair(sort_rule);
    }
 
public:
    SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) :
        Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv),
        m_sigma(sigma)
    {}
};
 
 
// Partial specialization for mode = GEigsMode::Cayley
template <typename OpType, typename BOpType>
class SymGEigsShiftSolver<OpType, BOpType, GEigsMode::Cayley> :
    public SymEigsBase<SymGEigsCayleyOp<OpType, BOpType>, BOpType>
{
private:
    using Scalar = typename OpType::Scalar;
    using Index = Eigen::Index;
    using Array = Eigen::Array<Scalar, Eigen::Dynamic, 1>;
 
    using ModeMatOp = SymGEigsCayleyOp<OpType, BOpType>;
    using Base = SymEigsBase<ModeMatOp, BOpType>;
    using Base::m_nev;
    using Base::m_ritz_val;
 
    const Scalar m_sigma;
 
    // Set shift and forward
    static ModeMatOp set_shift_and_move(ModeMatOp&& op, const Scalar& sigma)
    {
        if (sigma == Scalar(0))
            throw std::invalid_argument("SymGEigsShiftSolver: sigma cannot be zero in the Cayley mode");
        op.set_shift(sigma);
        return std::move(op);
    }
 
    // 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 = (lambda + sigma) / (lambda - sigma)
        // So the eigenvalues of the original problem is lambda = sigma * (nu + 1) / (nu - 1)
        m_ritz_val.head(m_nev).array() = m_sigma * (m_ritz_val.head(m_nev).array() + Scalar(1)) /
            (m_ritz_val.head(m_nev).array() - Scalar(1));
        Base::sort_ritzpair(sort_rule);
    }
 
public:
    SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) :
        Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv),
        m_sigma(sigma)
    {}
};
 
}  // namespace Spectra
 
#endif  // SPECTRA_SYM_GEIGS_SHIFT_SOLVER_H