Spectra  1.0.1
Header-only C++ Library for Large Scale Eigenvalue Problems
SparseSymShiftSolve.h
1 // Copyright (C) 2016-2022 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_SYM_SHIFT_SOLVE_H
8 #define SPECTRA_SPARSE_SYM_SHIFT_SOLVE_H
9 
10 #include <Eigen/Core>
11 #include <Eigen/SparseCore>
12 #include <Eigen/SparseLU>
13 #include <stdexcept>
14 
15 namespace Spectra {
16 
32 template <typename Scalar_, int Uplo = Eigen::Lower, int Flags = Eigen::ColMajor, typename StorageIndex = int>
34 {
35 public:
39  using Scalar = Scalar_;
40 
41 private:
42  using Index = Eigen::Index;
43  using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
44  using MapConstVec = Eigen::Map<const Vector>;
45  using MapVec = Eigen::Map<Vector>;
46  using SparseMatrix = Eigen::SparseMatrix<Scalar, Flags, StorageIndex>;
47  using ConstGenericSparseMatrix = const Eigen::Ref<const SparseMatrix>;
48 
49  ConstGenericSparseMatrix m_mat;
50  const Index m_n;
51  Eigen::SparseLU<SparseMatrix> m_solver;
52 
53 public:
61  template <typename Derived>
62  SparseSymShiftSolve(const Eigen::SparseMatrixBase<Derived>& mat) :
63  m_mat(mat), m_n(mat.rows())
64  {
65  static_assert(
66  static_cast<int>(Derived::PlainObject::IsRowMajor) == static_cast<int>(SparseMatrix::IsRowMajor),
67  "SparseSymShiftSolve: the \"Flags\" template parameter does not match the input matrix (Eigen::ColMajor/Eigen::RowMajor)");
68 
69  if (mat.rows() != mat.cols())
70  throw std::invalid_argument("SparseSymShiftSolve: matrix must be square");
71  }
72 
76  Index rows() const { return m_n; }
80  Index cols() const { return m_n; }
81 
85  void set_shift(const Scalar& sigma)
86  {
87  SparseMatrix mat = m_mat.template selfadjointView<Uplo>();
88  SparseMatrix identity(m_n, m_n);
89  identity.setIdentity();
90  mat = mat - sigma * identity;
91  m_solver.isSymmetric(true);
92  m_solver.compute(mat);
93  if (m_solver.info() != Eigen::Success)
94  throw std::invalid_argument("SparseSymShiftSolve: factorization failed with the given shift");
95  }
96 
103  // y_out = inv(A - sigma * I) * x_in
104  void perform_op(const Scalar* x_in, Scalar* y_out) const
105  {
106  MapConstVec x(x_in, m_n);
107  MapVec y(y_out, m_n);
108  y.noalias() = m_solver.solve(x);
109  }
110 };
111 
112 } // namespace Spectra
113 
114 #endif // SPECTRA_SPARSE_SYM_SHIFT_SOLVE_H
SparseSymShiftSolve(const Eigen::SparseMatrixBase< Derived > &mat)
void perform_op(const Scalar *x_in, Scalar *y_out) const
void set_shift(const Scalar &sigma)