Spectra  1.0.1
Header-only C++ Library for Large Scale Eigenvalue Problems
DenseSymShiftSolve.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_DENSE_SYM_SHIFT_SOLVE_H
8 #define SPECTRA_DENSE_SYM_SHIFT_SOLVE_H
9 
10 #include <Eigen/Core>
11 #include <stdexcept>
12 
13 #include "../LinAlg/BKLDLT.h"
14 #include "../Util/CompInfo.h"
15 
16 namespace Spectra {
17 
32 template <typename Scalar_, int Uplo = Eigen::Lower, int Flags = Eigen::ColMajor>
34 {
35 public:
39  using Scalar = Scalar_;
40 
41 private:
42  using Index = Eigen::Index;
43  using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Flags>;
44  using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
45  using MapConstVec = Eigen::Map<const Vector>;
46  using MapVec = Eigen::Map<Vector>;
47  using ConstGenericMatrix = const Eigen::Ref<const Matrix>;
48 
49  ConstGenericMatrix m_mat;
50  const Index m_n;
51  BKLDLT<Scalar> m_solver;
52 
53 public:
62  template <typename Derived>
63  DenseSymShiftSolve(const Eigen::MatrixBase<Derived>& mat) :
64  m_mat(mat), m_n(mat.rows())
65  {
66  static_assert(
67  static_cast<int>(Derived::PlainObject::IsRowMajor) == static_cast<int>(Matrix::IsRowMajor),
68  "DenseSymShiftSolve: the \"Flags\" template parameter does not match the input matrix (Eigen::ColMajor/Eigen::RowMajor)");
69 
70  if (m_n != mat.cols())
71  throw std::invalid_argument("DenseSymShiftSolve: matrix must be square");
72  }
73 
77  Index rows() const { return m_n; }
81  Index cols() const { return m_n; }
82 
86  void set_shift(const Scalar& sigma)
87  {
88  m_solver.compute(m_mat, Uplo, sigma);
89  if (m_solver.info() != CompInfo::Successful)
90  throw std::invalid_argument("DenseSymShiftSolve: factorization failed with the given shift");
91  }
92 
99  // y_out = inv(A - sigma * I) * x_in
100  void perform_op(const Scalar* x_in, Scalar* y_out) const
101  {
102  MapConstVec x(x_in, m_n);
103  MapVec y(y_out, m_n);
104  y.noalias() = m_solver.solve(x);
105  }
106 };
107 
108 } // namespace Spectra
109 
110 #endif // SPECTRA_DENSE_SYM_SHIFT_SOLVE_H
void perform_op(const Scalar *x_in, Scalar *y_out) const
DenseSymShiftSolve(const Eigen::MatrixBase< Derived > &mat)
void set_shift(const Scalar &sigma)
@ Successful
Computation was successful.