GenEigsComplexShiftSolver.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_GEN_EIGS_COMPLEX_SHIFT_SOLVER_H
#define SPECTRA_GEN_EIGS_COMPLEX_SHIFT_SOLVER_H
 
#include <Eigen/Core>
 
#include "GenEigsBase.h"
#include "Util/SelectionRule.h"
#include "MatOp/DenseGenComplexShiftSolve.h"
 
namespace Spectra {
 
template <typename OpType = DenseGenComplexShiftSolve<double>>
class GenEigsComplexShiftSolver : public GenEigsBase<OpType, IdentityBOp>
{
private:
    using Scalar = typename OpType::Scalar;
    using Index = Eigen::Index;
    using Complex = std::complex<Scalar>;
    using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
    using ComplexArray = Eigen::Array<Complex, Eigen::Dynamic, 1>;
 
    using Base = GenEigsBase<OpType, IdentityBOp>;
    using Base::m_op;
    using Base::m_n;
    using Base::m_nev;
    using Base::m_fac;
    using Base::m_ritz_val;
    using Base::m_ritz_vec;
 
    const Scalar m_sigmar;
    const Scalar m_sigmai;
 
    // First transform back the Ritz values, and then sort
    void sort_ritzpair(SortRule sort_rule) override
    {
        using std::abs;
        using std::sqrt;
        using std::norm;
 
        // The eigenvalues we get from the iteration is
        //     nu = 0.5 * (1 / (lambda - sigma) + 1 / (lambda - conj(sigma)))
        // So the eigenvalues of the original problem is
        //                       1 \pm sqrt(1 - 4 * nu^2 * sigmai^2)
        //     lambda = sigmar + -----------------------------------
        //                                     2 * nu
        // We need to pick the correct root
        // Let (lambdaj, vj) be the j-th eigen pair, then A * vj = lambdaj * vj
        // and inv(A - r * I) * vj = 1 / (lambdaj - r) * vj
        // where r is any shift value.
        // We can use this identity to determine lambdaj
        //
        // op(v) computes Re(inv(A - r * I) * v) for any real v
        // If r is real, then op(v) is also real. Let a = Re(vj), b = Im(vj),
        // then op(vj) = op(a) + op(b) * i
        // By comparing op(vj) and [1 / (lambdaj - r) * vj], we can determine
        // which one is the correct root
 
        // Select a random shift value
        SimpleRandom<Scalar> rng(0);
        const Scalar shiftr = rng.random() * m_sigmar + rng.random();
        const Complex shift = Complex(shiftr, Scalar(0));
        m_op.set_shift(shiftr, Scalar(0));
 
        // Calculate inv(A - r * I) * vj
        Vector v_real(m_n), v_imag(m_n), OPv_real(m_n), OPv_imag(m_n);
        constexpr Scalar eps = TypeTraits<Scalar>::epsilon();
        for (Index i = 0; i < m_nev; i++)
        {
            v_real.noalias() = m_fac.matrix_V() * m_ritz_vec.col(i).real();
            v_imag.noalias() = m_fac.matrix_V() * m_ritz_vec.col(i).imag();
            m_op.perform_op(v_real.data(), OPv_real.data());
            m_op.perform_op(v_imag.data(), OPv_imag.data());
 
            // Two roots computed from the quadratic equation
            const Complex nu = m_ritz_val[i];
            const Complex root_part1 = m_sigmar + Scalar(0.5) / nu;
            const Complex root_part2 = Scalar(0.5) * sqrt(Scalar(1) - Scalar(4) * m_sigmai * m_sigmai * (nu * nu)) / nu;
            const Complex root1 = root_part1 + root_part2;
            const Complex root2 = root_part1 - root_part2;
 
            // Test roots
            Scalar err1 = Scalar(0), err2 = Scalar(0);
            for (int k = 0; k < m_n; k++)
            {
                const Complex rhs1 = Complex(v_real[k], v_imag[k]) / (root1 - shift);
                const Complex rhs2 = Complex(v_real[k], v_imag[k]) / (root2 - shift);
                const Complex OPv = Complex(OPv_real[k], OPv_imag[k]);
                err1 += norm(OPv - rhs1);
                err2 += norm(OPv - rhs2);
            }
 
            const Complex lambdaj = (err1 < err2) ? root1 : root2;
            m_ritz_val[i] = lambdaj;
 
            if (abs(Eigen::numext::imag(lambdaj)) > eps)
            {
                m_ritz_val[i + 1] = Eigen::numext::conj(lambdaj);
                i++;
            }
            else
            {
                m_ritz_val[i] = Complex(Eigen::numext::real(lambdaj), Scalar(0));
            }
        }
 
        Base::sort_ritzpair(sort_rule);
    }
 
public:
    GenEigsComplexShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai) :
        Base(op, IdentityBOp(), nev, ncv),
        m_sigmar(sigmar), m_sigmai(sigmai)
    {
        op.set_shift(m_sigmar, m_sigmai);
    }
};
 
}  // namespace Spectra
 
#endif  // SPECTRA_GEN_EIGS_COMPLEX_SHIFT_SOLVER_H