1 // Copyright (C) 2020-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_SPARSE_GEN_COMPLEX_SHIFT_SOLVE_H 8 #define SPECTRA_SPARSE_GEN_COMPLEX_SHIFT_SOLVE_H 9 10 #include <Eigen/Core> 11 #include <Eigen/SparseCore> 12 #include <Eigen/SparseLU> 13 #include <stdexcept> 14 15 namespace Spectra { 16 17 /// 18 /// \ingroup MatOp 19 /// 20 /// This class defines the complex shift-solve operation on a sparse real matrix \f$A\f$, 21 /// i.e., calculating \f$y=\mathrm{Re}\{(A-\sigma I)^{-1}x\}\f$ for any complex-valued 22 /// \f$\sigma\f$ and real-valued vector \f$x\f$. It is mainly used in the 23 /// GenEigsComplexShiftSolver eigen solver. 24 /// 25 /// \tparam Scalar_ The element type of the matrix, for example, 26 /// `float`, `double`, and `long double`. 27 /// \tparam Flags Either `Eigen::ColMajor` or `Eigen::RowMajor`, indicating 28 /// the storage format of the input matrix. 29 /// \tparam StorageIndex The type of the indices for the sparse matrix. 30 /// 31 template <typename Scalar_, int Flags = Eigen::ColMajor, typename StorageIndex = int> 32 class SparseGenComplexShiftSolve 33 { 34 public: 35 /// 36 /// Element type of the matrix. 37 /// 38 using Scalar = Scalar_; 39 40 private: 41 using Index = Eigen::Index; 42 using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>; 43 using MapConstVec = Eigen::Map<const Vector>; 44 using MapVec = Eigen::Map<Vector>; 45 using SparseMatrix = Eigen::SparseMatrix<Scalar, Flags, StorageIndex>; 46 using ConstGenericSparseMatrix = const Eigen::Ref<const SparseMatrix>; 47 48 using Complex = std::complex<Scalar>; 49 using ComplexVector = Eigen::Matrix<Complex, Eigen::Dynamic, 1>; 50 using SparseComplexMatrix = Eigen::SparseMatrix<Complex, Flags, StorageIndex>; 51 52 using ComplexSolver = Eigen::SparseLU<SparseComplexMatrix>; 53 54 ConstGenericSparseMatrix m_mat; 55 const Index m_n; 56 ComplexSolver m_solver; 57 mutable ComplexVector m_x_cache; 58 59 public: 60 /// 61 /// Constructor to create the matrix operation object. 62 /// 63 /// \param mat An **Eigen** sparse matrix object, whose type can be 64 /// `Eigen::SparseMatrix<Scalar, ...>` or its mapped version 65 /// `Eigen::Map<Eigen::SparseMatrix<Scalar, ...> >`. 66 /// 67 template <typename Derived> SparseGenComplexShiftSolve(const Eigen::SparseMatrixBase<Derived> & mat)68 SparseGenComplexShiftSolve(const Eigen::SparseMatrixBase<Derived>& mat) : 69 m_mat(mat), m_n(mat.rows()) 70 { 71 static_assert( 72 static_cast<int>(Derived::PlainObject::IsRowMajor) == static_cast<int>(SparseMatrix::IsRowMajor), 73 "SparseGenComplexShiftSolve: the \"Flags\" template parameter does not match the input matrix (Eigen::ColMajor/Eigen::RowMajor)"); 74 75 if (mat.rows() != mat.cols()) 76 throw std::invalid_argument("SparseGenComplexShiftSolve: matrix must be square"); 77 } 78 79 /// 80 /// Return the number of rows of the underlying matrix. 81 /// rows()82 Index rows() const { return m_n; } 83 /// 84 /// Return the number of columns of the underlying matrix. 85 /// cols()86 Index cols() const { return m_n; } 87 88 /// 89 /// Set the complex shift \f$\sigma\f$. 90 /// 91 /// \param sigmar Real part of \f$\sigma\f$. 92 /// \param sigmai Imaginary part of \f$\sigma\f$. 93 /// set_shift(const Scalar & sigmar,const Scalar & sigmai)94 void set_shift(const Scalar& sigmar, const Scalar& sigmai) 95 { 96 // Create a sparse idendity matrix (1 + 0i on diagonal) 97 SparseComplexMatrix I(m_n, m_n); 98 I.setIdentity(); 99 // Sparse LU decomposition 100 m_solver.compute(m_mat.template cast<Complex>() - Complex(sigmar, sigmai) * I); 101 // Set cache to zero 102 m_x_cache.resize(m_n); 103 m_x_cache.setZero(); 104 } 105 106 /// 107 /// Perform the complex shift-solve operation 108 /// \f$y=\mathrm{Re}\{(A-\sigma I)^{-1}x\}\f$. 109 /// 110 /// \param x_in Pointer to the \f$x\f$ vector. 111 /// \param y_out Pointer to the \f$y\f$ vector. 112 /// 113 // y_out = Re( inv(A - sigma * I) * x_in ) perform_op(const Scalar * x_in,Scalar * y_out)114 void perform_op(const Scalar* x_in, Scalar* y_out) const 115 { 116 m_x_cache.real() = MapConstVec(x_in, m_n); 117 MapVec y(y_out, m_n); 118 y.noalias() = m_solver.solve(m_x_cache).real(); 119 } 120 }; 121 122 } // namespace Spectra 123 124 #endif // SPECTRA_SPARSE_GEN_COMPLEX_SHIFT_SOLVE_H 125