G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsSpectra.h
Go to the documentation of this file.
1 
24 #pragma once
25 
26 #include <gsCore/gsConfig.h>
27 #include <gsCore/gsLinearAlgebra.h>
28 
29 #define Eigen gsEigen
30 
31 #include <Spectra/include/Spectra/SymEigsSolver.h>
32 #include <Spectra/include/Spectra/SymEigsShiftSolver.h>
33 #include <Spectra/include/Spectra/SymGEigsSolver.h>
34 #include <Spectra/include/Spectra/SymGEigsShiftSolver.h>
35 #include <Spectra/include/Spectra/GenEigsSolver.h>
36 #include <Spectra/include/Spectra/MatOp/SparseSymShiftSolve.h>
37 
38 namespace gismo {
39 
40 //using Spectra::SELECT_EIGENVALUE; // forward
41 
42 // Product operation wrapper
43 template <class MatrixType> class SpectraMatProd
44 {
45 public:
46  typedef typename MatrixType::Scalar Scalar;
47  typedef typename MatrixType::Nested NestedMatrix;
48  NestedMatrix m_mat;
49 public:
50  SpectraMatProd(const MatrixType&& ) = delete;
51  SpectraMatProd(const MatrixType& mat) : m_mat(mat) {}
52  index_t rows() const { return m_mat.rows(); }
53  index_t cols() const { return m_mat.cols(); }
54  void perform_op(const Scalar* x_in, Scalar* y_out) const
55  {
56  GISMO_ASSERT(m_mat.rows()!=0 && m_mat.cols()!=0,"The matrix has zero rows or columns. Is the matrix a temporary (e.g. A-B)?");
57  gsAsVector<Scalar>(y_out, m_mat.rows()).noalias() =
58  m_mat * gsAsConstVector<Scalar>(x_in, m_mat.cols());
59  }
60 };
61 
62 // Product operation wrapper
63 template <class MatrixType> class SpectraSymMatProd
64 {
65 public:
66  typedef typename MatrixType::Scalar Scalar;
67  typedef typename MatrixType::Nested NestedMatrix;
68  NestedMatrix m_mat;
69 public:
70  SpectraSymMatProd(const MatrixType&& ) = delete;
71  SpectraSymMatProd(const MatrixType& mat) : m_mat(mat) {}
72  index_t rows() const { return m_mat.rows(); }
73  index_t cols() const { return m_mat.cols(); }
74  void perform_op(const Scalar* x_in, Scalar* y_out) const
75  {
76  GISMO_ASSERT(m_mat.rows()!=0 && m_mat.cols()!=0,"The matrix has zero rows or columns. Is the matrix a temporary (e.g. A-B)?");
77  gsAsVector<Scalar>(y_out, m_mat.rows()).noalias() =
78  m_mat.template selfadjointView<gsEigen::Lower>() * gsAsConstVector<Scalar>(x_in, m_mat.cols());
79  }
80 };
81 
82 // Shift operation wrapper
83 template <class MatrixType>
84 class SpectraMatShiftSolve
85 {
86 public:
87  typedef typename MatrixType::Scalar Scalar;
88  typedef typename MatrixType::Nested NestedMatrix;
89  NestedMatrix m_mat;
90  const index_t m_n;
91 #ifdef GISMO_WITH_PARDISO
92  typename gsSparseSolver<Scalar>::PardisoLU m_solver;
93 #elif GISMO_WITH_SUPERLU
94  typename gsSparseSolver<Scalar>::SuperLU m_solver;
95 #else
96  typename gsSparseSolver<Scalar>::LU m_solver;
97 #endif
98 public:
99  SpectraMatShiftSolve(const MatrixType&& ) = delete;
100  SpectraMatShiftSolve(const MatrixType& mat)
101  :
102  m_mat(mat), m_n(mat.rows())
103  {
104  GISMO_ASSERT(m_mat.rows() == m_mat.cols(),"Matrix must be square!");
105  }
106  index_t rows() const { return m_n; }
107  index_t cols() const { return m_n; }
108  void perform_op(const Scalar* x_in, Scalar* y_out) const
109  {
110  gsAsVector<Scalar>(y_out, m_n).noalias() = m_solver.solve( gsAsConstVector<Scalar>(x_in, m_n) );
111  }
112 
113  void set_shift(const Scalar& sigma)
114  {
115  MatrixType mat(m_mat.template selfadjointView<gsEigen::Lower>());
116  MatrixType identity(m_n,m_n);
117  identity.setIdentity();
118 
119  mat = mat - sigma * identity;
120 #ifndef GISMO_WITH_PARDISO
121 #ifndef GISMO_WITH_SUPERLU
122  m_solver.isSymmetric(true);
123 #endif
124 #endif
125  m_solver.compute(mat);
126 
127  GISMO_ASSERT(m_solver.info() == gsEigen::Success,"SpectraMatShiftSolve: factorization failed with the given shift");
128  }
129 };
130 
131 // Shift operation wrapper (dense matrix)
132 template <class T>
133 class SpectraMatShiftSolve<gsMatrix<T>>
134 {
135 public:
136  typedef typename gsMatrix<T>::Nested NestedMatrix;
137  typedef typename gsMatrix<T>::Scalar Scalar;
138  NestedMatrix m_mat;
139  const index_t m_n;
140  typename Spectra::BKLDLT<Scalar> m_solver;
141 
142 public:
143  SpectraMatShiftSolve(const gsMatrix<Scalar>&&) = delete;
144  SpectraMatShiftSolve(const gsMatrix<Scalar>& mat)
145  :
146  m_mat(mat), m_n(mat.rows())
147  {
148  GISMO_ASSERT(m_mat.rows() == m_mat.cols(),"Matrix must be square!");
149  }
150  index_t rows() const { return m_n; }
151  index_t cols() const { return m_n; }
152  void perform_op(const Scalar* x_in, Scalar* y_out) const
153  {
154  gsAsVector<Scalar>(y_out, m_n).noalias() = m_solver.solve( gsAsConstVector<Scalar>(x_in, m_n) );
155  }
156 
157  void set_shift(const Scalar& sigma)
158  {
159  m_solver.compute(m_mat, gsEigen::Lower, sigma);
160  GISMO_ASSERT(m_solver.info() == Spectra::CompInfo::Successful,"SpectraMatShiftSolve: factorization failed with the given shift");
161  }
162 };
163 
164 // Cholesky wrapper (see Spectra::SparseCholesky)
165 template <class MatrixType>
166 class SpectraCholesky
167 {
168 public:
169  typedef typename MatrixType::Scalar Scalar;
170  typedef typename MatrixType::Nested NestedMatrix;
171  NestedMatrix m_mat;
172  const index_t m_n;
173  // NOTE: Does not work for gsSparseSolver<Scalar>::SimplicialLLT!
174  gsEigen::SimplicialLLT<gsSparseMatrix<Scalar>, gsEigen::Lower> m_solver;
175 
176  Spectra::CompInfo m_info; // status of the decomposition
177 
178 public:
179  SpectraCholesky(const MatrixType&& ) = delete;
180  SpectraCholesky(const MatrixType& mat)
181  :
182  m_mat(mat), m_n(mat.rows())
183  {
184  GISMO_ASSERT(m_mat.rows() == m_mat.cols(),"Matrix A must be square!");
185 
186  m_solver.compute(mat);
187  m_info = (m_solver.info() == gsEigen::Success) ?
188  Spectra::CompInfo::Successful :
189  Spectra::CompInfo::NumericalIssue;
190  }
191  index_t rows() const { return m_n; }
192  index_t cols() const { return m_n; }
193 
195  Spectra::CompInfo info() const { return m_info; }
196 
198  void lower_triangular_solve(const Scalar* x_in, Scalar* y_out) const
199  {
200  gsAsConstVector<Scalar> x(x_in, m_n);
201  gsAsVector<Scalar> y(y_out, m_n);
202  y.noalias() = m_solver.permutationP() * x;
203  m_solver.matrixL().solveInPlace(y);
204  }
205 
207  void upper_triangular_solve(const Scalar* x_in, Scalar* y_out) const
208  {
209  gsAsConstVector<Scalar> x(x_in, m_n);
210  gsAsVector<Scalar> y(y_out, m_n);
211  y.noalias() = m_solver.matrixU().solve(x);
212  y = (m_solver.permutationPinv() * y).eval();
213  }
214 };
215 
216 // Regular inverse wrapper
217 template <class MatrixType>
218 class SpectraRegularInverse
219 {
220 public:
221  typedef typename MatrixType::Scalar Scalar;
222  typedef typename MatrixType::Nested NestedMatrix;
223  NestedMatrix m_mat;
224  const index_t m_n;
225 #ifdef GISMO_WITH_PARDISO
226  typename gsSparseSolver<Scalar>::PardisoLU m_solver;
227 #elif GISMO_WITH_SUPERLU
228  typename gsSparseSolver<Scalar>::SuperLU m_solver;
229 #else
230  typename gsSparseSolver<Scalar>::CGDiagonal m_solver;
231 #endif
232 
233  Spectra::CompInfo m_info; // status of the decomposition
234 
235 public:
236  SpectraRegularInverse(const MatrixType&& ) = delete;
237  SpectraRegularInverse(const MatrixType& mat)
238  :
239  m_mat(mat), m_n(mat.rows())
240  {
241  GISMO_ASSERT(m_mat.rows() == m_mat.cols(),"Matrix A must be square!");
242 
243  m_solver.compute(mat);
244  m_info = (m_solver.info() == gsEigen::Success) ?
245  Spectra::CompInfo::Successful :
246  Spectra::CompInfo::NumericalIssue;
247  }
248  index_t rows() const { return m_n; }
249  index_t cols() const { return m_n; }
250 
252  Spectra::CompInfo info() const { return m_info; }
253 
255  void solve(const Scalar* x_in, Scalar* y_out) const
256  {
257  gsAsVector<Scalar>(y_out, m_n).noalias() = m_solver.solve( gsAsConstVector<Scalar>(x_in, m_n) );
258  }
259 
261  void perform_op(const Scalar* x_in, Scalar* y_out) const
262  {
263  gsAsVector<Scalar>(y_out, m_n).noalias() = m_mat.template selfadjointView<gsEigen::Lower>() * gsAsConstVector<Scalar>(x_in, m_n);
264  }
265 };
266 
267 // Shift operation wrapper
268 template <class MatrixType>
269 class SpectraSymShiftInvert
270 {
271 public:
272  typedef typename MatrixType::Scalar Scalar;
273  typedef typename MatrixType::Nested NestedMatrix;
274  NestedMatrix m_A, m_B;
275  const index_t m_n;
276 #ifdef GISMO_WITH_PARDISO
277  typename gsSparseSolver<Scalar>::PardisoLU m_solver;
278 #elif GISMO_WITH_SUPERLU
279  typename gsSparseSolver<Scalar>::SuperLU m_solver;
280 #else
281  typename gsSparseSolver<Scalar>::LU m_solver;
282 #endif
283 public:
284  SpectraSymShiftInvert(const MatrixType&& ,const MatrixType&& ) = delete;
285  SpectraSymShiftInvert(const MatrixType& A ,const MatrixType& B)
286  :
287  m_A(A), m_B(B), m_n(A.rows())
288  {
289  GISMO_ASSERT(m_A.rows() == m_A.cols(),"Matrix A must be square!");
290  GISMO_ASSERT(m_B.rows() == m_B.cols(),"Matrix B must be square!");
291  GISMO_ASSERT(m_B.rows() == m_A.rows(),"Matrix A and B must be the same size!");
292  }
293  index_t rows() const { return m_n; }
294  index_t cols() const { return m_n; }
295  void perform_op(const Scalar* x_in, Scalar* y_out) const
296  {
297  gsAsVector<Scalar>(y_out, m_n).noalias() = m_solver.solve( gsAsConstVector<Scalar>(x_in, m_n) );
298  }
299 
300  void set_shift(const Scalar& sigma)
301  {
302  MatrixType matA = m_A.template selfadjointView<gsEigen::Lower>();
303  MatrixType matB = m_B.template selfadjointView<gsEigen::Lower>();
304  MatrixType mat = matA - sigma * matB;
305 
306  #ifndef GISMO_WITH_PARDISO
307  #ifndef GISMO_WITH_SUPERLU
308  m_solver.isSymmetric(true);
309  #endif
310  #endif
311  m_solver.compute(mat);
312  GISMO_ASSERT(m_solver.info()==gsEigen::Success,"Factorization failed");
313  }
314 };
315 
316 
332 template <class MatrixType>
333 class gsSpectraSolver : private SpectraMatProd<MatrixType>,
334  public Spectra::GenEigsSolver<SpectraMatProd<MatrixType> >
335 {
336  typedef SpectraMatProd<MatrixType> MatOp;
337  typedef Spectra::GenEigsSolver<MatOp> Base;
338 public:
339  gsSpectraSolver(const MatrixType && , index_t nev_, index_t ncv_) = delete;
340  gsSpectraSolver(const MatrixType & mat, index_t nev_, index_t ncv_) :
341  MatOp(mat), Base(*this, nev_, ncv_) { Base::init(); }
342 };
343 
345 template <class MatrixType>
346 class gsSpectraSymSolver : private SpectraSymMatProd<MatrixType>,
347  public Spectra::SymEigsSolver<SpectraSymMatProd<MatrixType> >
348 {
349  typedef SpectraSymMatProd<MatrixType> MatOp;
350  typedef Spectra::SymEigsSolver<MatOp> Base;
351 public:
352  gsSpectraSymSolver(const MatrixType && , index_t nev_, index_t ncv_) = delete;
353  gsSpectraSymSolver(const MatrixType & mat, index_t nev_, index_t ncv_) :
354  MatOp(mat), Base(*this, nev_, ncv_) { Base::init(); }
355 };
356 
358 template <class MatrixType>
359 class gsSpectraSymShiftSolver :
360  private SpectraMatShiftSolve<MatrixType>,
361  public Spectra::SymEigsShiftSolver<SpectraMatShiftSolve<MatrixType> >
362 {
363  typedef typename MatrixType::Scalar Scalar;
364  typedef SpectraMatShiftSolve<MatrixType> Op;
365  typedef Spectra::SymEigsShiftSolver<Op> Base;
366 public:
367  gsSpectraSymShiftSolver(const MatrixType && , index_t nev_, index_t ncv_, const Scalar& sigma) = delete;
368  gsSpectraSymShiftSolver(const MatrixType & mat, index_t nev_, index_t ncv_, const Scalar& sigma) :
369  Op(mat), Base(*this, nev_, ncv_,sigma) { Base::init(); }
370 };
371 
374 template <class MatrixType, Spectra::GEigsMode GEigsMode = Spectra::GEigsMode::Cholesky>
375 class SpectraOps
376 {
377 public:
378  typedef SpectraSymMatProd<MatrixType> ProdOp;
379  typedef SpectraCholesky<MatrixType> InvOp;
380  SpectraOps(const MatrixType & A, const MatrixType & B) : opA(A), opB(B) { }
381  ProdOp opA;
382  InvOp opB;
383 };
384 
385 //template<> //compilation fails with this
386 template <class T> class SpectraOps<gsMatrix<T> >
387 {
388 public:
389  typedef gsMatrix<T> MatrixType;
390  typedef Spectra::DenseCholesky<T> InvOp;
391  typedef SpectraMatProd<MatrixType> ProdOp;
392 protected:
393  SpectraOps(const MatrixType & A, const MatrixType & B) : opA(A), opB(B) { }
394  ProdOp opA;
395  InvOp opB;
396 };
397 
400 template <class MatrixType>
401 class SpectraOps<MatrixType,Spectra::GEigsMode::RegularInverse>
402 {
403 public:
404  typedef SpectraRegularInverse<MatrixType> InvOp;
405  typedef SpectraSymMatProd<MatrixType> ProdOp;
406  SpectraOps(const MatrixType & A, const MatrixType & B) : opA(A), opB(B) { }
407  ProdOp opA;
408  InvOp opB;
409 };
410 
413 template <class MatrixType, Spectra::GEigsMode GEigsMode = Spectra::GEigsMode::ShiftInvert>
414 class SpectraShiftOps
415 {
416 public:
417  typedef typename MatrixType::Scalar Scalar;
418  typedef SpectraSymShiftInvert<MatrixType> InvOp;
419  typedef SpectraSymMatProd<MatrixType> ProdOp;
420 
421  SpectraShiftOps(const MatrixType & A, const MatrixType & B) : opA(A,B), opB(B) { }
422  InvOp opA;
423  ProdOp opB;
424 };
425 
428 template <class MatrixType>
429 class SpectraShiftOps<MatrixType,Spectra::GEigsMode::Buckling>
430 {
431 public:
432  typedef typename MatrixType::Scalar Scalar;
433  typedef SpectraSymShiftInvert<MatrixType> InvOp;
434  typedef SpectraSymMatProd<MatrixType> ProdOp;
435 
436  SpectraShiftOps(const MatrixType & A, const MatrixType & B) : opA(A,B), opB(A) { }
437  InvOp opA;
438  ProdOp opB;
439 };
440 
443 template <class MatrixType, Spectra::GEigsMode GEigsMode = Spectra::GEigsMode::Cholesky>
444 class gsSpectraGenSymSolver :
445  private SpectraOps<MatrixType,GEigsMode>,
446  public Spectra::SymGEigsSolver<typename SpectraOps<MatrixType,GEigsMode>::ProdOp, typename SpectraOps<MatrixType,GEigsMode>::InvOp, GEigsMode>
447 {
448  typedef typename MatrixType::Scalar Scalar;
449  typedef SpectraOps<MatrixType,GEigsMode> OpType;
450  typedef typename OpType::ProdOp Op;
451  typedef typename OpType::InvOp BOp;
452 
453  typedef Spectra::SymGEigsSolver<Op, BOp,GEigsMode> Base;
454 public:
455  gsSpectraGenSymSolver(const MatrixType && , const MatrixType && , index_t nev_, index_t ncv_) = delete;
456  gsSpectraGenSymSolver(const MatrixType & Amat, const MatrixType & Bmat, index_t nev_, index_t ncv_)
457  : OpType(Amat,Bmat), Base(this->opA, this->opB, nev_, math::min(ncv_,Amat.rows()))
458  { Base::init(); }
459 };
460 
463 template <class MatrixType, Spectra::GEigsMode GEigsMode = Spectra::GEigsMode::ShiftInvert>
464 class gsSpectraGenSymShiftSolver :
465  private SpectraShiftOps<MatrixType,GEigsMode>,
466  public Spectra::SymGEigsShiftSolver<typename SpectraShiftOps<MatrixType,GEigsMode>::InvOp, typename SpectraShiftOps<MatrixType,GEigsMode>::ProdOp, GEigsMode>
467 {
468  typedef typename MatrixType::Scalar Scalar;
469  typedef SpectraShiftOps<MatrixType,GEigsMode> OpType;
470  typedef typename OpType::InvOp Op;
471  typedef typename OpType::ProdOp BOp;
472 
473  typedef Spectra::SymGEigsShiftSolver<Op, BOp,GEigsMode> Base;
474 public:
475  gsSpectraGenSymShiftSolver(const MatrixType && , const MatrixType && , index_t nev_, index_t ncv_, const Scalar& sigma) = delete;
476  gsSpectraGenSymShiftSolver(const MatrixType & Amat, const MatrixType & Bmat, index_t nev_, index_t ncv_, const Scalar& sigma)
477  : OpType(Amat,Bmat), Base(this->opA, this->opB, nev_, math::min(ncv_,Amat.rows()),sigma)
478  { Base::init(); }
479 };
480 
481 } //namespace gismo
482 
483 
484 #undef Eigen
485 
gismo::gsSpectraSymShiftSolver
Shifted Eigenvalue solver for real symmetric matrices.
Definition: gsSpectra.h:359
gismo::gsSpectraGenSymSolver
Definition: gsSpectra.h:444
index_t
#define index_t
Definition: gsConfig.h:32
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition: gsMatrix.h:38
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gsConfig.h
Provides preprocessor directives configuration of G+Smo.
gismo::SpectraOps
Definition: gsSpectra.h:375
gismo::gsSpectraSymSolver
Eigenvalue solver for real symmetric matrices.
Definition: gsSpectra.h:346
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gismo::gsSpectraGenSymShiftSolver
Definition: gsSpectra.h:464
gismo::SpectraShiftOps
Definition: gsSpectra.h:414
gismo::gsSpectraSolver
Eigenvalue solver for general real matrices.
Definition: gsSpectra.h:333