1 // Copyright (C) 2016-2021 Yixuan Qiu <yixuan.qiu@cos.name>
2 //
3 // This Source Code Form is subject to the terms of the Mozilla
4 // Public License v. 2.0. If a copy of the MPL was not distributed
5 // with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
6 
7 #ifndef SPECTRA_GEN_EIGS_REAL_SHIFT_SOLVER_H
8 #define SPECTRA_GEN_EIGS_REAL_SHIFT_SOLVER_H
9 
10 #include <Eigen/Core>
11 
12 #include "GenEigsBase.h"
13 #include "Util/SelectionRule.h"
14 #include "MatOp/DenseGenRealShiftSolve.h"
15 
16 namespace Spectra {
17 
18 ///
19 /// \ingroup EigenSolver
20 ///
21 /// This class implements the eigen solver for general real matrices with
22 /// a real shift value in the **shift-and-invert mode**. The background
23 /// knowledge of the shift-and-invert mode can be found in the documentation
24 /// of the SymEigsShiftSolver class.
25 ///
26 /// \tparam OpType  The name of the matrix operation class. Users could either
27 ///                 use the wrapper classes such as DenseGenRealShiftSolve and
28 ///                 SparseGenRealShiftSolve, or define their own that implements the type
29 ///                 definition `Scalar` and all the public member functions as in
30 ///                 DenseGenRealShiftSolve.
31 ///
32 template <typename OpType = DenseGenRealShiftSolve<double>>
33 class GenEigsRealShiftSolver : public GenEigsBase<OpType, IdentityBOp>
34 {
35 private:
36     using Scalar = typename OpType::Scalar;
37     using Index = Eigen::Index;
38     using Complex = std::complex<Scalar>;
39     using ComplexArray = Eigen::Array<Complex, Eigen::Dynamic, 1>;
40 
41     using Base = GenEigsBase<OpType, IdentityBOp>;
42     using Base::m_nev;
43     using Base::m_ritz_val;
44 
45     const Scalar m_sigma;
46 
47     // First transform back the Ritz values, and then sort
sort_ritzpair(SortRule sort_rule)48     void sort_ritzpair(SortRule sort_rule) override
49     {
50         // The eigenvalues we get from the iteration is nu = 1 / (lambda - sigma)
51         // So the eigenvalues of the original problem is lambda = 1 / nu + sigma
52         m_ritz_val.head(m_nev) = Scalar(1) / m_ritz_val.head(m_nev).array() + m_sigma;
53         Base::sort_ritzpair(sort_rule);
54     }
55 
56 public:
57     ///
58     /// Constructor to create a eigen solver object using the shift-and-invert mode.
59     ///
60     /// \param op     The matrix operation object that implements
61     ///               the shift-solve operation of \f$A\f$: calculating
62     ///               \f$(A-\sigma I)^{-1}v\f$ for any vector \f$v\f$. Users could either
63     ///               create the object from the wrapper class such as DenseGenRealShiftSolve, or
64     ///               define their own that implements all the public members
65     ///               as in DenseGenRealShiftSolve.
66     /// \param nev    Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-2\f$,
67     ///               where \f$n\f$ is the size of matrix.
68     /// \param ncv    Parameter that controls the convergence speed of the algorithm.
69     ///               Typically a larger `ncv` means faster convergence, but it may
70     ///               also result in greater memory use and more matrix operations
71     ///               in each iteration. This parameter must satisfy \f$nev+2 \le ncv \le n\f$,
72     ///               and is advised to take \f$ncv \ge 2\cdot nev + 1\f$.
73     /// \param sigma  The real-valued shift.
74     ///
GenEigsRealShiftSolver(OpType & op,Index nev,Index ncv,const Scalar & sigma)75     GenEigsRealShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma) :
76         Base(op, IdentityBOp(), nev, ncv),
77         m_sigma(sigma)
78     {
79         op.set_shift(m_sigma);
80     }
81 };
82 
83 }  // namespace Spectra
84 
85 #endif  // SPECTRA_GEN_EIGS_REAL_SHIFT_SOLVER_H
86