G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
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
38namespace gismo {
39
40//using Spectra::SELECT_EIGENVALUE; // forward
41
42// Product operation wrapper
43template <class MatrixType> class SpectraMatProd
44{
45public:
46 typedef typename MatrixType::Scalar Scalar;
47 typedef typename MatrixType::Nested NestedMatrix;
48 NestedMatrix m_mat;
49public:
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
63template <class MatrixType> class SpectraSymMatProd
64{
65public:
66 typedef typename MatrixType::Scalar Scalar;
67 typedef typename MatrixType::Nested NestedMatrix;
68 NestedMatrix m_mat;
69public:
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
83template <class MatrixType>
84class SpectraMatShiftSolve
85{
86public:
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
98public:
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)
132template <class T>
133class SpectraMatShiftSolve<gsMatrix<T>>
134{
135public:
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
142public:
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)
165template <class MatrixType>
166class SpectraCholesky
167{
168public:
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
178public:
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
217template <class MatrixType>
218class SpectraRegularInverse
219{
220public:
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
235public:
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
268template <class MatrixType>
269class SpectraSymShiftInvert
270{
271public:
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
283public:
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
332template <class MatrixType>
333class gsSpectraSolver : private SpectraMatProd<MatrixType>,
334 public Spectra::GenEigsSolver<SpectraMatProd<MatrixType> >
335{
336 typedef SpectraMatProd<MatrixType> MatOp;
337 typedef Spectra::GenEigsSolver<MatOp> Base;
338public:
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
345template <class MatrixType>
346class gsSpectraSymSolver : private SpectraSymMatProd<MatrixType>,
347 public Spectra::SymEigsSolver<SpectraSymMatProd<MatrixType> >
348{
349 typedef SpectraSymMatProd<MatrixType> MatOp;
350 typedef Spectra::SymEigsSolver<MatOp> Base;
351public:
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
358template <class MatrixType>
359class 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;
366public:
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
374template <class MatrixType, Spectra::GEigsMode GEigsMode = Spectra::GEigsMode::Cholesky>
375class SpectraOps
376{
377public:
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
386template <class T> class SpectraOps<gsMatrix<T> >
387{
388public:
389 typedef gsMatrix<T> MatrixType;
390 typedef Spectra::DenseCholesky<T> InvOp;
391 typedef SpectraMatProd<MatrixType> ProdOp;
392protected:
393 SpectraOps(const MatrixType & A, const MatrixType & B) : opA(A), opB(B) { }
394 ProdOp opA;
395 InvOp opB;
396};
397
400template <class MatrixType>
401class SpectraOps<MatrixType,Spectra::GEigsMode::RegularInverse>
402{
403public:
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
413template <class MatrixType, Spectra::GEigsMode GEigsMode = Spectra::GEigsMode::ShiftInvert>
414class SpectraShiftOps
415{
416public:
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
428template <class MatrixType>
429class SpectraShiftOps<MatrixType,Spectra::GEigsMode::Buckling>
430{
431public:
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
443template <class MatrixType, Spectra::GEigsMode GEigsMode = Spectra::GEigsMode::Cholesky>
444class 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;
454public:
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
463template <class MatrixType, Spectra::GEigsMode GEigsMode = Spectra::GEigsMode::ShiftInvert>
464class 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;
474public:
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::SpectraOps
Definition gsSpectra.h:376
gismo::SpectraShiftOps
Definition gsSpectra.h:415
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsSpectraGenSymShiftSolver
Definition gsSpectra.h:467
gismo::gsSpectraGenSymSolver
Definition gsSpectra.h:447
gismo::gsSpectraSolver
Eigenvalue solver for general real matrices.
Definition gsSpectra.h:335
gismo::gsSpectraSymShiftSolver
Shifted Eigenvalue solver for real symmetric matrices.
Definition gsSpectra.h:362
gismo::gsSpectraSymSolver
Eigenvalue solver for real symmetric matrices.
Definition gsSpectra.h:348
gsConfig.h
Provides preprocessor directives configuration of G+Smo.
index_t
#define index_t
Definition gsConfig.h:32
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gismo
The G+Smo namespace, containing all definitions for the library.