DenseGenComplexShiftSolve.h
Go to the documentation of this file.
1 // Copyright (C) 2016-2019 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 DENSE_GEN_COMPLEX_SHIFT_SOLVE_H
8 #define DENSE_GEN_COMPLEX_SHIFT_SOLVE_H
9 
10 #include <Eigen/Core>
11 #include <Eigen/LU>
12 #include <stdexcept>
13 
14 namespace Spectra {
15 
24 template <typename Scalar>
26 {
27 private:
34 
35  typedef std::complex<Scalar> Complex;
38 
40 
41  ConstGenericMatrix m_mat;
42  const Index m_n;
43  ComplexSolver m_solver;
44  ComplexVector m_x_cache;
45 
46 public:
55  DenseGenComplexShiftSolve(ConstGenericMatrix& mat) :
56  m_mat(mat), m_n(mat.rows())
57  {
58  if (mat.rows() != mat.cols())
59  throw std::invalid_argument("DenseGenComplexShiftSolve: matrix must be square");
60  }
61 
65  Index rows() const { return m_n; }
69  Index cols() const { return m_n; }
70 
77  void set_shift(Scalar sigmar, Scalar sigmai)
78  {
79  m_solver.compute(m_mat.template cast<Complex>() - Complex(sigmar, sigmai) * ComplexMatrix::Identity(m_n, m_n));
80  m_x_cache.resize(m_n);
81  m_x_cache.setZero();
82  }
83 
91  // y_out = Re( inv(A - sigma * I) * x_in )
92  void perform_op(const Scalar* x_in, Scalar* y_out)
93  {
94  m_x_cache.real() = MapConstVec(x_in, m_n);
95  MapVec y(y_out, m_n);
96  y.noalias() = m_solver.solve(m_x_cache).real();
97  }
98 };
99 
100 } // namespace Spectra
101 
102 #endif // DENSE_GEN_COMPLEX_SHIFT_SOLVE_H
DenseGenComplexShiftSolve(ConstGenericMatrix &mat)
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SolverBase.h:106
SCALAR Scalar
Definition: bench_gemm.cpp:46
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
PartialPivLU & compute(const EigenBase< InputType > &matrix)
Definition: PartialPivLU.h:131
Scalar * y
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
const Eigen::Ref< const Matrix > ConstGenericMatrix
void set_shift(Scalar sigmar, Scalar sigmai)
void perform_op(const Scalar *x_in, Scalar *y_out)
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > Matrix
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > Vector
Eigen::Matrix< Complex, Eigen::Dynamic, Eigen::Dynamic > ComplexMatrix
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Eigen::PartialPivLU< ComplexMatrix > ComplexSolver
Eigen::Matrix< Complex, Eigen::Dynamic, 1 > ComplexVector


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:09