1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
|
// 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_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 {
///
/// \ingroup EigenSolver
///
/// This class implements the eigen solver for general real matrices with
/// a complex shift value in the **shift-and-invert mode**. The background
/// knowledge of the shift-and-invert mode can be found in the documentation
/// of the SymEigsShiftSolver class.
///
/// \tparam OpType The name of the matrix operation class. Users could either
/// use the wrapper classes such as DenseGenComplexShiftSolve and
/// SparseGenComplexShiftSolve, or define their own that implements the type
/// definition `Scalar` and all the public member functions as in
/// DenseGenComplexShiftSolve.
///
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);
const 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:
///
/// Constructor to create a eigen solver object using the shift-and-invert mode.
///
/// \param op The matrix operation object that implements
/// the complex shift-solve operation of \f$A\f$: calculating
/// \f$\mathrm{Re}\{(A-\sigma I)^{-1}v\}\f$ for any vector \f$v\f$. Users could either
/// create the object from the wrapper class such as DenseGenComplexShiftSolve, or
/// define their own that implements all the public members
/// as in DenseGenComplexShiftSolve.
/// \param nev Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-2\f$,
/// where \f$n\f$ is the size of matrix.
/// \param ncv Parameter that controls the convergence speed of the algorithm.
/// Typically a larger `ncv` means faster convergence, but it may
/// also result in greater memory use and more matrix operations
/// in each iteration. This parameter must satisfy \f$nev+2 \le ncv \le n\f$,
/// and is advised to take \f$ncv \ge 2\cdot nev + 1\f$.
/// \param sigmar The real part of the shift.
/// \param sigmai The imaginary part of the shift.
///
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
|
