Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "slepc.hpp"
5 :
6 : #if defined(PALACE_WITH_SLEPC)
7 :
8 : #include <algorithm>
9 : #include <petsc.h>
10 : #include <slepc.h>
11 : #include <mfem.hpp>
12 : #include "linalg/divfree.hpp"
13 : #include "linalg/nleps.hpp"
14 : #include "linalg/rap.hpp"
15 : #include "utils/communication.hpp"
16 :
17 : static PetscErrorCode __mat_apply_EPS_A0(Mat, Vec, Vec);
18 : static PetscErrorCode __mat_apply_EPS_A1(Mat, Vec, Vec);
19 : static PetscErrorCode __mat_apply_EPS_B(Mat, Vec, Vec);
20 : static PetscErrorCode __pc_apply_EPS(PC, Vec, Vec);
21 : static PetscErrorCode __mat_apply_PEPLinear_L0(Mat, Vec, Vec);
22 : static PetscErrorCode __mat_apply_PEPLinear_L1(Mat, Vec, Vec);
23 : static PetscErrorCode __mat_apply_PEPLinear_B(Mat, Vec, Vec);
24 : static PetscErrorCode __pc_apply_PEPLinear(PC, Vec, Vec);
25 : static PetscErrorCode __mat_apply_PEP_A0(Mat, Vec, Vec);
26 : static PetscErrorCode __mat_apply_PEP_A1(Mat, Vec, Vec);
27 : static PetscErrorCode __mat_apply_PEP_A2(Mat, Vec, Vec);
28 : static PetscErrorCode __mat_apply_PEP_B(Mat, Vec, Vec);
29 : static PetscErrorCode __pc_apply_PEP(PC, Vec, Vec);
30 : // for NEP
31 : static PetscErrorCode __mat_apply_NEP_A(Mat, Vec, Vec);
32 : static PetscErrorCode __mat_apply_NEP_J(Mat, Vec, Vec);
33 : static PetscErrorCode __mat_apply_NEP_B(Mat, Vec, Vec);
34 : static PetscErrorCode __pc_apply_NEP(PC, Vec, Vec);
35 : static PetscErrorCode __form_NEP_function(NEP, PetscScalar, Mat, Mat, void *);
36 : static PetscErrorCode __form_NEP_jacobian(NEP, PetscScalar, Mat, void *);
37 :
38 : using namespace std::complex_literals;
39 :
40 : namespace
41 : {
42 :
43 0 : inline PetscErrorCode FromPetscVec(Vec x, palace::ComplexVector &y, int block = 0,
44 : int nblocks = 1)
45 : {
46 : PetscInt n;
47 : const PetscScalar *px;
48 : PetscMemType mtype;
49 0 : PetscCall(VecGetLocalSize(x, &n));
50 : MFEM_ASSERT(y.Size() * nblocks == n,
51 : "Invalid size mismatch for PETSc vector conversion!");
52 0 : PetscCall(VecGetArrayReadAndMemType(x, &px, &mtype));
53 0 : y.Set(px + block * n / nblocks, n / nblocks, PetscMemTypeDevice(mtype));
54 0 : PetscCall(VecRestoreArrayReadAndMemType(x, &px));
55 : return PETSC_SUCCESS;
56 : }
57 :
58 0 : inline PetscErrorCode FromPetscVec(Vec x, palace::ComplexVector &y1,
59 : palace::ComplexVector &y2)
60 : {
61 : PetscInt n;
62 : const PetscScalar *px;
63 : PetscMemType mtype;
64 0 : PetscCall(VecGetLocalSize(x, &n));
65 : MFEM_ASSERT(y1.Size() == n / 2 && y2.Size() == n / 2,
66 : "Invalid size mismatch for PETSc vector conversion!");
67 0 : PetscCall(VecGetArrayReadAndMemType(x, &px, &mtype));
68 0 : y1.Set(px, n / 2, PetscMemTypeDevice(mtype));
69 0 : y2.Set(px + n / 2, n / 2, PetscMemTypeDevice(mtype));
70 0 : PetscCall(VecRestoreArrayReadAndMemType(x, &px));
71 : return PETSC_SUCCESS;
72 : }
73 :
74 0 : inline PetscErrorCode ToPetscVec(const palace::ComplexVector &x, Vec y, int block = 0,
75 : int nblocks = 1)
76 : {
77 : PetscInt n;
78 : PetscScalar *py;
79 : PetscMemType mtype;
80 0 : PetscCall(VecGetLocalSize(y, &n));
81 : MFEM_ASSERT(x.Size() * nblocks == n,
82 : "Invalid size mismatch for PETSc vector conversion!");
83 0 : PetscCall(VecGetArrayWriteAndMemType(y, &py, &mtype));
84 0 : x.Get(py + block * n / nblocks, n / nblocks, PetscMemTypeDevice(mtype));
85 0 : PetscCall(VecRestoreArrayWriteAndMemType(y, &py));
86 : return PETSC_SUCCESS;
87 : }
88 :
89 0 : inline PetscErrorCode ToPetscVec(const palace::ComplexVector &x1,
90 : const palace::ComplexVector &x2, Vec y)
91 : {
92 : PetscInt n;
93 : PetscScalar *py;
94 : PetscMemType mtype;
95 0 : PetscCall(VecGetLocalSize(y, &n));
96 : MFEM_ASSERT(x1.Size() == n / 2 && x2.Size() == n / 2,
97 : "Invalid size mismatch for PETSc vector conversion!");
98 0 : PetscCall(VecGetArrayWriteAndMemType(y, &py, &mtype));
99 0 : x1.Get(py, n / 2, PetscMemTypeDevice(mtype));
100 0 : x2.Get(py + n / 2, n / 2, PetscMemTypeDevice(mtype));
101 0 : PetscCall(VecRestoreArrayWriteAndMemType(y, &py));
102 : return PETSC_SUCCESS;
103 : }
104 :
105 : } // namespace
106 :
107 : namespace palace::slepc
108 : {
109 :
110 : namespace
111 : {
112 :
113 0 : inline PetscErrorCode ConfigurePetscDevice()
114 : {
115 : // Tell PETSc to use the same CUDA or HIP device as MFEM.
116 0 : if (mfem::Device::Allows(mfem::Backend::CUDA_MASK))
117 : {
118 0 : PetscCall(PetscOptionsSetValue(NULL, "-use_gpu_aware_mpi",
119 : mfem::Device::GetGPUAwareMPI() ? "1" : "0"));
120 0 : PetscCall(PetscOptionsSetValue(NULL, "-device_select_cuda",
121 : std::to_string(mfem::Device::GetId()).c_str()));
122 : // PetscCall(PetscOptionsSetValue(NULL, "-bv_type", "svec"));
123 : // PetscCall(PetscOptionsSetValue(NULL, "-vec_type", "cuda"));
124 : }
125 0 : if (mfem::Device::Allows(mfem::Backend::HIP_MASK))
126 : {
127 0 : PetscCall(PetscOptionsSetValue(NULL, "-use_gpu_aware_mpi",
128 : mfem::Device::GetGPUAwareMPI() ? "1" : "0"));
129 0 : PetscCall(PetscOptionsSetValue(NULL, "-device_select_hip",
130 : std::to_string(mfem::Device::GetId()).c_str()));
131 : // PetscCall(PetscOptionsSetValue(NULL, "-bv_type", "svec"));
132 : // PetscCall(PetscOptionsSetValue(NULL, "-vec_type", "hip"));
133 : }
134 : return PETSC_SUCCESS;
135 : }
136 :
137 : inline VecType PetscVecType()
138 : {
139 0 : if (mfem::Device::Allows(mfem::Backend::CUDA_MASK))
140 : {
141 : return VECCUDA;
142 : }
143 0 : if (mfem::Device::Allows(mfem::Backend::HIP_MASK))
144 : {
145 0 : return VECHIP;
146 : }
147 : return VECSTANDARD;
148 : }
149 :
150 : struct MatShellContext
151 : {
152 : const ComplexOperator &A;
153 : ComplexVector &x, &y;
154 : };
155 :
156 0 : PetscErrorCode __mat_apply_shell(Mat A, Vec x, Vec y)
157 : {
158 : PetscFunctionBeginUser;
159 : MatShellContext *ctx;
160 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
161 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
162 :
163 0 : PetscCall(FromPetscVec(x, ctx->x));
164 0 : ctx->A.Mult(ctx->x, ctx->y);
165 0 : PetscCall(ToPetscVec(ctx->y, y));
166 :
167 : PetscFunctionReturn(PETSC_SUCCESS);
168 : }
169 :
170 0 : PetscErrorCode __mat_apply_transpose_shell(Mat A, Vec x, Vec y)
171 : {
172 : PetscFunctionBeginUser;
173 : MatShellContext *ctx;
174 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
175 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
176 :
177 0 : PetscCall(FromPetscVec(x, ctx->x));
178 0 : ctx->A.MultTranspose(ctx->x, ctx->y);
179 0 : PetscCall(ToPetscVec(ctx->y, y));
180 :
181 : PetscFunctionReturn(PETSC_SUCCESS);
182 : }
183 :
184 0 : PetscErrorCode __mat_apply_hermitian_transpose_shell(Mat A, Vec x, Vec y)
185 : {
186 : PetscFunctionBeginUser;
187 : MatShellContext *ctx;
188 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
189 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
190 :
191 0 : PetscCall(FromPetscVec(x, ctx->x));
192 0 : ctx->A.MultHermitianTranspose(ctx->x, ctx->y);
193 0 : PetscCall(ToPetscVec(ctx->y, y));
194 :
195 : PetscFunctionReturn(PETSC_SUCCESS);
196 : };
197 :
198 0 : inline void ConfigurePCShell(ST st, void *ctx, PetscErrorCode (*__pc_apply)(PC, Vec, Vec))
199 : {
200 : KSP ksp;
201 : PC pc;
202 0 : PalacePetscCall(STGetKSP(st, &ksp));
203 0 : PalacePetscCall(KSPGetPC(ksp, &pc));
204 0 : PalacePetscCall(PCSetType(pc, PCSHELL));
205 0 : PalacePetscCall(PCShellSetContext(pc, ctx));
206 0 : PalacePetscCall(PCShellSetApply(pc, __pc_apply));
207 0 : }
208 :
209 0 : inline void ConfigureRG(RG rg, PetscReal lr, PetscReal ur, PetscReal li, PetscReal ui,
210 : bool complement = false)
211 : {
212 0 : PalacePetscCall(RGSetType(rg, RGINTERVAL));
213 0 : PalacePetscCall(RGIntervalSetEndpoints(rg, lr, ur, li, ui));
214 0 : if (complement)
215 : {
216 0 : PalacePetscCall(RGSetComplement(rg, PETSC_TRUE));
217 : }
218 0 : }
219 :
220 : } // namespace
221 :
222 0 : void Initialize(int &argc, char **&argv, const char rc_file[], const char help[])
223 : {
224 0 : ConfigurePetscDevice();
225 0 : PalacePetscCall(SlepcInitialize(&argc, &argv, rc_file, help));
226 :
227 : // Remove default PETSc signal handling, since it can be confusing when the errors occur
228 : // not from within SLEPc/PETSc.
229 0 : PalacePetscCall(PetscPopSignalHandler());
230 0 : }
231 :
232 0 : void Initialize()
233 : {
234 0 : ConfigurePetscDevice();
235 0 : PalacePetscCall(SlepcInitializeNoArguments());
236 :
237 : // Remove default PETSc signal handling, since it can be confusing when the errors occur
238 : // not from within SLEPc/PETSc.
239 0 : PalacePetscCall(PetscPopSignalHandler());
240 0 : }
241 :
242 0 : void Finalize()
243 : {
244 0 : PalacePetscCall(SlepcFinalize());
245 0 : }
246 :
247 0 : PetscReal GetMaxSingularValue(MPI_Comm comm, const ComplexOperator &A, bool herm,
248 : PetscReal tol, PetscInt max_it)
249 : {
250 : // This method assumes the provided operator has the required operations for SLEPc's EPS
251 : // or SVD solvers, namely MATOP_MULT and MATOP_MULT_HERMITIAN_TRANSPOSE (if the matrix
252 : // is not Hermitian).
253 0 : MFEM_VERIFY(A.Height() == A.Width(), "Spectral norm requires a square matrix!");
254 0 : const PetscInt n = A.Height();
255 0 : ComplexVector x(n), y(n);
256 0 : x.UseDevice(true);
257 0 : y.UseDevice(true);
258 0 : MatShellContext ctx = {A, x, y};
259 : Mat A0;
260 0 : PalacePetscCall(
261 : MatCreateShell(comm, n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)&ctx, &A0));
262 0 : PalacePetscCall(MatShellSetOperation(A0, MATOP_MULT, (void (*)(void))__mat_apply_shell));
263 0 : PalacePetscCall(MatShellSetVecType(A0, PetscVecType()));
264 0 : if (herm)
265 : {
266 : EPS eps;
267 : PetscInt num_conv;
268 0 : PetscScalar eig;
269 0 : PalacePetscCall(EPSCreate(comm, &eps));
270 0 : PalacePetscCall(EPSSetOperators(eps, A0, nullptr));
271 0 : PalacePetscCall(EPSSetProblemType(eps, EPS_HEP));
272 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE));
273 0 : PalacePetscCall(EPSSetDimensions(eps, 1, PETSC_DEFAULT, PETSC_DEFAULT));
274 0 : PalacePetscCall(EPSSetTolerances(eps, tol, max_it));
275 0 : PalacePetscCall(EPSSolve(eps));
276 0 : PalacePetscCall(EPSGetConverged(eps, &num_conv));
277 0 : if (num_conv < 1)
278 : {
279 0 : Mpi::Warning(comm, "SLEPc EPS solve did not converge for maximum singular value!\n");
280 : eig = 0.0;
281 : }
282 : else
283 : {
284 0 : PalacePetscCall(EPSGetEigenvalue(eps, 0, &eig, nullptr));
285 0 : MFEM_VERIFY(PetscImaginaryPart(eig) == 0.0,
286 : "Unexpected complex eigenvalue for Hermitian matrix (λ = " << eig
287 : << ")!");
288 : }
289 0 : PalacePetscCall(EPSDestroy(&eps));
290 0 : PalacePetscCall(MatDestroy(&A0));
291 0 : return PetscAbsScalar(eig);
292 : }
293 : else
294 : {
295 0 : PalacePetscCall(MatShellSetOperation(A0, MATOP_MULT_TRANSPOSE,
296 : (void (*)(void))__mat_apply_transpose_shell));
297 0 : PalacePetscCall(
298 : MatShellSetOperation(A0, MATOP_MULT_HERMITIAN_TRANSPOSE,
299 : (void (*)(void))__mat_apply_hermitian_transpose_shell));
300 : SVD svd;
301 : PetscInt num_conv;
302 : PetscReal sigma;
303 0 : PalacePetscCall(SVDCreate(comm, &svd));
304 0 : PalacePetscCall(SVDSetOperators(svd, A0, nullptr));
305 0 : PalacePetscCall(SVDSetProblemType(svd, SVD_STANDARD));
306 0 : PalacePetscCall(SVDSetWhichSingularTriplets(svd, SVD_LARGEST));
307 0 : PalacePetscCall(SVDSetDimensions(svd, 1, PETSC_DEFAULT, PETSC_DEFAULT));
308 0 : PalacePetscCall(SVDSetTolerances(svd, tol, max_it));
309 0 : PalacePetscCall(SVDSolve(svd));
310 0 : PalacePetscCall(SVDGetConverged(svd, &num_conv));
311 0 : if (num_conv < 1)
312 : {
313 0 : Mpi::Warning(comm, "SLEPc SVD solve did not converge for maximum singular value!\n");
314 0 : sigma = 0.0;
315 : }
316 : else
317 : {
318 0 : PalacePetscCall(SVDGetSingularTriplet(svd, 0, &sigma, nullptr, nullptr));
319 : }
320 0 : PalacePetscCall(SVDDestroy(&svd));
321 0 : PalacePetscCall(MatDestroy(&A0));
322 0 : return sigma;
323 : }
324 : }
325 :
326 : // Eigensolver base class methods.
327 :
328 0 : SlepcEigenvalueSolver::SlepcEigenvalueSolver(int print) : print(print)
329 : {
330 0 : sinvert = false;
331 0 : region = true;
332 : sigma = 0.0;
333 0 : gamma = delta = 1.0;
334 :
335 0 : opInv = nullptr;
336 0 : opProj = nullptr;
337 0 : opB = nullptr;
338 :
339 0 : B0 = nullptr;
340 0 : v0 = nullptr;
341 :
342 0 : cl_custom = false;
343 0 : }
344 :
345 0 : SlepcEigenvalueSolver::~SlepcEigenvalueSolver()
346 : {
347 0 : PalacePetscCall(MatDestroy(&B0));
348 0 : PalacePetscCall(VecDestroy(&v0));
349 0 : }
350 :
351 0 : void SlepcEigenvalueSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &M,
352 : EigenvalueSolver::ScaleType type)
353 : {
354 0 : MFEM_ABORT("SetOperators not defined for base class SlepcEigenvalueSolver!");
355 : }
356 :
357 0 : void SlepcEigenvalueSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &C,
358 : const ComplexOperator &M,
359 : EigenvalueSolver::ScaleType type)
360 : {
361 0 : MFEM_ABORT("SetOperators not defined for base class SlepcEigenvalueSolver!");
362 : }
363 :
364 0 : void SlepcEigenvalueSolver::SetLinearSolver(ComplexKspSolver &ksp)
365 : {
366 0 : opInv = &ksp;
367 0 : }
368 :
369 0 : void SlepcEigenvalueSolver::SetDivFreeProjector(const DivFreeSolver<ComplexVector> &divfree)
370 : {
371 0 : opProj = &divfree;
372 0 : }
373 :
374 0 : void SlepcEigenvalueSolver::SetBMat(const Operator &B)
375 : {
376 0 : opB = &B;
377 0 : }
378 :
379 0 : void SlepcEigenvalueSolver::SetShiftInvert(std::complex<double> s, bool precond)
380 : {
381 0 : ST st = GetST();
382 0 : if (precond)
383 : {
384 0 : PalacePetscCall(STSetType(st, STPRECOND));
385 : }
386 : else
387 : {
388 0 : PalacePetscCall(STSetType(st, STSINVERT));
389 : }
390 0 : PalacePetscCall(STSetTransform(st, PETSC_TRUE));
391 0 : PalacePetscCall(STSetMatMode(st, ST_MATMODE_SHELL));
392 0 : sigma = s; // Wait until solve time to call EPS/PEPSetTarget
393 0 : sinvert = true;
394 0 : }
395 :
396 0 : void SlepcEigenvalueSolver::SetOrthogonalization(bool mgs, bool cgs2)
397 : {
398 : // The SLEPc default is CGS with refinement if needed.
399 0 : if (mgs || cgs2)
400 : {
401 0 : BV bv = GetBV();
402 : BVOrthogType type;
403 : BVOrthogRefineType refine;
404 0 : if (mgs)
405 : {
406 : type = BV_ORTHOG_MGS;
407 : refine = BV_ORTHOG_REFINE_NEVER;
408 : }
409 : else // cgs2
410 : {
411 : type = BV_ORTHOG_CGS;
412 : refine = BV_ORTHOG_REFINE_ALWAYS;
413 : }
414 0 : PalacePetscCall(BVSetOrthogonalization(bv, type, refine, 1.0, BV_ORTHOG_BLOCK_GS));
415 : }
416 0 : }
417 :
418 0 : void SlepcEigenvalueSolver::Customize()
419 : {
420 : // Configure the KSP object for non-preconditioned spectral transformations.
421 : PetscBool precond;
422 0 : ST st = GetST();
423 0 : PalacePetscCall(
424 : PetscObjectTypeCompare(reinterpret_cast<PetscObject>(st), STPRECOND, &precond));
425 0 : if (!precond)
426 : {
427 : KSP ksp;
428 0 : PalacePetscCall(STGetKSP(st, &ksp));
429 0 : PalacePetscCall(KSPSetType(ksp, KSPPREONLY));
430 : }
431 :
432 : // Configure the region based on the given target if necessary.
433 0 : if (sinvert && region)
434 : {
435 0 : if (PetscImaginaryPart(sigma) == 0.0)
436 : {
437 : PetscReal sr = PetscRealPart(sigma);
438 0 : if (sr > 0.0)
439 : {
440 0 : ConfigureRG(GetRG(), sr / gamma, mfem::infinity(), -mfem::infinity(),
441 : mfem::infinity());
442 : }
443 0 : else if (sr < 0.0)
444 : {
445 0 : ConfigureRG(GetRG(), -mfem::infinity(), sr / gamma, -mfem::infinity(),
446 : mfem::infinity());
447 : }
448 : }
449 0 : else if (PetscRealPart(sigma) == 0.0)
450 : {
451 : PetscReal si = PetscImaginaryPart(sigma);
452 0 : if (si > 0.0)
453 : {
454 0 : ConfigureRG(GetRG(), -mfem::infinity(), mfem::infinity(), si / gamma,
455 : mfem::infinity());
456 : }
457 0 : else if (si < 0.0)
458 : {
459 0 : ConfigureRG(GetRG(), -mfem::infinity(), mfem::infinity(), -mfem::infinity(),
460 0 : si / gamma);
461 : }
462 : }
463 : else
464 : {
465 0 : MFEM_ABORT("Shift-and-invert with general complex eigenvalue target is unsupported!");
466 : }
467 : }
468 0 : }
469 :
470 0 : PetscReal SlepcEigenvalueSolver::GetEigenvectorNorm(const ComplexVector &x,
471 : ComplexVector &Bx) const
472 : {
473 0 : if (opB)
474 : {
475 0 : return linalg::Norml2(GetComm(), x, *opB, Bx);
476 : }
477 : else
478 : {
479 0 : return linalg::Norml2(GetComm(), x);
480 : }
481 : }
482 :
483 0 : PetscReal SlepcEigenvalueSolver::GetError(int i, EigenvalueSolver::ErrorType type) const
484 : {
485 0 : switch (type)
486 : {
487 0 : case ErrorType::ABSOLUTE:
488 0 : return res.get()[i];
489 0 : case ErrorType::RELATIVE:
490 0 : return res.get()[i] / PetscAbsScalar(GetEigenvalue(i));
491 0 : case ErrorType::BACKWARD:
492 0 : return res.get()[i] / GetBackwardScaling(GetEigenvalue(i));
493 : }
494 : return 0.0;
495 : }
496 :
497 0 : void SlepcEigenvalueSolver::RescaleEigenvectors(int num_eig)
498 : {
499 0 : res = std::make_unique<PetscReal[]>(num_eig);
500 0 : xscale = std::make_unique<PetscReal[]>(num_eig);
501 0 : for (int i = 0; i < num_eig; i++)
502 : {
503 0 : xscale.get()[i] = 0.0;
504 0 : GetEigenvector(i, x1);
505 0 : xscale.get()[i] = 1.0 / GetEigenvectorNorm(x1, y1);
506 0 : res.get()[i] =
507 0 : GetResidualNorm(GetEigenvalue(i), x1, y1) / linalg::Norml2(GetComm(), x1);
508 : }
509 0 : }
510 :
511 : // EPS specific methods.
512 :
513 0 : SlepcEPSSolverBase::SlepcEPSSolverBase(MPI_Comm comm, int print, const std::string &prefix)
514 0 : : SlepcEigenvalueSolver(print)
515 : {
516 0 : PalacePetscCall(EPSCreate(comm, &eps));
517 0 : PalacePetscCall(EPSSetOptionsPrefix(eps, prefix.c_str()));
518 0 : if (print > 0)
519 : {
520 0 : std::string opts = "-eps_monitor";
521 0 : if (print > 2)
522 : {
523 0 : opts.append(" -eps_view");
524 : }
525 0 : if (prefix.length() > 0)
526 : {
527 0 : PetscOptionsPrefixPush(nullptr, prefix.c_str());
528 : }
529 0 : PetscOptionsInsertString(nullptr, opts.c_str());
530 0 : if (prefix.length() > 0)
531 : {
532 0 : PetscOptionsPrefixPop(nullptr);
533 : }
534 : }
535 0 : A0 = A1 = nullptr;
536 0 : }
537 :
538 0 : SlepcEPSSolverBase::~SlepcEPSSolverBase()
539 : {
540 0 : PalacePetscCall(EPSDestroy(&eps));
541 0 : PalacePetscCall(MatDestroy(&A0));
542 0 : PalacePetscCall(MatDestroy(&A1));
543 0 : }
544 :
545 0 : void SlepcEPSSolverBase::SetNumModes(int num_eig, int num_vec)
546 : {
547 0 : PalacePetscCall(EPSSetDimensions(eps, num_eig, (num_vec > 0) ? num_vec : PETSC_DEFAULT,
548 : PETSC_DEFAULT));
549 0 : }
550 :
551 0 : void SlepcEPSSolverBase::SetTol(PetscReal tol)
552 : {
553 0 : PalacePetscCall(EPSSetTolerances(eps, tol, PETSC_DEFAULT));
554 0 : PalacePetscCall(EPSSetConvergenceTest(eps, EPS_CONV_REL));
555 : // PalacePetscCall(EPSSetTrackAll(eps, PETSC_TRUE));
556 : // PalacePetscCall(EPSSetTrueResidual(eps, PETSC_TRUE));
557 0 : }
558 :
559 0 : void SlepcEPSSolverBase::SetMaxIter(int max_it)
560 : {
561 0 : PalacePetscCall(
562 : EPSSetTolerances(eps, PETSC_DEFAULT, (max_it > 0) ? max_it : PETSC_DEFAULT));
563 0 : }
564 :
565 0 : void SlepcEPSSolverBase::SetWhichEigenpairs(EigenvalueSolver::WhichType type)
566 : {
567 0 : switch (type)
568 : {
569 0 : case WhichType::LARGEST_MAGNITUDE:
570 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE));
571 0 : region = false;
572 0 : break;
573 0 : case WhichType::SMALLEST_MAGNITUDE:
574 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE));
575 0 : region = false;
576 0 : break;
577 0 : case WhichType::LARGEST_REAL:
578 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_LARGEST_REAL));
579 : break;
580 0 : case WhichType::SMALLEST_REAL:
581 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL));
582 : break;
583 0 : case WhichType::LARGEST_IMAGINARY:
584 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_LARGEST_IMAGINARY));
585 : break;
586 0 : case WhichType::SMALLEST_IMAGINARY:
587 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_IMAGINARY));
588 : break;
589 0 : case WhichType::TARGET_MAGNITUDE:
590 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE));
591 0 : region = false;
592 0 : break;
593 0 : case WhichType::TARGET_REAL:
594 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_TARGET_REAL));
595 : break;
596 0 : case WhichType::TARGET_IMAGINARY:
597 0 : PalacePetscCall(EPSSetWhichEigenpairs(eps, EPS_TARGET_IMAGINARY));
598 : break;
599 : }
600 0 : }
601 :
602 0 : void SlepcEPSSolverBase::SetProblemType(SlepcEigenvalueSolver::ProblemType type)
603 : {
604 0 : switch (type)
605 : {
606 0 : case ProblemType::HERMITIAN:
607 0 : PalacePetscCall(EPSSetProblemType(eps, EPS_HEP));
608 : break;
609 0 : case ProblemType::NON_HERMITIAN:
610 0 : PalacePetscCall(EPSSetProblemType(eps, EPS_NHEP));
611 : break;
612 0 : case ProblemType::GEN_HERMITIAN:
613 0 : PalacePetscCall(EPSSetProblemType(eps, EPS_GHEP));
614 : break;
615 0 : case ProblemType::GEN_INDEFINITE:
616 0 : PalacePetscCall(EPSSetProblemType(eps, EPS_GHIEP));
617 : break;
618 0 : case ProblemType::GEN_NON_HERMITIAN:
619 0 : PalacePetscCall(EPSSetProblemType(eps, EPS_GNHEP));
620 : // PalacePetscCall(EPSSetProblemType(eps, EPS_PGNHEP)); // If B is SPD
621 : break;
622 0 : case ProblemType::HYPERBOLIC:
623 : case ProblemType::GYROSCOPIC:
624 : case ProblemType::GENERAL:
625 0 : MFEM_ABORT("Problem type not implemented!");
626 : break;
627 : }
628 0 : }
629 :
630 0 : void SlepcEPSSolverBase::SetType(SlepcEigenvalueSolver::Type type)
631 : {
632 0 : switch (type)
633 : {
634 0 : case Type::KRYLOVSCHUR:
635 0 : PalacePetscCall(EPSSetType(eps, EPSKRYLOVSCHUR));
636 : break;
637 0 : case Type::POWER:
638 0 : PalacePetscCall(EPSSetType(eps, EPSPOWER));
639 : break;
640 0 : case Type::SUBSPACE:
641 0 : PalacePetscCall(EPSSetType(eps, EPSSUBSPACE));
642 : break;
643 0 : case Type::JD:
644 0 : PalacePetscCall(EPSSetType(eps, EPSJD));
645 0 : region = false;
646 0 : break;
647 0 : case Type::TOAR:
648 : case Type::STOAR:
649 : case Type::QARNOLDI:
650 : case Type::SLP:
651 : case Type::NLEIGS:
652 0 : MFEM_ABORT("Eigenvalue solver type not implemented!");
653 : break;
654 : }
655 0 : }
656 :
657 0 : void SlepcEPSSolverBase::SetInitialSpace(const ComplexVector &v)
658 : {
659 0 : MFEM_VERIFY(
660 : A0 && A1,
661 : "Must call SetOperators before using SetInitialSpace for SLEPc eigenvalue solver!");
662 0 : if (!v0)
663 : {
664 0 : PalacePetscCall(MatCreateVecs(A0, nullptr, &v0));
665 : }
666 0 : PalacePetscCall(ToPetscVec(v, v0));
667 0 : Vec is[1] = {v0};
668 0 : PalacePetscCall(EPSSetInitialSpace(eps, 1, is));
669 0 : }
670 :
671 0 : void SlepcEPSSolverBase::Customize()
672 : {
673 0 : SlepcEigenvalueSolver::Customize();
674 0 : PalacePetscCall(EPSSetTarget(eps, sigma / gamma));
675 0 : if (!cl_custom)
676 : {
677 0 : PalacePetscCall(EPSSetFromOptions(eps));
678 0 : if (print > 0)
679 : {
680 0 : PetscOptionsView(nullptr, PETSC_VIEWER_STDOUT_(GetComm()));
681 0 : Mpi::Print(GetComm(), "\n");
682 : }
683 0 : cl_custom = true;
684 : }
685 0 : }
686 :
687 0 : int SlepcEPSSolverBase::Solve()
688 : {
689 0 : MFEM_VERIFY(A0 && A1 && opInv, "Operators are not set for SlepcEPSSolverBase!");
690 :
691 : // Solve the eigenvalue problem.
692 : PetscInt num_conv;
693 0 : Customize();
694 0 : PalacePetscCall(EPSSolve(eps));
695 0 : PalacePetscCall(EPSGetConverged(eps, &num_conv));
696 0 : if (print > 0)
697 : {
698 0 : Mpi::Print(GetComm(), "\n");
699 0 : PalacePetscCall(EPSConvergedReasonView(eps, PETSC_VIEWER_STDOUT_(GetComm())));
700 0 : Mpi::Print(GetComm(),
701 : " Total number of linear systems solved: {:d}\n"
702 : " Total number of linear solver iterations: {:d}\n",
703 0 : opInv->NumTotalMult(), opInv->NumTotalMultIterations());
704 : }
705 :
706 : // Compute and store the eigenpair residuals.
707 0 : RescaleEigenvectors(num_conv);
708 0 : return (int)num_conv;
709 : }
710 :
711 0 : std::complex<double> SlepcEPSSolverBase::GetEigenvalue(int i) const
712 : {
713 0 : PetscScalar l;
714 0 : PalacePetscCall(EPSGetEigenvalue(eps, i, &l, nullptr));
715 0 : return l * gamma;
716 : }
717 :
718 0 : void SlepcEPSSolverBase::GetEigenvector(int i, ComplexVector &x) const
719 : {
720 0 : MFEM_VERIFY(
721 : v0,
722 : "Must call SetOperators before using GetEigenvector for SLEPc eigenvalue solver!");
723 0 : PalacePetscCall(EPSGetEigenvector(eps, i, v0, nullptr));
724 0 : PalacePetscCall(FromPetscVec(v0, x));
725 0 : if (xscale.get()[i] > 0.0)
726 : {
727 0 : x *= xscale.get()[i];
728 : }
729 0 : }
730 :
731 0 : BV SlepcEPSSolverBase::GetBV() const
732 : {
733 : BV bv;
734 0 : PalacePetscCall(EPSGetBV(eps, &bv));
735 0 : return bv;
736 : }
737 :
738 0 : ST SlepcEPSSolverBase::GetST() const
739 : {
740 : ST st;
741 0 : PalacePetscCall(EPSGetST(eps, &st));
742 0 : return st;
743 : }
744 :
745 0 : RG SlepcEPSSolverBase::GetRG() const
746 : {
747 : RG rg;
748 0 : PalacePetscCall(EPSGetRG(eps, &rg));
749 0 : return rg;
750 : }
751 :
752 0 : SlepcEPSSolver::SlepcEPSSolver(MPI_Comm comm, int print, const std::string &prefix)
753 0 : : SlepcEPSSolverBase(comm, print, prefix)
754 : {
755 0 : opK = opM = nullptr;
756 0 : normK = normM = 0.0;
757 0 : }
758 :
759 0 : void SlepcEPSSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &M,
760 : EigenvalueSolver::ScaleType type)
761 : {
762 : // Construct shell matrices for the scaled operators which define the generalized
763 : // eigenvalue problem.
764 0 : const bool first = (opK == nullptr);
765 0 : opK = &K;
766 0 : opM = &M;
767 :
768 0 : if (first)
769 : {
770 0 : const PetscInt n = opK->Height();
771 0 : PalacePetscCall(
772 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &A0));
773 0 : PalacePetscCall(
774 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &A1));
775 0 : PalacePetscCall(
776 : MatShellSetOperation(A0, MATOP_MULT, (void (*)(void))__mat_apply_EPS_A0));
777 0 : PalacePetscCall(
778 : MatShellSetOperation(A1, MATOP_MULT, (void (*)(void))__mat_apply_EPS_A1));
779 0 : PalacePetscCall(MatShellSetVecType(A0, PetscVecType()));
780 0 : PalacePetscCall(MatShellSetVecType(A1, PetscVecType()));
781 0 : PalacePetscCall(EPSSetOperators(eps, A0, A1));
782 : }
783 :
784 0 : if (first && type != ScaleType::NONE)
785 : {
786 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
787 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
788 0 : MFEM_VERIFY(normK >= 0.0 && normM >= 0.0, "Invalid matrix norms for EPS scaling!");
789 0 : if (normK > 0 && normM > 0.0)
790 : {
791 0 : gamma = normK / normM; // Store γ² for linear problem
792 0 : delta = 2.0 / normK;
793 : }
794 : }
795 :
796 : // Set up workspace.
797 0 : if (!v0)
798 : {
799 0 : PalacePetscCall(MatCreateVecs(A0, nullptr, &v0));
800 : }
801 0 : x1.SetSize(opK->Height());
802 0 : y1.SetSize(opK->Height());
803 0 : x1.UseDevice(true);
804 0 : y1.UseDevice(true);
805 :
806 : // Configure linear solver for generalized problem or spectral transformation. This also
807 : // allows use of the divergence-free projector as a linear solve side-effect.
808 0 : if (first)
809 : {
810 0 : ConfigurePCShell(GetST(), (void *)this, __pc_apply_EPS);
811 : }
812 0 : }
813 :
814 0 : void SlepcEPSSolver::SetBMat(const Operator &B)
815 : {
816 0 : SlepcEigenvalueSolver::SetBMat(B);
817 :
818 0 : const PetscInt n = B.Height();
819 0 : PalacePetscCall(
820 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &B0));
821 0 : PalacePetscCall(MatShellSetOperation(B0, MATOP_MULT, (void (*)(void))__mat_apply_EPS_B));
822 0 : PalacePetscCall(MatShellSetVecType(B0, PetscVecType()));
823 :
824 0 : BV bv = GetBV();
825 0 : PalacePetscCall(BVSetMatrix(bv, B0, PETSC_FALSE));
826 0 : }
827 :
828 0 : PetscReal SlepcEPSSolver::GetResidualNorm(PetscScalar l, const ComplexVector &x,
829 : ComplexVector &r) const
830 : {
831 : // Compute the i-th eigenpair residual: || (K - λ M) x ||₂ for eigenvalue λ.
832 0 : opK->Mult(x, r);
833 0 : opM->AddMult(x, r, -l);
834 0 : return linalg::Norml2(GetComm(), r);
835 : }
836 :
837 0 : PetscReal SlepcEPSSolver::GetBackwardScaling(PetscScalar l) const
838 : {
839 : // Make sure not to use norms from scaling as this can be confusing if they are different.
840 : // Note that SLEPc typically uses ||.||∞, not the 2-norm.
841 0 : if (normK <= 0.0)
842 : {
843 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
844 : }
845 0 : if (normM <= 0.0)
846 : {
847 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
848 : }
849 0 : return normK + PetscAbsScalar(l) * normM;
850 : }
851 :
852 0 : SlepcPEPLinearSolver::SlepcPEPLinearSolver(MPI_Comm comm, int print,
853 0 : const std::string &prefix)
854 0 : : SlepcEPSSolverBase(comm, print, prefix)
855 : {
856 0 : opK = opC = opM = nullptr;
857 0 : normK = normC = normM = 0.0;
858 0 : }
859 :
860 0 : void SlepcPEPLinearSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &C,
861 : const ComplexOperator &M,
862 : EigenvalueSolver::ScaleType type)
863 : {
864 : // Construct shell matrices for the scaled linearized operators which define the block 2x2
865 : // eigenvalue problem.
866 0 : const bool first = (opK == nullptr);
867 0 : opK = &K;
868 0 : opC = &C;
869 0 : opM = &M;
870 0 : if (first)
871 : {
872 0 : const PetscInt n = opK->Height();
873 0 : PalacePetscCall(MatCreateShell(GetComm(), 2 * n, 2 * n, PETSC_DECIDE, PETSC_DECIDE,
874 : (void *)this, &A0));
875 0 : PalacePetscCall(MatCreateShell(GetComm(), 2 * n, 2 * n, PETSC_DECIDE, PETSC_DECIDE,
876 : (void *)this, &A1));
877 0 : PalacePetscCall(
878 : MatShellSetOperation(A0, MATOP_MULT, (void (*)(void))__mat_apply_PEPLinear_L0));
879 0 : PalacePetscCall(
880 : MatShellSetOperation(A1, MATOP_MULT, (void (*)(void))__mat_apply_PEPLinear_L1));
881 0 : PalacePetscCall(MatShellSetVecType(A0, PetscVecType()));
882 0 : PalacePetscCall(MatShellSetVecType(A1, PetscVecType()));
883 0 : PalacePetscCall(EPSSetOperators(eps, A0, A1));
884 : }
885 :
886 0 : if (first && type != ScaleType::NONE)
887 : {
888 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
889 0 : normC = linalg::SpectralNorm(GetComm(), *opC, opC->IsReal());
890 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
891 0 : MFEM_VERIFY(normK >= 0.0 && normC >= 0.0 && normM >= 0.0,
892 : "Invalid matrix norms for PEP scaling!");
893 0 : if (normK > 0 && normC >= 0.0 && normM > 0.0)
894 : {
895 0 : gamma = std::sqrt(normK / normM);
896 0 : delta = 2.0 / (normK + gamma * normC);
897 : }
898 : }
899 :
900 : // Set up workspace.
901 0 : if (!v0)
902 : {
903 0 : PalacePetscCall(MatCreateVecs(A0, nullptr, &v0));
904 : }
905 0 : x1.SetSize(opK->Height());
906 0 : x2.SetSize(opK->Height());
907 0 : y1.SetSize(opK->Height());
908 0 : y2.SetSize(opK->Height());
909 0 : x1.UseDevice(true);
910 0 : x2.UseDevice(true);
911 0 : y1.UseDevice(true);
912 0 : y2.UseDevice(true);
913 :
914 : // Configure linear solver.
915 0 : if (first)
916 : {
917 0 : ConfigurePCShell(GetST(), (void *)this, __pc_apply_PEPLinear);
918 : }
919 0 : }
920 :
921 0 : void SlepcPEPLinearSolver::SetBMat(const Operator &B)
922 : {
923 0 : SlepcEigenvalueSolver::SetBMat(B);
924 :
925 0 : const PetscInt n = B.Height();
926 0 : PalacePetscCall(MatCreateShell(GetComm(), 2 * n, 2 * n, PETSC_DECIDE, PETSC_DECIDE,
927 : (void *)this, &B0));
928 0 : PalacePetscCall(
929 : MatShellSetOperation(B0, MATOP_MULT, (void (*)(void))__mat_apply_PEPLinear_B));
930 0 : PalacePetscCall(MatShellSetVecType(B0, PetscVecType()));
931 :
932 0 : BV bv = GetBV();
933 0 : PalacePetscCall(BVSetMatrix(bv, B0, PETSC_FALSE));
934 0 : }
935 :
936 0 : void SlepcPEPLinearSolver::SetInitialSpace(const ComplexVector &v)
937 : {
938 0 : MFEM_VERIFY(
939 : A0 && A1,
940 : "Must call SetOperators before using SetInitialSpace for SLEPc eigenvalue solver!");
941 0 : if (!v0)
942 : {
943 0 : PalacePetscCall(MatCreateVecs(A0, nullptr, &v0));
944 : }
945 0 : PalacePetscCall(ToPetscVec(v, v0, 0, 2));
946 0 : Vec is[1] = {v0};
947 0 : PalacePetscCall(EPSSetInitialSpace(eps, 1, is));
948 0 : }
949 :
950 0 : void SlepcPEPLinearSolver::GetEigenvector(int i, ComplexVector &x) const
951 : {
952 : // Select the most accurate x for y = [x₁; x₂] from the linearized eigenvalue problem. Or,
953 : // just take x = x₁.
954 0 : MFEM_VERIFY(
955 : v0,
956 : "Must call SetOperators before using GetEigenvector for SLEPc eigenvalue solver!");
957 0 : PalacePetscCall(EPSGetEigenvector(eps, i, v0, nullptr));
958 0 : PalacePetscCall(FromPetscVec(v0, x, 0, 2));
959 0 : if (xscale.get()[i] > 0.0)
960 : {
961 0 : x *= xscale.get()[i];
962 : }
963 0 : }
964 :
965 0 : PetscReal SlepcPEPLinearSolver::GetResidualNorm(PetscScalar l, const ComplexVector &x,
966 : ComplexVector &r) const
967 : {
968 : // Compute the i-th eigenpair residual: || P(λ) x ||₂ = || (K + λ C + λ² M) x ||₂ for
969 : // eigenvalue λ.
970 0 : opK->Mult(x, r);
971 0 : if (opC)
972 : {
973 0 : opC->AddMult(x, r, l);
974 : }
975 0 : opM->AddMult(x, r, l * l);
976 0 : return linalg::Norml2(GetComm(), r);
977 : }
978 :
979 0 : PetscReal SlepcPEPLinearSolver::GetBackwardScaling(PetscScalar l) const
980 : {
981 : // Make sure not to use norms from scaling as this can be confusing if they are different.
982 : // Note that SLEPc typically uses ||.||∞, not the 2-norm.
983 0 : if (normK <= 0.0)
984 : {
985 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
986 : }
987 0 : if (normC <= 0.0 && opC)
988 : {
989 0 : normC = linalg::SpectralNorm(GetComm(), *opC, opC->IsReal());
990 : }
991 0 : if (normM <= 0.0)
992 : {
993 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
994 : }
995 0 : PetscReal t = PetscAbsScalar(l);
996 0 : return normK + t * normC + t * t * normM;
997 : }
998 :
999 : // PEP specific methods.
1000 :
1001 0 : SlepcPEPSolverBase::SlepcPEPSolverBase(MPI_Comm comm, int print, const std::string &prefix)
1002 0 : : SlepcEigenvalueSolver(print)
1003 : {
1004 0 : PalacePetscCall(PEPCreate(comm, &pep));
1005 0 : PalacePetscCall(PEPSetOptionsPrefix(pep, prefix.c_str()));
1006 0 : if (print > 0)
1007 : {
1008 0 : std::string opts = "-pep_monitor";
1009 0 : if (print > 2)
1010 : {
1011 0 : opts.append(" -pep_view");
1012 : }
1013 0 : if (prefix.length() > 0)
1014 : {
1015 0 : PetscOptionsPrefixPush(nullptr, prefix.c_str());
1016 : }
1017 0 : PetscOptionsInsertString(nullptr, opts.c_str());
1018 0 : if (prefix.length() > 0)
1019 : {
1020 0 : PetscOptionsPrefixPop(nullptr);
1021 : }
1022 : }
1023 0 : A0 = A1 = A2 = nullptr;
1024 0 : }
1025 :
1026 0 : SlepcPEPSolverBase::~SlepcPEPSolverBase()
1027 : {
1028 0 : PalacePetscCall(PEPDestroy(&pep));
1029 0 : PalacePetscCall(MatDestroy(&A0));
1030 0 : PalacePetscCall(MatDestroy(&A1));
1031 0 : PalacePetscCall(MatDestroy(&A2));
1032 0 : }
1033 :
1034 0 : void SlepcPEPSolverBase::SetNumModes(int num_eig, int num_vec)
1035 : {
1036 0 : PalacePetscCall(PEPSetDimensions(pep, num_eig, (num_vec > 0) ? num_vec : PETSC_DEFAULT,
1037 : PETSC_DEFAULT));
1038 0 : }
1039 :
1040 0 : void SlepcPEPSolverBase::SetTol(PetscReal tol)
1041 : {
1042 0 : PalacePetscCall(PEPSetTolerances(pep, tol, PETSC_DEFAULT));
1043 0 : PalacePetscCall(PEPSetConvergenceTest(pep, PEP_CONV_REL));
1044 : // PalacePetscCall(PEPSetTrackAll(pep, PETSC_TRUE));
1045 0 : }
1046 :
1047 0 : void SlepcPEPSolverBase::SetMaxIter(int max_it)
1048 : {
1049 0 : PalacePetscCall(
1050 : PEPSetTolerances(pep, PETSC_DEFAULT, (max_it > 0) ? max_it : PETSC_DEFAULT));
1051 0 : }
1052 :
1053 0 : void SlepcPEPSolverBase::SetWhichEigenpairs(EigenvalueSolver::WhichType type)
1054 : {
1055 0 : switch (type)
1056 : {
1057 0 : case WhichType::LARGEST_MAGNITUDE:
1058 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_LARGEST_MAGNITUDE));
1059 0 : region = false;
1060 0 : break;
1061 0 : case WhichType::SMALLEST_MAGNITUDE:
1062 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_SMALLEST_MAGNITUDE));
1063 0 : region = false;
1064 0 : break;
1065 0 : case WhichType::LARGEST_REAL:
1066 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_LARGEST_REAL));
1067 : break;
1068 0 : case WhichType::SMALLEST_REAL:
1069 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_SMALLEST_REAL));
1070 : break;
1071 0 : case WhichType::LARGEST_IMAGINARY:
1072 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_LARGEST_IMAGINARY));
1073 : break;
1074 0 : case WhichType::SMALLEST_IMAGINARY:
1075 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_SMALLEST_IMAGINARY));
1076 : break;
1077 0 : case WhichType::TARGET_MAGNITUDE:
1078 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE));
1079 0 : region = false;
1080 0 : break;
1081 0 : case WhichType::TARGET_REAL:
1082 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_TARGET_REAL));
1083 : break;
1084 0 : case WhichType::TARGET_IMAGINARY:
1085 0 : PalacePetscCall(PEPSetWhichEigenpairs(pep, PEP_TARGET_IMAGINARY));
1086 : break;
1087 : }
1088 0 : }
1089 :
1090 0 : void SlepcPEPSolverBase::SetProblemType(SlepcEigenvalueSolver::ProblemType type)
1091 : {
1092 0 : switch (type)
1093 : {
1094 0 : case ProblemType::HERMITIAN:
1095 : case ProblemType::GEN_HERMITIAN:
1096 0 : PalacePetscCall(PEPSetProblemType(pep, PEP_HERMITIAN));
1097 : break;
1098 0 : case ProblemType::NON_HERMITIAN:
1099 : case ProblemType::GEN_INDEFINITE:
1100 : case ProblemType::GEN_NON_HERMITIAN:
1101 : case ProblemType::GENERAL:
1102 0 : PalacePetscCall(PEPSetProblemType(pep, PEP_GENERAL));
1103 : break;
1104 0 : case ProblemType::HYPERBOLIC:
1105 0 : PalacePetscCall(PEPSetProblemType(pep, PEP_HYPERBOLIC));
1106 : break;
1107 0 : case ProblemType::GYROSCOPIC:
1108 0 : PalacePetscCall(PEPSetProblemType(pep, PEP_GYROSCOPIC));
1109 : break;
1110 : }
1111 0 : }
1112 :
1113 0 : void SlepcPEPSolverBase::SetType(SlepcEigenvalueSolver::Type type)
1114 : {
1115 0 : switch (type)
1116 : {
1117 0 : case Type::TOAR:
1118 0 : PalacePetscCall(PEPSetType(pep, PEPTOAR));
1119 : break;
1120 0 : case Type::STOAR:
1121 0 : PalacePetscCall(PEPSetType(pep, PEPSTOAR));
1122 : break;
1123 0 : case Type::QARNOLDI:
1124 0 : PalacePetscCall(PEPSetType(pep, PEPQARNOLDI));
1125 : break;
1126 0 : case Type::JD:
1127 0 : PalacePetscCall(PEPSetType(pep, PEPJD));
1128 0 : region = false;
1129 0 : break;
1130 0 : case Type::KRYLOVSCHUR:
1131 : case Type::POWER:
1132 : case Type::SUBSPACE:
1133 : case Type::SLP:
1134 : case Type::NLEIGS:
1135 0 : MFEM_ABORT("Eigenvalue solver type not implemented!");
1136 : break;
1137 : }
1138 0 : }
1139 :
1140 0 : void SlepcPEPSolverBase::SetInitialSpace(const ComplexVector &v)
1141 : {
1142 0 : MFEM_VERIFY(
1143 : A0 && A1 && A2,
1144 : "Must call SetOperators before using SetInitialSpace for SLEPc eigenvalue solver!");
1145 0 : if (!v0)
1146 : {
1147 0 : PalacePetscCall(MatCreateVecs(A0, nullptr, &v0));
1148 : }
1149 0 : PalacePetscCall(ToPetscVec(v, v0));
1150 0 : Vec is[1] = {v0};
1151 0 : PalacePetscCall(PEPSetInitialSpace(pep, 1, is));
1152 0 : }
1153 :
1154 0 : void SlepcPEPSolverBase::Customize()
1155 : {
1156 0 : SlepcEigenvalueSolver::Customize();
1157 0 : PalacePetscCall(PEPSetTarget(pep, sigma / gamma));
1158 0 : if (!cl_custom)
1159 : {
1160 0 : PalacePetscCall(PEPSetFromOptions(pep));
1161 0 : if (print > 0)
1162 : {
1163 0 : PetscOptionsView(nullptr, PETSC_VIEWER_STDOUT_(GetComm()));
1164 0 : Mpi::Print(GetComm(), "\n");
1165 : }
1166 0 : cl_custom = true;
1167 : }
1168 0 : }
1169 :
1170 0 : int SlepcPEPSolverBase::Solve()
1171 : {
1172 0 : MFEM_VERIFY(A0 && A1 && A2 && opInv, "Operators are not set for SlepcPEPSolverBase!");
1173 :
1174 : // Solve the eigenvalue problem.
1175 : PetscInt num_conv;
1176 0 : Customize();
1177 0 : PalacePetscCall(PEPSolve(pep));
1178 0 : PalacePetscCall(PEPGetConverged(pep, &num_conv));
1179 0 : if (print > 0)
1180 : {
1181 0 : Mpi::Print(GetComm(), "\n");
1182 0 : PalacePetscCall(PEPConvergedReasonView(pep, PETSC_VIEWER_STDOUT_(GetComm())));
1183 0 : Mpi::Print(GetComm(),
1184 : " Total number of linear systems solved: {:d}\n"
1185 : " Total number of linear solver iterations: {:d}\n",
1186 0 : opInv->NumTotalMult(), opInv->NumTotalMultIterations());
1187 : }
1188 :
1189 : // Compute and store the eigenpair residuals.
1190 0 : RescaleEigenvectors(num_conv);
1191 0 : return (int)num_conv;
1192 : }
1193 :
1194 0 : std::complex<double> SlepcPEPSolverBase::GetEigenvalue(int i) const
1195 : {
1196 0 : PetscScalar l;
1197 0 : PalacePetscCall(PEPGetEigenpair(pep, i, &l, nullptr, nullptr, nullptr));
1198 0 : return l * gamma;
1199 : }
1200 :
1201 0 : void SlepcPEPSolverBase::GetEigenvector(int i, ComplexVector &x) const
1202 : {
1203 0 : MFEM_VERIFY(
1204 : v0,
1205 : "Must call SetOperators before using GetEigenvector for SLEPc eigenvalue solver!");
1206 0 : PalacePetscCall(PEPGetEigenpair(pep, i, nullptr, nullptr, v0, nullptr));
1207 0 : PalacePetscCall(FromPetscVec(v0, x));
1208 0 : if (xscale.get()[i] > 0.0)
1209 : {
1210 0 : x *= xscale.get()[i];
1211 : }
1212 0 : }
1213 :
1214 0 : BV SlepcPEPSolverBase::GetBV() const
1215 : {
1216 : BV bv;
1217 0 : PalacePetscCall(PEPGetBV(pep, &bv));
1218 0 : return bv;
1219 : }
1220 :
1221 0 : ST SlepcPEPSolverBase::GetST() const
1222 : {
1223 : ST st;
1224 0 : PalacePetscCall(PEPGetST(pep, &st));
1225 0 : return st;
1226 : }
1227 :
1228 0 : RG SlepcPEPSolverBase::GetRG() const
1229 : {
1230 : RG rg;
1231 0 : PalacePetscCall(PEPGetRG(pep, &rg));
1232 0 : return rg;
1233 : }
1234 :
1235 0 : SlepcPEPSolver::SlepcPEPSolver(MPI_Comm comm, int print, const std::string &prefix)
1236 0 : : SlepcPEPSolverBase(comm, print, prefix)
1237 : {
1238 0 : opK = opC = opM = nullptr;
1239 0 : normK = normC = normM = 0.0;
1240 0 : }
1241 :
1242 0 : void SlepcPEPSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &C,
1243 : const ComplexOperator &M,
1244 : EigenvalueSolver::ScaleType type)
1245 : {
1246 : // Construct shell matrices for the scaled operators which define the quadratic polynomial
1247 : // eigenvalue problem.
1248 0 : const bool first = (opK == nullptr);
1249 0 : opK = &K;
1250 0 : opC = &C;
1251 0 : opM = &M;
1252 :
1253 0 : if (first)
1254 : {
1255 0 : const PetscInt n = opK->Height();
1256 0 : PalacePetscCall(
1257 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &A0));
1258 0 : PalacePetscCall(
1259 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &A1));
1260 0 : PalacePetscCall(
1261 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &A2));
1262 0 : PalacePetscCall(
1263 : MatShellSetOperation(A0, MATOP_MULT, (void (*)(void))__mat_apply_PEP_A0));
1264 0 : PalacePetscCall(
1265 : MatShellSetOperation(A1, MATOP_MULT, (void (*)(void))__mat_apply_PEP_A1));
1266 0 : PalacePetscCall(
1267 : MatShellSetOperation(A2, MATOP_MULT, (void (*)(void))__mat_apply_PEP_A2));
1268 0 : PalacePetscCall(MatShellSetVecType(A0, PetscVecType()));
1269 0 : PalacePetscCall(MatShellSetVecType(A1, PetscVecType()));
1270 0 : PalacePetscCall(MatShellSetVecType(A2, PetscVecType()));
1271 0 : Mat A[3] = {A0, A1, A2};
1272 0 : PalacePetscCall(PEPSetOperators(pep, 3, A));
1273 : }
1274 :
1275 0 : if (first && type != ScaleType::NONE)
1276 : {
1277 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
1278 0 : normC = linalg::SpectralNorm(GetComm(), *opC, opC->IsReal());
1279 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
1280 0 : MFEM_VERIFY(normK >= 0.0 && normC >= 0.0 && normM >= 0.0,
1281 : "Invalid matrix norms for PEP scaling!");
1282 0 : if (normK > 0 && normC >= 0.0 && normM > 0.0)
1283 : {
1284 0 : gamma = std::sqrt(normK / normM);
1285 0 : delta = 2.0 / (normK + gamma * normC);
1286 : }
1287 : }
1288 :
1289 : // Set up workspace.
1290 0 : if (!v0)
1291 : {
1292 0 : PalacePetscCall(MatCreateVecs(A0, nullptr, &v0));
1293 : }
1294 0 : x1.SetSize(opK->Height());
1295 0 : y1.SetSize(opK->Height());
1296 :
1297 : // Configure linear solver.
1298 0 : if (first)
1299 : {
1300 0 : ConfigurePCShell(GetST(), (void *)this, __pc_apply_PEP);
1301 : }
1302 0 : }
1303 :
1304 0 : void SlepcPEPSolver::SetBMat(const Operator &B)
1305 : {
1306 0 : SlepcEigenvalueSolver::SetBMat(B);
1307 :
1308 0 : const PetscInt n = B.Height();
1309 0 : PalacePetscCall(
1310 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &B0));
1311 0 : PalacePetscCall(MatShellSetOperation(B0, MATOP_MULT, (void (*)(void))__mat_apply_PEP_B));
1312 0 : PalacePetscCall(MatShellSetVecType(B0, PetscVecType()));
1313 :
1314 0 : BV bv = GetBV();
1315 0 : PalacePetscCall(BVSetMatrix(bv, B0, PETSC_FALSE));
1316 0 : }
1317 :
1318 0 : PetscReal SlepcPEPSolver::GetResidualNorm(PetscScalar l, const ComplexVector &x,
1319 : ComplexVector &r) const
1320 : {
1321 : // Compute the i-th eigenpair residual: || P(λ) x ||₂ = || (K + λ C + λ² M) x ||₂ for
1322 : // eigenvalue λ.
1323 0 : opK->Mult(x, r);
1324 0 : if (opC)
1325 : {
1326 0 : opC->AddMult(x, r, l);
1327 : }
1328 0 : opM->AddMult(x, r, l * l);
1329 0 : return linalg::Norml2(GetComm(), r);
1330 : }
1331 :
1332 0 : PetscReal SlepcPEPSolver::GetBackwardScaling(PetscScalar l) const
1333 : {
1334 : // Make sure not to use norms from scaling as this can be confusing if they are different.
1335 : // Note that SLEPc typically uses ||.||∞, not Frobenius.
1336 0 : if (normK <= 0.0)
1337 : {
1338 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
1339 : }
1340 0 : if (normC <= 0.0 && opC)
1341 : {
1342 0 : normC = linalg::SpectralNorm(GetComm(), *opC, opC->IsReal());
1343 : }
1344 0 : if (normM <= 0.0)
1345 : {
1346 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
1347 : }
1348 0 : PetscReal t = PetscAbsScalar(l);
1349 0 : return normK + t * normC + t * t * normM;
1350 : }
1351 :
1352 : // NEP specific methods.
1353 :
1354 0 : SlepcNEPSolverBase::SlepcNEPSolverBase(MPI_Comm comm, int print, const std::string &prefix)
1355 0 : : SlepcEigenvalueSolver(print)
1356 : {
1357 0 : PalacePetscCall(NEPCreate(comm, &nep));
1358 0 : PalacePetscCall(NEPSetOptionsPrefix(nep, prefix.c_str()));
1359 0 : if (print > 0)
1360 : {
1361 0 : std::string opts = "-nep_monitor";
1362 0 : if (print > 2)
1363 : {
1364 0 : opts.append(" -nep_view");
1365 : }
1366 0 : if (prefix.length() > 0)
1367 : {
1368 0 : PetscOptionsPrefixPush(nullptr, prefix.c_str());
1369 : }
1370 0 : PetscOptionsInsertString(nullptr, opts.c_str());
1371 0 : if (prefix.length() > 0)
1372 : {
1373 0 : PetscOptionsPrefixPop(nullptr);
1374 : }
1375 : }
1376 0 : A = J = nullptr;
1377 0 : }
1378 :
1379 0 : SlepcNEPSolverBase::~SlepcNEPSolverBase()
1380 : {
1381 0 : PalacePetscCall(NEPDestroy(&nep));
1382 0 : PalacePetscCall(MatDestroy(&A));
1383 0 : PalacePetscCall(MatDestroy(&J));
1384 0 : }
1385 :
1386 0 : void SlepcNEPSolverBase::SetNumModes(int num_eig, int num_vec)
1387 : {
1388 0 : PalacePetscCall(NEPSetDimensions(nep, num_eig, (num_vec > 0) ? num_vec : PETSC_DEFAULT,
1389 : PETSC_DEFAULT));
1390 0 : }
1391 :
1392 0 : void SlepcNEPSolverBase::SetTol(PetscReal tol)
1393 : {
1394 0 : PalacePetscCall(NEPSetTolerances(nep, tol, PETSC_DEFAULT));
1395 0 : PalacePetscCall(NEPSetConvergenceTest(nep, NEP_CONV_REL));
1396 0 : }
1397 :
1398 0 : void SlepcNEPSolverBase::SetMaxIter(int max_it)
1399 : {
1400 0 : PalacePetscCall(
1401 : NEPSetTolerances(nep, PETSC_DEFAULT, (max_it > 0) ? max_it : PETSC_DEFAULT));
1402 0 : }
1403 :
1404 0 : void SlepcNEPSolverBase::SetShiftInvert(std::complex<double> s, bool precond)
1405 : {
1406 0 : sigma = s; // Wait until solve time to call NEPSetTarget
1407 0 : sinvert = false;
1408 0 : }
1409 :
1410 0 : void SlepcNEPSolverBase::SetWhichEigenpairs(EigenvalueSolver::WhichType type)
1411 : {
1412 0 : switch (type)
1413 : {
1414 0 : case WhichType::LARGEST_MAGNITUDE:
1415 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_LARGEST_MAGNITUDE));
1416 0 : region = false;
1417 0 : break;
1418 0 : case WhichType::SMALLEST_MAGNITUDE:
1419 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_SMALLEST_MAGNITUDE));
1420 0 : region = false;
1421 0 : break;
1422 0 : case WhichType::LARGEST_REAL:
1423 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_LARGEST_REAL));
1424 : break;
1425 0 : case WhichType::SMALLEST_REAL:
1426 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_SMALLEST_REAL));
1427 : break;
1428 0 : case WhichType::LARGEST_IMAGINARY:
1429 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_LARGEST_IMAGINARY));
1430 : break;
1431 0 : case WhichType::SMALLEST_IMAGINARY:
1432 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_SMALLEST_IMAGINARY));
1433 : break;
1434 0 : case WhichType::TARGET_MAGNITUDE:
1435 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_TARGET_MAGNITUDE));
1436 0 : region = false;
1437 0 : break;
1438 0 : case WhichType::TARGET_REAL:
1439 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_TARGET_REAL));
1440 : break;
1441 0 : case WhichType::TARGET_IMAGINARY:
1442 0 : PalacePetscCall(NEPSetWhichEigenpairs(nep, NEP_TARGET_IMAGINARY));
1443 : break;
1444 : }
1445 0 : }
1446 :
1447 0 : void SlepcNEPSolverBase::SetProblemType(SlepcEigenvalueSolver::ProblemType type)
1448 : {
1449 0 : switch (type)
1450 : {
1451 0 : case ProblemType::GENERAL:
1452 0 : PalacePetscCall(NEPSetProblemType(nep, NEP_GENERAL));
1453 : break;
1454 0 : case ProblemType::HERMITIAN:
1455 : case ProblemType::GEN_HERMITIAN:
1456 : case ProblemType::NON_HERMITIAN:
1457 : case ProblemType::GEN_INDEFINITE:
1458 : case ProblemType::GEN_NON_HERMITIAN:
1459 : case ProblemType::HYPERBOLIC:
1460 : case ProblemType::GYROSCOPIC:
1461 0 : MFEM_ABORT("Problem type not implemented!");
1462 : break;
1463 : }
1464 0 : }
1465 :
1466 0 : void SlepcNEPSolverBase::SetType(SlepcEigenvalueSolver::Type type)
1467 : {
1468 0 : switch (type)
1469 : {
1470 0 : case Type::SLP:
1471 0 : PalacePetscCall(NEPSetType(nep, NEPSLP));
1472 : break;
1473 0 : case Type::NLEIGS:
1474 : case Type::KRYLOVSCHUR:
1475 : case Type::POWER:
1476 : case Type::SUBSPACE:
1477 : case Type::JD:
1478 : case Type::TOAR:
1479 : case Type::STOAR:
1480 : case Type::QARNOLDI:
1481 0 : MFEM_ABORT("Eigenvalue solver type not implemented!");
1482 : break;
1483 : }
1484 0 : }
1485 :
1486 0 : void SlepcNEPSolverBase::SetInitialSpace(const ComplexVector &v)
1487 : {
1488 0 : MFEM_VERIFY(
1489 : A && J,
1490 : "Must call SetOperators before using SetInitialSpace for SLEPc eigenvalue solver!");
1491 0 : if (!v0)
1492 : {
1493 0 : PalacePetscCall(MatCreateVecs(A, nullptr, &v0));
1494 : }
1495 0 : PalacePetscCall(ToPetscVec(v, v0));
1496 0 : Vec is[1] = {v0};
1497 0 : PalacePetscCall(NEPSetInitialSpace(nep, 1, is));
1498 0 : }
1499 :
1500 0 : void SlepcNEPSolverBase::Customize()
1501 : {
1502 : // Configure the region based on the given target if necessary.
1503 0 : PalacePetscCall(NEPSetTarget(nep, sigma));
1504 0 : if (!cl_custom)
1505 : {
1506 0 : PalacePetscCall(NEPSetFromOptions(nep));
1507 0 : if (print > 0)
1508 : {
1509 0 : PetscOptionsView(nullptr, PETSC_VIEWER_STDOUT_(GetComm()));
1510 0 : Mpi::Print(GetComm(), "\n");
1511 : }
1512 0 : cl_custom = true;
1513 : }
1514 0 : }
1515 :
1516 0 : int SlepcNEPSolverBase::Solve()
1517 : {
1518 0 : MFEM_VERIFY(A && J && opInv, "Operators are not set for SlepcNEPSolverBase!");
1519 :
1520 : // Solve the eigenvalue problem.
1521 : perm.reset();
1522 : PetscInt num_conv;
1523 0 : Customize();
1524 0 : PalacePetscCall(NEPSolve(nep));
1525 0 : PalacePetscCall(NEPGetConverged(nep, &num_conv));
1526 0 : if (print > 0)
1527 : {
1528 0 : Mpi::Print(GetComm(), "\n");
1529 0 : PalacePetscCall(NEPConvergedReasonView(nep, PETSC_VIEWER_STDOUT_(GetComm())));
1530 0 : Mpi::Print(GetComm(),
1531 : " Total number of linear systems solved: {:d}\n"
1532 : " Total number of linear solver iterations: {:d}\n",
1533 0 : opInv->NumTotalMult(), opInv->NumTotalMultIterations());
1534 : }
1535 :
1536 : // Compute and store the ordered eigenpair residuals.
1537 0 : const int nev = (int)num_conv;
1538 0 : perm = std::make_unique<int[]>(nev);
1539 0 : std::vector<std::complex<double>> eig(nev);
1540 0 : for (int i = 0; i < nev; i++)
1541 : {
1542 0 : PetscScalar l;
1543 0 : PalacePetscCall(NEPGetEigenpair(nep, i, &l, nullptr, nullptr, nullptr));
1544 0 : eig[i] = l;
1545 0 : perm[i] = i;
1546 : }
1547 : // Sort by ascending imaginary component.
1548 0 : std::sort(perm.get(), perm.get() + nev,
1549 0 : [&eig](auto l, auto r) { return eig[l].imag() < eig[r].imag(); });
1550 0 : RescaleEigenvectors(nev);
1551 0 : return nev;
1552 : }
1553 :
1554 0 : std::complex<double> SlepcNEPSolverBase::GetEigenvalue(int i) const
1555 : {
1556 0 : PetscScalar l;
1557 0 : const int &j = perm.get()[i];
1558 0 : PalacePetscCall(NEPGetEigenpair(nep, j, &l, nullptr, nullptr, nullptr));
1559 0 : return l;
1560 : }
1561 :
1562 0 : void SlepcNEPSolverBase::GetEigenvector(int i, ComplexVector &x) const
1563 : {
1564 0 : MFEM_VERIFY(
1565 : v0,
1566 : "Must call SetOperators before using GetEigenvector for SLEPc eigenvalue solver!");
1567 0 : const int &j = perm.get()[i];
1568 0 : PalacePetscCall(NEPGetEigenpair(nep, j, nullptr, nullptr, v0, nullptr));
1569 0 : PalacePetscCall(FromPetscVec(v0, x));
1570 0 : if (xscale.get()[i] > 0.0)
1571 : {
1572 0 : x *= xscale.get()[i];
1573 : }
1574 0 : }
1575 :
1576 0 : BV SlepcNEPSolverBase::GetBV() const
1577 : {
1578 : BV bv;
1579 0 : PalacePetscCall(NEPGetBV(nep, &bv));
1580 0 : return bv;
1581 : }
1582 :
1583 0 : ST SlepcNEPSolverBase::GetST() const
1584 : {
1585 : ST st = nullptr;
1586 : // NEPGetST does not exist.
1587 0 : return st;
1588 : }
1589 :
1590 0 : RG SlepcNEPSolverBase::GetRG() const
1591 : {
1592 : RG rg;
1593 0 : PalacePetscCall(NEPGetRG(nep, &rg));
1594 0 : return rg;
1595 : }
1596 :
1597 0 : SlepcNEPSolver::SlepcNEPSolver(MPI_Comm comm, int print, const std::string &prefix)
1598 0 : : SlepcNEPSolverBase(comm, print, prefix)
1599 : {
1600 0 : opK = opC = opM = nullptr;
1601 0 : normK = normC = normM = 0.0;
1602 0 : }
1603 :
1604 0 : void SlepcNEPSolver::SetExtraSystemMatrix(
1605 : std::function<std::unique_ptr<ComplexOperator>(double)> A2)
1606 : {
1607 0 : funcA2 = A2;
1608 0 : }
1609 :
1610 0 : void SlepcNEPSolver::SetPreconditionerUpdate(
1611 : std::function<std::unique_ptr<ComplexOperator>(
1612 : std::complex<double>, std::complex<double>, std::complex<double>, double)>
1613 : P)
1614 : {
1615 0 : funcP = P;
1616 0 : }
1617 :
1618 0 : void SlepcNEPSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &M,
1619 : EigenvalueSolver::ScaleType type)
1620 : {
1621 : // Construct shell matrices for the scaled operators which define the quadratic polynomial
1622 : // eigenvalue problem.
1623 0 : const bool first = (opK == nullptr);
1624 0 : opK = &K;
1625 0 : opM = &M;
1626 :
1627 0 : if (first)
1628 : {
1629 0 : const PetscInt n = opK->Height();
1630 0 : PalacePetscCall(
1631 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &A));
1632 0 : PalacePetscCall(
1633 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &J));
1634 0 : PalacePetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))__mat_apply_NEP_A));
1635 0 : PalacePetscCall(MatShellSetOperation(J, MATOP_MULT, (void (*)(void))__mat_apply_NEP_J));
1636 0 : PalacePetscCall(MatShellSetVecType(A, PetscVecType()));
1637 0 : PalacePetscCall(MatShellSetVecType(J, PetscVecType()));
1638 0 : PalacePetscCall(NEPSetFunction(nep, A, A, __form_NEP_function, NULL));
1639 0 : PalacePetscCall(NEPSetJacobian(nep, J, __form_NEP_jacobian, NULL));
1640 : }
1641 :
1642 0 : if (first && type != ScaleType::NONE)
1643 : {
1644 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
1645 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
1646 0 : MFEM_VERIFY(normK >= 0.0 && normM >= 0.0, "Invalid matrix norms for NEP scaling!");
1647 0 : if (normK > 0 && normM > 0.0)
1648 : {
1649 0 : gamma = std::sqrt(normK / normM);
1650 0 : delta = 2.0 / normK;
1651 : }
1652 : }
1653 :
1654 : // Set up workspace.
1655 0 : if (!v0)
1656 : {
1657 0 : PalacePetscCall(MatCreateVecs(A, nullptr, &v0));
1658 : }
1659 0 : x1.SetSize(opK->Height());
1660 0 : y1.SetSize(opK->Height());
1661 :
1662 : // Configure linear solver.
1663 0 : if (first)
1664 : {
1665 : // SLP.
1666 : PC pc;
1667 : KSP ksp;
1668 : EPS eps;
1669 0 : PalacePetscCall(NEPSLPGetKSP(nep, &ksp));
1670 0 : PalacePetscCall(KSPSetType(ksp, KSPPREONLY));
1671 0 : PalacePetscCall(NEPSLPGetEPS(nep, &eps));
1672 0 : PalacePetscCall(EPSSetType(eps, EPSKRYLOVSCHUR));
1673 0 : PalacePetscCall(KSPGetPC(ksp, &pc));
1674 0 : PalacePetscCall(PCSetType(pc, PCSHELL));
1675 0 : PalacePetscCall(PCShellSetContext(pc, (void *)this));
1676 0 : PalacePetscCall(PCShellSetApply(pc, __pc_apply_NEP));
1677 : }
1678 0 : }
1679 :
1680 0 : void SlepcNEPSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &C,
1681 : const ComplexOperator &M,
1682 : EigenvalueSolver::ScaleType type)
1683 : {
1684 : // Construct shell matrices for the scaled operators which define the quadratic polynomial
1685 : // eigenvalue problem.
1686 0 : const bool first = (opK == nullptr);
1687 0 : opK = &K;
1688 0 : opC = &C;
1689 0 : opM = &M;
1690 :
1691 0 : if (first)
1692 : {
1693 0 : const PetscInt n = opK->Height();
1694 0 : PalacePetscCall(
1695 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &A));
1696 0 : PalacePetscCall(
1697 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &J));
1698 0 : PalacePetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))__mat_apply_NEP_A));
1699 0 : PalacePetscCall(MatShellSetOperation(J, MATOP_MULT, (void (*)(void))__mat_apply_NEP_J));
1700 0 : PalacePetscCall(MatShellSetVecType(A, PetscVecType()));
1701 0 : PalacePetscCall(MatShellSetVecType(J, PetscVecType()));
1702 0 : PalacePetscCall(NEPSetFunction(nep, A, A, __form_NEP_function, NULL));
1703 0 : PalacePetscCall(NEPSetJacobian(nep, J, __form_NEP_jacobian, NULL));
1704 : }
1705 :
1706 0 : if (first && type != ScaleType::NONE)
1707 : {
1708 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
1709 0 : normC = linalg::SpectralNorm(GetComm(), *opC, opC->IsReal());
1710 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
1711 0 : MFEM_VERIFY(normK >= 0.0 && normC >= 0.0 && normM >= 0.0,
1712 : "Invalid matrix norms for NEP scaling!");
1713 0 : if (normK > 0 && normC >= 0.0 && normM > 0.0)
1714 : {
1715 0 : gamma = std::sqrt(normK / normM);
1716 0 : delta = 2.0 / (normK + gamma * normC);
1717 : }
1718 : }
1719 :
1720 : // Set up workspace.
1721 0 : if (!v0)
1722 : {
1723 0 : PalacePetscCall(MatCreateVecs(A, nullptr, &v0));
1724 : }
1725 0 : x1.SetSize(opK->Height());
1726 0 : y1.SetSize(opK->Height());
1727 :
1728 : // Configure linear solver.
1729 0 : if (first)
1730 : {
1731 : // SLP.
1732 : PC pc;
1733 : KSP ksp;
1734 : EPS eps;
1735 0 : PalacePetscCall(NEPSLPGetKSP(nep, &ksp));
1736 0 : PalacePetscCall(KSPSetType(ksp, KSPPREONLY));
1737 0 : PalacePetscCall(NEPSLPGetEPS(nep, &eps));
1738 0 : PalacePetscCall(EPSSetType(eps, EPSKRYLOVSCHUR));
1739 0 : PalacePetscCall(KSPGetPC(ksp, &pc));
1740 0 : PalacePetscCall(PCSetType(pc, PCSHELL));
1741 0 : PalacePetscCall(PCShellSetContext(pc, (void *)this));
1742 0 : PalacePetscCall(PCShellSetApply(pc, __pc_apply_NEP));
1743 : }
1744 0 : }
1745 :
1746 0 : void SlepcNEPSolver::SetBMat(const Operator &B)
1747 : {
1748 0 : SlepcEigenvalueSolver::SetBMat(B);
1749 :
1750 0 : const PetscInt n = B.Height();
1751 0 : PalacePetscCall(
1752 : MatCreateShell(GetComm(), n, n, PETSC_DECIDE, PETSC_DECIDE, (void *)this, &B0));
1753 0 : PalacePetscCall(MatShellSetOperation(B0, MATOP_MULT, (void (*)(void))__mat_apply_NEP_B));
1754 0 : PalacePetscCall(MatShellSetVecType(B0, PetscVecType()));
1755 :
1756 0 : BV bv = GetBV();
1757 0 : PalacePetscCall(BVSetMatrix(bv, B0, PETSC_FALSE));
1758 0 : }
1759 :
1760 0 : PetscReal SlepcNEPSolver::GetResidualNorm(PetscScalar l, const ComplexVector &x,
1761 : ComplexVector &r) const
1762 : {
1763 : // Compute the i-th eigenpair residual: || P(λ) x ||₂ = || (K + λ C + λ² M) x ||₂ for
1764 : // eigenvalue λ.
1765 0 : opK->Mult(x, r);
1766 0 : if (opC)
1767 : {
1768 0 : opC->AddMult(x, r, l);
1769 : }
1770 0 : opM->AddMult(x, r, l * l);
1771 0 : if (funcA2)
1772 : {
1773 0 : auto A2 = (*funcA2)(std::abs(l.imag()));
1774 0 : A2->AddMult(x, r, 1.0 + 0.0i);
1775 : }
1776 0 : return linalg::Norml2(GetComm(), r);
1777 : }
1778 :
1779 0 : PetscReal SlepcNEPSolver::GetBackwardScaling(PetscScalar l) const
1780 : {
1781 : // Make sure not to use norms from scaling as this can be confusing if they are different.
1782 : // Note that SLEPc typically uses ||.||∞, not Frobenius.
1783 0 : if (normK <= 0.0)
1784 : {
1785 0 : normK = linalg::SpectralNorm(GetComm(), *opK, opK->IsReal());
1786 : }
1787 0 : if (normC <= 0.0 && opC)
1788 : {
1789 0 : normC = linalg::SpectralNorm(GetComm(), *opC, opC->IsReal());
1790 : }
1791 0 : if (normM <= 0.0)
1792 : {
1793 0 : normM = linalg::SpectralNorm(GetComm(), *opM, opM->IsReal());
1794 : }
1795 0 : PetscReal t = PetscAbsScalar(l);
1796 0 : return normK + t * normC + t * t * normM;
1797 : }
1798 :
1799 : } // namespace palace::slepc
1800 :
1801 0 : PetscErrorCode __mat_apply_EPS_A0(Mat A, Vec x, Vec y)
1802 : {
1803 : PetscFunctionBeginUser;
1804 : palace::slepc::SlepcEPSSolver *ctx;
1805 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
1806 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
1807 :
1808 0 : PetscCall(FromPetscVec(x, ctx->x1));
1809 0 : ctx->opK->Mult(ctx->x1, ctx->y1);
1810 0 : ctx->y1 *= ctx->delta;
1811 0 : PetscCall(ToPetscVec(ctx->y1, y));
1812 :
1813 : PetscFunctionReturn(PETSC_SUCCESS);
1814 : }
1815 :
1816 0 : PetscErrorCode __mat_apply_EPS_A1(Mat A, Vec x, Vec y)
1817 : {
1818 : PetscFunctionBeginUser;
1819 : palace::slepc::SlepcEPSSolver *ctx;
1820 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
1821 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
1822 :
1823 0 : PetscCall(FromPetscVec(x, ctx->x1));
1824 0 : ctx->opM->Mult(ctx->x1, ctx->y1);
1825 0 : ctx->y1 *= ctx->delta * ctx->gamma;
1826 0 : PetscCall(ToPetscVec(ctx->y1, y));
1827 :
1828 : PetscFunctionReturn(PETSC_SUCCESS);
1829 : }
1830 :
1831 0 : PetscErrorCode __mat_apply_EPS_B(Mat A, Vec x, Vec y)
1832 : {
1833 : PetscFunctionBeginUser;
1834 : palace::slepc::SlepcEPSSolver *ctx;
1835 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
1836 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
1837 :
1838 0 : PetscCall(FromPetscVec(x, ctx->x1));
1839 0 : ctx->opB->Mult(ctx->x1.Real(), ctx->y1.Real());
1840 0 : ctx->opB->Mult(ctx->x1.Imag(), ctx->y1.Imag());
1841 0 : ctx->y1 *= ctx->delta * ctx->gamma;
1842 0 : PetscCall(ToPetscVec(ctx->y1, y));
1843 :
1844 : PetscFunctionReturn(PETSC_SUCCESS);
1845 : }
1846 :
1847 0 : PetscErrorCode __pc_apply_EPS(PC pc, Vec x, Vec y)
1848 : {
1849 : // Solve the linear system associated with the generalized eigenvalue problem: y =
1850 : // M⁻¹ x, or shift-and-invert spectral transformation: y = (K - σ M)⁻¹ x . Enforces the
1851 : // divergence-free constraint using the supplied projector.
1852 : PetscFunctionBeginUser;
1853 : palace::slepc::SlepcEPSSolver *ctx;
1854 0 : PetscCall(PCShellGetContext(pc, (void **)&ctx));
1855 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell PC context for SLEPc!");
1856 :
1857 0 : PetscCall(FromPetscVec(x, ctx->x1));
1858 0 : ctx->opInv->Mult(ctx->x1, ctx->y1);
1859 0 : if (!ctx->sinvert)
1860 : {
1861 0 : ctx->y1 *= 1.0 / (ctx->delta * ctx->gamma);
1862 : }
1863 : else
1864 : {
1865 0 : ctx->y1 *= 1.0 / ctx->delta;
1866 : }
1867 0 : if (ctx->opProj)
1868 : {
1869 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
1870 0 : ctx->opProj->Mult(ctx->y1);
1871 : // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
1872 : }
1873 0 : PetscCall(ToPetscVec(ctx->y1, y));
1874 :
1875 : PetscFunctionReturn(PETSC_SUCCESS);
1876 : }
1877 :
1878 0 : PetscErrorCode __mat_apply_PEPLinear_L0(Mat A, Vec x, Vec y)
1879 : {
1880 : // Apply the linearized operator L₀ = [ 0 I ]
1881 : // [ -K -C ] .
1882 : PetscFunctionBeginUser;
1883 : palace::slepc::SlepcPEPLinearSolver *ctx;
1884 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
1885 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
1886 0 : PetscCall(FromPetscVec(x, ctx->x1, ctx->x2));
1887 0 : ctx->y1 = ctx->x2;
1888 0 : if (ctx->opC)
1889 : {
1890 0 : ctx->opC->Mult(ctx->x2, ctx->y2);
1891 : }
1892 : else
1893 : {
1894 0 : ctx->y2 = 0.0;
1895 : }
1896 0 : ctx->y2 *= ctx->gamma;
1897 0 : ctx->opK->AddMult(ctx->x1, ctx->y2, std::complex<double>(1.0, 0.0));
1898 0 : ctx->y2 *= -ctx->delta;
1899 0 : PetscCall(ToPetscVec(ctx->y1, ctx->y2, y));
1900 :
1901 : PetscFunctionReturn(PETSC_SUCCESS);
1902 : }
1903 :
1904 0 : PetscErrorCode __mat_apply_PEPLinear_L1(Mat A, Vec x, Vec y)
1905 : {
1906 : // Apply the linearized operator L₁ = [ I 0 ]
1907 : // [ 0 M ] .
1908 : PetscFunctionBeginUser;
1909 : palace::slepc::SlepcPEPLinearSolver *ctx;
1910 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
1911 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
1912 :
1913 0 : PetscCall(FromPetscVec(x, ctx->x1, ctx->x2));
1914 0 : ctx->y1 = ctx->x1;
1915 0 : ctx->opM->Mult(ctx->x2, ctx->y2);
1916 0 : ctx->y2 *= ctx->delta * ctx->gamma * ctx->gamma;
1917 0 : PetscCall(ToPetscVec(ctx->y1, ctx->y2, y));
1918 :
1919 : PetscFunctionReturn(PETSC_SUCCESS);
1920 : }
1921 :
1922 0 : PetscErrorCode __mat_apply_PEPLinear_B(Mat A, Vec x, Vec y)
1923 : {
1924 : PetscFunctionBeginUser;
1925 : palace::slepc::SlepcPEPLinearSolver *ctx;
1926 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
1927 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
1928 :
1929 0 : PetscCall(FromPetscVec(x, ctx->x1, ctx->x2));
1930 0 : ctx->opB->Mult(ctx->x1.Real(), ctx->y1.Real());
1931 0 : ctx->opB->Mult(ctx->x1.Imag(), ctx->y1.Imag());
1932 0 : ctx->opB->Mult(ctx->x2.Real(), ctx->y2.Real());
1933 0 : ctx->opB->Mult(ctx->x2.Imag(), ctx->y2.Imag());
1934 0 : ctx->y1 *= ctx->delta * ctx->gamma * ctx->gamma;
1935 0 : ctx->y2 *= ctx->delta * ctx->gamma * ctx->gamma;
1936 0 : PetscCall(ToPetscVec(ctx->y1, ctx->y2, y));
1937 :
1938 : PetscFunctionReturn(PETSC_SUCCESS);
1939 : }
1940 :
1941 0 : PetscErrorCode __pc_apply_PEPLinear(PC pc, Vec x, Vec y)
1942 : {
1943 : // Solve the linear system associated with the generalized eigenvalue problem after
1944 : // linearization: y = L₁⁻¹ x, or with the shift-and-invert spectral transformation:
1945 : // y = (L₀ - σ L₁)⁻¹ x, with:
1946 : // L₀ = [ 0 I ] L₁ = [ I 0 ]
1947 : // [ -K -C ] , [ 0 M ] .
1948 : // Enforces the divergence-free constraint using the supplied projector.
1949 : PetscFunctionBeginUser;
1950 : palace::slepc::SlepcPEPLinearSolver *ctx;
1951 0 : PetscCall(PCShellGetContext(pc, (void **)&ctx));
1952 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell PC context for SLEPc!");
1953 :
1954 0 : PetscCall(FromPetscVec(x, ctx->x1, ctx->x2));
1955 0 : if (!ctx->sinvert)
1956 : {
1957 0 : ctx->y1 = ctx->x1;
1958 0 : if (ctx->opProj)
1959 : {
1960 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
1961 0 : ctx->opProj->Mult(ctx->y1);
1962 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
1963 : }
1964 :
1965 0 : ctx->opInv->Mult(ctx->x2, ctx->y2);
1966 0 : ctx->y2 *= 1.0 / (ctx->delta * ctx->gamma * ctx->gamma);
1967 0 : if (ctx->opProj)
1968 : {
1969 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y2));
1970 0 : ctx->opProj->Mult(ctx->y2);
1971 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y2));
1972 : }
1973 : }
1974 : else
1975 : {
1976 0 : ctx->y1.AXPBY(-ctx->sigma / (ctx->delta * ctx->gamma), ctx->x2, 0.0); // Temporarily
1977 0 : ctx->opK->AddMult(ctx->x1, ctx->y1, std::complex<double>(1.0, 0.0));
1978 0 : ctx->opInv->Mult(ctx->y1, ctx->y2);
1979 0 : if (ctx->opProj)
1980 : {
1981 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y2));
1982 0 : ctx->opProj->Mult(ctx->y2);
1983 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y2));
1984 : }
1985 :
1986 0 : ctx->y1.AXPBYPCZ(ctx->gamma / ctx->sigma, ctx->y2, -ctx->gamma / ctx->sigma, ctx->x1,
1987 : 0.0);
1988 0 : if (ctx->opProj)
1989 : {
1990 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
1991 0 : ctx->opProj->Mult(ctx->y1);
1992 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
1993 : }
1994 : }
1995 0 : PetscCall(ToPetscVec(ctx->y1, ctx->y2, y));
1996 :
1997 : PetscFunctionReturn(PETSC_SUCCESS);
1998 : }
1999 :
2000 0 : PetscErrorCode __mat_apply_PEP_A0(Mat A, Vec x, Vec y)
2001 : {
2002 : PetscFunctionBeginUser;
2003 : palace::slepc::SlepcPEPSolver *ctx;
2004 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
2005 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
2006 :
2007 0 : PetscCall(FromPetscVec(x, ctx->x1));
2008 0 : ctx->opK->Mult(ctx->x1, ctx->y1);
2009 0 : PetscCall(ToPetscVec(ctx->y1, y));
2010 :
2011 : PetscFunctionReturn(PETSC_SUCCESS);
2012 : }
2013 :
2014 0 : PetscErrorCode __mat_apply_PEP_A1(Mat A, Vec x, Vec y)
2015 : {
2016 : PetscFunctionBeginUser;
2017 : palace::slepc::SlepcPEPSolver *ctx;
2018 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
2019 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
2020 :
2021 0 : PetscCall(FromPetscVec(x, ctx->x1));
2022 0 : if (ctx->opC)
2023 : {
2024 0 : ctx->opC->Mult(ctx->x1, ctx->y1);
2025 : }
2026 : else
2027 : {
2028 0 : ctx->y1 = 0.0;
2029 : }
2030 0 : PetscCall(ToPetscVec(ctx->y1, y));
2031 :
2032 : PetscFunctionReturn(PETSC_SUCCESS);
2033 : }
2034 :
2035 0 : PetscErrorCode __mat_apply_PEP_A2(Mat A, Vec x, Vec y)
2036 : {
2037 : PetscFunctionBeginUser;
2038 : palace::slepc::SlepcPEPSolver *ctx;
2039 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
2040 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
2041 :
2042 0 : PetscCall(FromPetscVec(x, ctx->x1));
2043 0 : ctx->opM->Mult(ctx->x1, ctx->y1);
2044 0 : PetscCall(ToPetscVec(ctx->y1, y));
2045 :
2046 : PetscFunctionReturn(PETSC_SUCCESS);
2047 : }
2048 :
2049 0 : PetscErrorCode __mat_apply_PEP_B(Mat A, Vec x, Vec y)
2050 : {
2051 : PetscFunctionBeginUser;
2052 : palace::slepc::SlepcPEPSolver *ctx;
2053 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
2054 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
2055 :
2056 0 : PetscCall(FromPetscVec(x, ctx->x1));
2057 0 : ctx->opB->Mult(ctx->x1.Real(), ctx->y1.Real());
2058 0 : ctx->opB->Mult(ctx->x1.Imag(), ctx->y1.Imag());
2059 0 : ctx->y1 *= ctx->delta * ctx->gamma;
2060 0 : PetscCall(ToPetscVec(ctx->y1, y));
2061 :
2062 : PetscFunctionReturn(PETSC_SUCCESS);
2063 : }
2064 :
2065 0 : PetscErrorCode __pc_apply_PEP(PC pc, Vec x, Vec y)
2066 : {
2067 : // Solve the linear system associated with the generalized eigenvalue problem: y = M⁻¹ x,
2068 : // or shift-and-invert spectral transformation: y = P(σ)⁻¹ x . Enforces the divergence-
2069 : // free constraint using the supplied projector.
2070 : PetscFunctionBeginUser;
2071 : palace::slepc::SlepcPEPSolver *ctx;
2072 0 : PetscCall(PCShellGetContext(pc, (void **)&ctx));
2073 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell PC context for SLEPc!");
2074 :
2075 0 : PetscCall(FromPetscVec(x, ctx->x1));
2076 0 : ctx->opInv->Mult(ctx->x1, ctx->y1);
2077 0 : if (!ctx->sinvert)
2078 : {
2079 0 : ctx->y1 *= 1.0 / (ctx->delta * ctx->gamma * ctx->gamma);
2080 : }
2081 : else
2082 : {
2083 0 : ctx->y1 *= 1.0 / ctx->delta;
2084 : }
2085 0 : if (ctx->opProj)
2086 : {
2087 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
2088 0 : ctx->opProj->Mult(ctx->y1);
2089 : // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
2090 : }
2091 0 : PetscCall(ToPetscVec(ctx->y1, y));
2092 :
2093 : PetscFunctionReturn(PETSC_SUCCESS);
2094 : }
2095 :
2096 0 : PetscErrorCode __mat_apply_NEP_A(Mat A, Vec x, Vec y)
2097 : {
2098 : PetscFunctionBeginUser;
2099 : palace::slepc::SlepcNEPSolver *ctx;
2100 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
2101 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
2102 0 : PetscCall(FromPetscVec(x, ctx->x1));
2103 0 : ctx->opA->Mult(ctx->x1, ctx->y1);
2104 0 : PetscCall(ToPetscVec(ctx->y1, y));
2105 : PetscFunctionReturn(PETSC_SUCCESS);
2106 : }
2107 :
2108 0 : PetscErrorCode __mat_apply_NEP_J(Mat J, Vec x, Vec y)
2109 : {
2110 : PetscFunctionBeginUser;
2111 : palace::slepc::SlepcNEPSolver *ctx;
2112 0 : PetscCall(MatShellGetContext(J, (void **)&ctx));
2113 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
2114 0 : PetscCall(FromPetscVec(x, ctx->x1));
2115 0 : ctx->opJ->Mult(ctx->x1, ctx->y1);
2116 0 : PetscCall(ToPetscVec(ctx->y1, y));
2117 : PetscFunctionReturn(PETSC_SUCCESS);
2118 : }
2119 :
2120 0 : PetscErrorCode __mat_apply_NEP_B(Mat A, Vec x, Vec y)
2121 : {
2122 : PetscFunctionBeginUser;
2123 : palace::slepc::SlepcNEPSolver *ctx;
2124 0 : PetscCall(MatShellGetContext(A, (void **)&ctx));
2125 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell matrix context for SLEPc!");
2126 0 : PetscCall(FromPetscVec(x, ctx->x1));
2127 0 : ctx->opB->Mult(ctx->x1.Real(), ctx->y1.Real());
2128 0 : ctx->opB->Mult(ctx->x1.Imag(), ctx->y1.Imag());
2129 0 : PetscCall(ToPetscVec(ctx->y1, y));
2130 : PetscFunctionReturn(PETSC_SUCCESS);
2131 : }
2132 :
2133 0 : PetscErrorCode __pc_apply_NEP(PC pc, Vec x, Vec y)
2134 : {
2135 : PetscFunctionBeginUser;
2136 : palace::slepc::SlepcNEPSolver *ctx;
2137 0 : PetscCall(PCShellGetContext(pc, (void **)&ctx));
2138 0 : MFEM_VERIFY(ctx, "Invalid PETSc shell PC context for SLEPc!");
2139 0 : PetscCall(FromPetscVec(x, ctx->x1));
2140 : // Updating PC for new λ is needed for SLP, but should not be done for NLEIGS.
2141 0 : if (ctx->new_lambda && !ctx->first_pc)
2142 : {
2143 0 : if (ctx->lambda.imag() == 0.0)
2144 0 : ctx->lambda = ctx->sigma;
2145 0 : ctx->opA2_pc = (*ctx->funcA2)(std::abs(ctx->lambda.imag()));
2146 0 : ctx->opA_pc = palace::BuildParSumOperator(
2147 0 : {1.0 + 0.0i, ctx->lambda, ctx->lambda * ctx->lambda, 1.0 + 0.0i},
2148 0 : {ctx->opK, ctx->opC, ctx->opM, ctx->opA2_pc.get()}, true);
2149 0 : ctx->opP_pc = (*ctx->funcP)(std::complex<double>(1.0, 0.0), ctx->lambda,
2150 : ctx->lambda * ctx->lambda, ctx->lambda.imag());
2151 0 : ctx->opInv->SetOperators(*ctx->opA_pc, *ctx->opP_pc);
2152 0 : ctx->new_lambda = false;
2153 : }
2154 0 : else if (ctx->first_pc)
2155 : {
2156 0 : ctx->first_pc = false;
2157 0 : ctx->new_lambda = false;
2158 : }
2159 0 : ctx->opInv->Mult(ctx->x1, ctx->y1);
2160 0 : if (ctx->opProj)
2161 : {
2162 : // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
2163 0 : ctx->opProj->Mult(ctx->y1);
2164 : // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1));
2165 : }
2166 0 : PetscCall(ToPetscVec(ctx->y1, y));
2167 : PetscFunctionReturn(PETSC_SUCCESS);
2168 : }
2169 :
2170 0 : PetscErrorCode __form_NEP_function(NEP nep, PetscScalar lambda, Mat fun, Mat B, void *ctx)
2171 : {
2172 : PetscFunctionBeginUser;
2173 : palace::slepc::SlepcNEPSolver *ctxF;
2174 0 : PetscCall(MatShellGetContext(fun, (void **)&ctxF));
2175 : // A(λ) = K + λ C + λ² M + A2(Im{λ}).
2176 0 : ctxF->opA2 = (*ctxF->funcA2)(std::abs(lambda.imag()));
2177 0 : ctxF->opA = palace::BuildParSumOperator(
2178 0 : {1.0 + 0.0i, lambda, lambda * lambda, 1.0 + 0.0i},
2179 0 : {ctxF->opK, ctxF->opC, ctxF->opM, ctxF->opA2.get()}, true);
2180 0 : ctxF->lambda = lambda;
2181 0 : ctxF->new_lambda = true; // flag to update the preconditioner in SLP
2182 0 : PetscFunctionReturn(PETSC_SUCCESS);
2183 : }
2184 :
2185 0 : PetscErrorCode __form_NEP_jacobian(NEP nep, PetscScalar lambda, Mat fun, void *ctx)
2186 : {
2187 : PetscFunctionBeginUser;
2188 : palace::slepc::SlepcNEPSolver *ctxF;
2189 0 : PetscCall(MatShellGetContext(fun, (void **)&ctxF));
2190 : // A(λ) = K + λ C + λ² M + A2(Im{λ}).
2191 : // J(λ) = C + 2 λ M + A2'(Im{λ}).
2192 0 : ctxF->opA2 = (*ctxF->funcA2)(std::abs(lambda.imag()));
2193 : const auto eps = std::sqrt(std::numeric_limits<double>::epsilon());
2194 0 : ctxF->opA2p = (*ctxF->funcA2)(std::abs(lambda.imag()) * (1.0 + eps));
2195 0 : std::complex<double> denom = std::complex<double>(0.0, eps * std::abs(lambda.imag()));
2196 0 : ctxF->opAJ = palace::BuildParSumOperator({1.0 / denom, -1.0 / denom},
2197 0 : {ctxF->opA2p.get(), ctxF->opA2.get()}, true);
2198 0 : ctxF->opJ = palace::BuildParSumOperator(
2199 0 : {0.0 + 0.0i, 1.0 + 0.0i, 2.0 * lambda, 1.0 + 0.0i},
2200 0 : {ctxF->opK, ctxF->opC, ctxF->opM, ctxF->opAJ.get()}, true);
2201 0 : PetscFunctionReturn(PETSC_SUCCESS);
2202 : }
2203 :
2204 : #endif
|
