Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #ifndef PALACE_LINALG_SLEPC_HPP
5 : #define PALACE_LINALG_SLEPC_HPP
6 :
7 : #if defined(PALACE_WITH_SLEPC)
8 :
9 : #include "linalg/petsc.hpp"
10 :
11 : #if !defined(PETSC_USE_COMPLEX)
12 : #error "SLEPc interface requires PETSc compiled with complex scalars!"
13 : #endif
14 :
15 : #include <complex>
16 : #include <memory>
17 : #include <optional>
18 : #include <string>
19 : #include <mpi.h>
20 : #include "linalg/eps.hpp"
21 : #include "linalg/ksp.hpp"
22 : #include "linalg/operator.hpp"
23 : #include "linalg/vector.hpp"
24 :
25 : // Forward declarations of SLEPc objects.
26 : typedef struct _p_EPS *EPS;
27 : typedef struct _p_PEP *PEP;
28 : typedef struct _p_NEP *NEP;
29 : typedef struct _p_BV *BV;
30 : typedef struct _p_ST *ST;
31 : typedef struct _p_RG *RG;
32 :
33 : namespace palace
34 : {
35 :
36 : namespace slepc
37 : {
38 :
39 : // Wrappers for SlepcInitialize/SlepcInitializeNoArguments/SlepcFinalize.
40 : void Initialize(int &argc, char **&argv, const char rc_file[], const char help[]);
41 : void Initialize();
42 : void Finalize();
43 :
44 : // Compute and return the maximum singular value of the given operator, σₙ² = λₙ(Aᴴ A) .
45 : PetscReal GetMaxSingularValue(MPI_Comm comm, const ComplexOperator &A, bool herm = false,
46 : PetscReal tol = PETSC_DEFAULT,
47 : PetscInt max_it = PETSC_DEFAULT);
48 :
49 : //
50 : // A wrapper for the SLEPc library for generalized linear eigenvalue problems or quadratic
51 : // polynomial eigenvalue problems. Shift-and-invert spectral transformations can be used to
52 : // compute interior eigenvalues.
53 : //
54 : class SlepcEigenvalueSolver : public EigenvalueSolver
55 : {
56 : public:
57 : enum class ProblemType
58 : {
59 : HERMITIAN,
60 : NON_HERMITIAN,
61 : GEN_HERMITIAN,
62 : GEN_INDEFINITE,
63 : GEN_NON_HERMITIAN,
64 : HYPERBOLIC,
65 : GYROSCOPIC,
66 : GENERAL,
67 : };
68 :
69 : enum class Type
70 : {
71 : KRYLOVSCHUR,
72 : POWER,
73 : SUBSPACE,
74 : TOAR,
75 : STOAR,
76 : QARNOLDI,
77 : JD,
78 : SLP,
79 : NLEIGS
80 : };
81 :
82 : // Workspace vector for operator applications.
83 : mutable ComplexVector x1, y1;
84 :
85 : // References to matrices defining the (possibly quadratic) eigenvalue problem
86 : // (not owned).
87 : const ComplexOperator *opK, *opC, *opM;
88 :
89 : protected:
90 : // Control print level for debugging.
91 : int print;
92 :
93 : // Variables for scaling, from Higham et al., IJNME 2008.
94 : PetscReal gamma, delta;
95 :
96 : // Parameters defining the spectral transformation.
97 : PetscScalar sigma;
98 : bool sinvert, region;
99 :
100 : // Storage for computed residual norms and eigenvector normalizations.
101 : std::unique_ptr<PetscReal[]> res, xscale;
102 :
103 : // Reference to linear solver used for operator action for M⁻¹ (with no spectral
104 : // transformation) or (K - σ M)⁻¹ (generalized EVP with shift-and- invert) or P(σ)⁻¹
105 : // (polynomial with shift-and-invert) (not owned).
106 : ComplexKspSolver *opInv;
107 :
108 : // Reference to solver for projecting an intermediate vector onto a divergence-free space
109 : // (not owned).
110 : const DivFreeSolver<ComplexVector> *opProj;
111 :
112 : // Reference to matrix used for weighted inner products (not owned). May be nullptr, in
113 : // which case identity is used.
114 : const Operator *opB;
115 :
116 : // Workspace objects for eigenvalue calculations.
117 : Mat B0;
118 : Vec v0;
119 :
120 : // Boolean to handle SetFromOptions calls.
121 : mutable bool cl_custom;
122 :
123 : // Customize object with command line options set.
124 : virtual void Customize();
125 :
126 : // Helper routine for computing the eigenvector normalization.
127 : PetscReal GetEigenvectorNorm(const ComplexVector &x, ComplexVector &Bx) const;
128 :
129 : // Helper routine for computing the eigenpair residual.
130 : virtual PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x,
131 : ComplexVector &r) const = 0;
132 :
133 : // Helper routine for computing the backward error.
134 : virtual PetscReal GetBackwardScaling(PetscScalar l) const = 0;
135 :
136 : public:
137 : SlepcEigenvalueSolver(int print);
138 : ~SlepcEigenvalueSolver() override;
139 :
140 : // Set operators for the generalized eigenvalue problem, the quadratic polynomial
141 : // eigenvalue problem, or the nonlinear eigenvalue problem.
142 : void SetOperators(const ComplexOperator &K, const ComplexOperator &M,
143 : ScaleType type) override;
144 : void SetOperators(const ComplexOperator &K, const ComplexOperator &C,
145 : const ComplexOperator &M, ScaleType type) override;
146 :
147 : // For the linear generalized case, the linear solver should be configured to compute the
148 : // action of M⁻¹ (with no spectral transformation) or (K - σ M)⁻¹. For the quadratic
149 : // case, the linear solver should be configured to compute the action of M⁻¹ (with no
150 : // spectral transformation) or P(σ)⁻¹.
151 : void SetLinearSolver(ComplexKspSolver &ksp) override;
152 :
153 : // Set the projection operator for enforcing the divergence-free constraint.
154 : void SetDivFreeProjector(const DivFreeSolver<ComplexVector> &divfree) override;
155 :
156 : // Set optional B matrix used for weighted inner products. This must be set explicitly
157 : // even for generalized problems, otherwise the identity will be used.
158 : void SetBMat(const Operator &B) override;
159 :
160 : // Get scaling factors used by the solver.
161 0 : PetscReal GetScalingGamma() const override { return gamma; }
162 0 : PetscReal GetScalingDelta() const override { return delta; }
163 :
164 : // Set shift-and-invert spectral transformation.
165 : void SetShiftInvert(std::complex<double> s, bool precond = false) override;
166 :
167 : // Set problem type.
168 : virtual void SetProblemType(ProblemType type) = 0;
169 :
170 : // Set eigenvalue solver.
171 : virtual void SetType(Type type) = 0;
172 :
173 : // Configure the basis vectors object associated with the eigenvalue solver.
174 : void SetOrthogonalization(bool mgs, bool cgs2);
175 :
176 : // Get the corresponding eigenpair error.
177 : PetscReal GetError(int i, ErrorType type) const override;
178 :
179 : // Re-normalize the given number of eigenvectors, for example if the matrix B for weighted
180 : // inner products has changed. This does not perform re-orthogonalization with respect to
181 : // the new matrix, only normalization.
182 : void RescaleEigenvectors(int num_eig) override;
183 :
184 : // Get the basis vectors object.
185 : virtual BV GetBV() const = 0;
186 :
187 : // Get the spectral transformation object.
188 : virtual ST GetST() const = 0;
189 :
190 : // Get the filtering region object.
191 : virtual RG GetRG() const = 0;
192 :
193 : // Get the associated MPI communicator.
194 : virtual MPI_Comm GetComm() const = 0;
195 :
196 : // Conversion function to PetscObject.
197 : virtual operator PetscObject() const = 0;
198 : };
199 :
200 : // Base class for SLEPc's EPS problem type.
201 : class SlepcEPSSolverBase : public SlepcEigenvalueSolver
202 : {
203 : protected:
204 : // SLEPc eigensolver object. Polynomial problems are handled using linearization.
205 : EPS eps;
206 :
207 : // Shell matrices for the generalized eigenvalue problem.
208 : Mat A0, A1;
209 :
210 : void Customize() override;
211 :
212 : public:
213 : // Calls SLEPc's EPSCreate. Expects SLEPc to be initialized/finalized externally.
214 : SlepcEPSSolverBase(MPI_Comm comm, int print, const std::string &prefix = std::string());
215 :
216 : // Call's SLEPc's EPSDestroy.
217 : ~SlepcEPSSolverBase() override;
218 :
219 : // Conversion function to SLEPc's EPS type.
220 : operator EPS() const { return eps; }
221 :
222 : void SetNumModes(int num_eig, int num_vec = 0) override;
223 :
224 : void SetTol(PetscReal tol) override;
225 :
226 : void SetMaxIter(int max_it) override;
227 :
228 : void SetWhichEigenpairs(WhichType type) override;
229 :
230 : void SetProblemType(ProblemType type) override;
231 :
232 : void SetType(Type type) override;
233 :
234 : void SetInitialSpace(const ComplexVector &v) override;
235 :
236 : int Solve() override;
237 :
238 : std::complex<double> GetEigenvalue(int i) const override;
239 :
240 : void GetEigenvector(int i, ComplexVector &x) const override;
241 :
242 : BV GetBV() const override;
243 :
244 : ST GetST() const override;
245 :
246 : RG GetRG() const override;
247 :
248 0 : MPI_Comm GetComm() const override
249 : {
250 0 : return eps ? PetscObjectComm(reinterpret_cast<PetscObject>(eps)) : MPI_COMM_NULL;
251 : }
252 :
253 0 : operator PetscObject() const override { return reinterpret_cast<PetscObject>(eps); };
254 : };
255 :
256 : // Generalized eigenvalue problem solver: K x = λ M x .
257 : class SlepcEPSSolver : public SlepcEPSSolverBase
258 : {
259 : public:
260 : using SlepcEigenvalueSolver::delta;
261 : using SlepcEigenvalueSolver::gamma;
262 : using SlepcEigenvalueSolver::opB;
263 : using SlepcEigenvalueSolver::opInv;
264 : using SlepcEigenvalueSolver::opProj;
265 : using SlepcEigenvalueSolver::sigma;
266 : using SlepcEigenvalueSolver::sinvert;
267 :
268 : private:
269 : // Operator norms for scaling.
270 : mutable PetscReal normK, normM;
271 :
272 : protected:
273 : PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x,
274 : ComplexVector &r) const override;
275 :
276 : PetscReal GetBackwardScaling(PetscScalar l) const override;
277 :
278 : public:
279 : SlepcEPSSolver(MPI_Comm comm, int print, const std::string &prefix = std::string());
280 :
281 : using SlepcEigenvalueSolver::SetOperators;
282 : void SetOperators(const ComplexOperator &K, const ComplexOperator &M,
283 : ScaleType type) override;
284 :
285 : void SetBMat(const Operator &B) override;
286 : };
287 :
288 : // Quadratic eigenvalue problem solver: P(λ) x = (K + λ C + λ² M) x = 0 , solved via
289 : // linearization: L₀ y = λ L₁ y .
290 : class SlepcPEPLinearSolver : public SlepcEPSSolverBase
291 : {
292 : public:
293 : using SlepcEigenvalueSolver::delta;
294 : using SlepcEigenvalueSolver::gamma;
295 : using SlepcEigenvalueSolver::opB;
296 : using SlepcEigenvalueSolver::opInv;
297 : using SlepcEigenvalueSolver::opProj;
298 : using SlepcEigenvalueSolver::sigma;
299 : using SlepcEigenvalueSolver::sinvert;
300 :
301 : // Workspace vectors for operator applications.
302 : mutable ComplexVector x2, y2;
303 :
304 : private:
305 : // Operator norms for scaling.
306 : mutable PetscReal normK, normC, normM;
307 :
308 : protected:
309 : PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x,
310 : ComplexVector &r) const override;
311 :
312 : PetscReal GetBackwardScaling(PetscScalar l) const override;
313 :
314 : public:
315 : SlepcPEPLinearSolver(MPI_Comm comm, int print, const std::string &prefix = std::string());
316 :
317 : using SlepcEigenvalueSolver::SetOperators;
318 : void SetOperators(const ComplexOperator &K, const ComplexOperator &C,
319 : const ComplexOperator &M, ScaleType type) override;
320 :
321 : void SetBMat(const Operator &B) override;
322 :
323 : void SetInitialSpace(const ComplexVector &v) override;
324 :
325 : void GetEigenvector(int i, ComplexVector &x) const override;
326 : };
327 :
328 : // Base class for SLEPc's PEP problem type.
329 : class SlepcPEPSolverBase : public SlepcEigenvalueSolver
330 : {
331 : protected:
332 : // SLEPc eigensolver object.
333 : PEP pep;
334 :
335 : // Shell matrices for the quadratic polynomial eigenvalue problem.
336 : Mat A0, A1, A2;
337 :
338 : void Customize() override;
339 :
340 : public:
341 : // Calls SLEPc's PEPCreate. Expects SLEPc to be initialized/finalized externally.
342 : SlepcPEPSolverBase(MPI_Comm comm, int print, const std::string &prefix = std::string());
343 :
344 : // Call's SLEPc's PEPDestroy.
345 : ~SlepcPEPSolverBase() override;
346 :
347 : // Conversion function to SLEPc's PEP type.
348 : operator PEP() const { return pep; }
349 :
350 : void SetNumModes(int num_eig, int num_vec = 0) override;
351 :
352 : void SetTol(PetscReal tol) override;
353 :
354 : void SetMaxIter(int max_it) override;
355 :
356 : void SetWhichEigenpairs(WhichType type) override;
357 :
358 : void SetProblemType(ProblemType type) override;
359 :
360 : void SetType(Type type) override;
361 :
362 : void SetInitialSpace(const ComplexVector &v) override;
363 :
364 : int Solve() override;
365 :
366 : std::complex<double> GetEigenvalue(int i) const override;
367 :
368 : void GetEigenvector(int i, ComplexVector &x) const override;
369 :
370 : BV GetBV() const override;
371 :
372 : ST GetST() const override;
373 :
374 : RG GetRG() const override;
375 :
376 0 : MPI_Comm GetComm() const override
377 : {
378 0 : return pep ? PetscObjectComm(reinterpret_cast<PetscObject>(pep)) : MPI_COMM_NULL;
379 : }
380 :
381 0 : operator PetscObject() const override { return reinterpret_cast<PetscObject>(pep); };
382 : };
383 :
384 : // Quadratic eigenvalue problem solver: P(λ) x = (K + λ C + λ² M) x = 0 .
385 : class SlepcPEPSolver : public SlepcPEPSolverBase
386 : {
387 : public:
388 : using SlepcEigenvalueSolver::delta;
389 : using SlepcEigenvalueSolver::gamma;
390 : using SlepcEigenvalueSolver::opB;
391 : using SlepcEigenvalueSolver::opInv;
392 : using SlepcEigenvalueSolver::opProj;
393 : using SlepcEigenvalueSolver::sigma;
394 : using SlepcEigenvalueSolver::sinvert;
395 :
396 : private:
397 : // Operator norms for scaling.
398 : mutable PetscReal normK, normC, normM;
399 :
400 : protected:
401 : PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x,
402 : ComplexVector &r) const override;
403 :
404 : PetscReal GetBackwardScaling(PetscScalar l) const override;
405 :
406 : public:
407 : SlepcPEPSolver(MPI_Comm comm, int print, const std::string &prefix = std::string());
408 :
409 : using SlepcEigenvalueSolver::SetOperators;
410 : void SetOperators(const ComplexOperator &K, const ComplexOperator &C,
411 : const ComplexOperator &M, ScaleType type) override;
412 : void SetBMat(const Operator &B) override;
413 : };
414 :
415 : // Base class for SLEPc's NEP problem type
416 : class SlepcNEPSolverBase : public SlepcEigenvalueSolver
417 : {
418 : protected:
419 : // SLEPc eigensolver object.
420 : NEP nep;
421 :
422 : // Shell matrices for the nonlinear eigenvalue problem.
423 : Mat A, J;
424 :
425 : // Order of sorted eigenvalues.
426 : std::unique_ptr<int[]> perm;
427 :
428 : void Customize() override;
429 :
430 : public:
431 : // Calls SLEPc's NEPCreate. Expects SLEPc to be initialized/finalized externally.
432 : SlepcNEPSolverBase(MPI_Comm comm, int print, const std::string &prefix = std::string());
433 :
434 : // Call's SLEPc's NEPDestroy.
435 : ~SlepcNEPSolverBase() override;
436 :
437 : // Conversion function to SLEPc's PEP type.
438 : operator NEP() const { return nep; }
439 :
440 : void SetNumModes(int num_eig, int num_vec = 0) override;
441 :
442 : void SetTol(PetscReal tol) override;
443 :
444 : void SetMaxIter(int max_it) override;
445 :
446 : void SetWhichEigenpairs(WhichType type) override;
447 :
448 : void SetShiftInvert(std::complex<double> s, bool precond = false) override;
449 :
450 : void SetProblemType(ProblemType type) override;
451 :
452 : void SetType(Type type) override;
453 :
454 : void SetInitialSpace(const ComplexVector &v) override;
455 :
456 : int Solve() override;
457 :
458 : std::complex<double> GetEigenvalue(int i) const override;
459 :
460 : void GetEigenvector(int i, ComplexVector &x) const override;
461 :
462 : BV GetBV() const override;
463 :
464 : ST GetST() const override;
465 :
466 : RG GetRG() const override;
467 :
468 0 : MPI_Comm GetComm() const override
469 : {
470 0 : return nep ? PetscObjectComm(reinterpret_cast<PetscObject>(nep)) : MPI_COMM_NULL;
471 : }
472 :
473 0 : operator PetscObject() const override { return reinterpret_cast<PetscObject>(nep); };
474 : };
475 :
476 : // Nonlinear eigenvalue problem solver: T(λ) x = (K + λ C + λ² M + A2(λ)) x = 0.
477 : class SlepcNEPSolver : public SlepcNEPSolverBase
478 : {
479 : public:
480 : using SlepcEigenvalueSolver::delta;
481 : using SlepcEigenvalueSolver::gamma;
482 : using SlepcEigenvalueSolver::opB;
483 : using SlepcEigenvalueSolver::opInv;
484 : using SlepcEigenvalueSolver::opProj;
485 : using SlepcEigenvalueSolver::sigma;
486 : using SlepcEigenvalueSolver::sinvert;
487 :
488 : // Operators for the nonlinear eigenvalue problem.
489 : std::unique_ptr<ComplexOperator> opA2, opA2p, opJ, opA, opAJ, opA2_pc, opA_pc, opP_pc;
490 :
491 : // Function to compute the A2 operator.
492 : std::optional<std::function<std::unique_ptr<ComplexOperator>(double)>> funcA2;
493 :
494 : // Function to compute the preconditioner matrix.
495 : std::optional<std::function<std::unique_ptr<ComplexOperator>(
496 : std::complex<double>, std::complex<double>, std::complex<double>, double)>>
497 : funcP;
498 :
499 : // Eigenvalue estimate at current iteration.
500 : PetscScalar lambda;
501 :
502 : // Boolean flag to identify new λ estimate requiring a preconditioner update.
503 : bool new_lambda = true;
504 :
505 : // Boolean flag to avoid modifying an unused preconditioner.
506 : bool first_pc = true;
507 :
508 : private:
509 : // Operator norms for scaling.
510 : mutable PetscReal normK, normC, normM;
511 :
512 : protected:
513 : PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x,
514 : ComplexVector &r) const override;
515 :
516 : PetscReal GetBackwardScaling(PetscScalar l) const override;
517 :
518 : public:
519 : SlepcNEPSolver(MPI_Comm comm, int print, const std::string &prefix = std::string());
520 :
521 : using SlepcEigenvalueSolver::SetOperators;
522 : void SetOperators(const ComplexOperator &K, const ComplexOperator &M,
523 : ScaleType type) override;
524 : void SetOperators(const ComplexOperator &K, const ComplexOperator &C,
525 : const ComplexOperator &M, ScaleType type) override;
526 : void SetBMat(const Operator &B) override;
527 :
528 : // Set the frequency-dependent A2 matrix function.
529 : void SetExtraSystemMatrix(
530 : std::function<std::unique_ptr<ComplexOperator>(double)>) override;
531 :
532 : // Set the preconditioner update function.
533 : void SetPreconditionerUpdate(std::function<std::unique_ptr<ComplexOperator>(
534 : std::complex<double>, std::complex<double>,
535 : std::complex<double>, double)>) override;
536 : };
537 :
538 : } // namespace slepc
539 :
540 : } // namespace palace
541 :
542 : #endif
543 :
544 : #endif // PALACE_LINALG_SLEPC_HPP
|
