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 "nleps.hpp"
5 :
6 : #include <algorithm>
7 : #include <string>
8 : #include <Eigen/Dense>
9 : #include <mfem.hpp>
10 : #include "linalg/divfree.hpp"
11 : #include "linalg/rap.hpp"
12 : #include "utils/communication.hpp"
13 :
14 : namespace palace
15 : {
16 :
17 : using namespace std::complex_literals;
18 : // Base class methods.
19 :
20 0 : NonLinearEigenvalueSolver::NonLinearEigenvalueSolver(MPI_Comm comm, int print)
21 0 : : comm(comm), print(print)
22 : {
23 : // Initialization.
24 0 : nleps_it = 0;
25 :
26 0 : gamma = delta = 1.0;
27 : sigma = 0.0;
28 :
29 0 : opInv = nullptr;
30 0 : opProj = nullptr;
31 0 : opB = nullptr;
32 0 : }
33 :
34 0 : void NonLinearEigenvalueSolver::SetLinearSolver(ComplexKspSolver &ksp)
35 : {
36 0 : opInv = &ksp;
37 0 : }
38 :
39 0 : void NonLinearEigenvalueSolver::SetDivFreeProjector(
40 : const DivFreeSolver<ComplexVector> &divfree)
41 : {
42 0 : opProj = &divfree;
43 0 : }
44 :
45 0 : void NonLinearEigenvalueSolver::SetBMat(const Operator &B)
46 : {
47 0 : opB = &B;
48 0 : }
49 :
50 0 : void NonLinearEigenvalueSolver::SetNumModes(int num_eig, int num_vec)
51 : {
52 0 : if (nev > 0 && num_eig != nev)
53 : {
54 : res.reset();
55 : xscale.reset();
56 : perm.reset();
57 : }
58 0 : nev = num_eig;
59 0 : }
60 :
61 0 : void NonLinearEigenvalueSolver::SetTol(double tol)
62 : {
63 0 : rtol = tol;
64 0 : }
65 :
66 0 : void NonLinearEigenvalueSolver::SetMaxIter(int max_it)
67 : {
68 0 : nleps_it = max_it;
69 0 : }
70 :
71 0 : void NonLinearEigenvalueSolver::SetWhichEigenpairs(EigenvalueSolver::WhichType type)
72 : {
73 0 : which_type = type;
74 0 : }
75 :
76 0 : void NonLinearEigenvalueSolver::SetShiftInvert(std::complex<double> s, bool precond)
77 : {
78 0 : MFEM_VERIFY(!precond, "Nonlinear eigenvalue solver does not support preconditioned "
79 : "spectral transformation option!");
80 0 : sigma = s;
81 0 : sinvert = true;
82 0 : }
83 :
84 0 : void NonLinearEigenvalueSolver::SetInitialSpace(const ComplexVector &v)
85 : {
86 0 : MFEM_ABORT("SetInitialSpace not defined for base class NonLinearEigenvalueSolver!");
87 : }
88 :
89 0 : std::complex<double> NonLinearEigenvalueSolver::GetEigenvalue(int i) const
90 : {
91 0 : MFEM_VERIFY(i >= 0 && i < nev,
92 : "Out of range eigenpair requested (i = " << i << ", nev = " << nev << ")!");
93 0 : const int &j = perm.get()[i];
94 0 : return eigenvalues[j];
95 : }
96 :
97 0 : void NonLinearEigenvalueSolver::GetEigenvector(int i, ComplexVector &x) const
98 : {
99 0 : MFEM_VERIFY(i >= 0 && i < nev,
100 : "Out of range eigenpair requested (i = " << i << ", nev = " << nev << ")!");
101 0 : MFEM_VERIFY(x.Size() == n, "Invalid size mismatch for provided eigenvector!");
102 0 : const int &j = perm.get()[i];
103 0 : x = eigenvectors[j];
104 0 : if (xscale.get()[j] > 0.0)
105 : {
106 0 : x *= xscale.get()[j];
107 : }
108 0 : }
109 :
110 0 : double NonLinearEigenvalueSolver::GetEigenvectorNorm(const ComplexVector &x,
111 : ComplexVector &Bx) const
112 : {
113 0 : if (opB)
114 : {
115 0 : return linalg::Norml2(comm, x, *opB, Bx);
116 : }
117 : else
118 : {
119 0 : return linalg::Norml2(comm, x);
120 : }
121 : }
122 :
123 0 : double NonLinearEigenvalueSolver::GetError(int i, EigenvalueSolver::ErrorType type) const
124 : {
125 0 : MFEM_VERIFY(i >= 0 && i < nev,
126 : "Out of range eigenpair requested (i = " << i << ", nev = " << nev << ")!");
127 0 : const int &j = perm.get()[i];
128 0 : switch (type)
129 : {
130 0 : case ErrorType::ABSOLUTE:
131 0 : return res.get()[j];
132 0 : case ErrorType::RELATIVE:
133 0 : return res.get()[j] / std::abs(eigenvalues[j]);
134 0 : case ErrorType::BACKWARD:
135 0 : return res.get()[j] / GetBackwardScaling(eigenvalues[j]);
136 : }
137 : return 0.0;
138 : }
139 :
140 0 : void NonLinearEigenvalueSolver::RescaleEigenvectors(int num_eig)
141 : {
142 0 : res = std::make_unique<double[]>(num_eig);
143 0 : xscale = std::make_unique<double[]>(num_eig);
144 0 : for (int i = 0; i < num_eig; i++)
145 : {
146 0 : x1 = eigenvectors[i];
147 0 : xscale.get()[i] = 1.0 / GetEigenvectorNorm(x1, y1);
148 0 : res.get()[i] = GetResidualNorm(eigenvalues[i], x1, y1) / linalg::Norml2(comm, x1);
149 : }
150 0 : }
151 :
152 : // Quasi-Newton specific methods.
153 0 : QuasiNewtonSolver::QuasiNewtonSolver(MPI_Comm comm,
154 : std::unique_ptr<EigenvalueSolver> linear_eigensolver,
155 0 : int num_conv, int print, bool refine)
156 : : NonLinearEigenvalueSolver(comm, print),
157 0 : linear_eigensolver_(std::move(linear_eigensolver)), nev_linear(num_conv),
158 0 : refine_nonlinear(refine)
159 : {
160 0 : opK = opC = opM = nullptr;
161 0 : normK = normC = normM = 0.0;
162 0 : }
163 :
164 : // Set the update frequency of the preconditioner.
165 0 : void QuasiNewtonSolver::SetPreconditionerLag(int preconditioner_update_freq,
166 : double preconditioner_update_tol)
167 : {
168 0 : preconditioner_lag = preconditioner_update_freq;
169 0 : preconditioner_tol = preconditioner_update_tol;
170 0 : }
171 :
172 : // Set the maximum number of restarts with the same initial guess.
173 0 : void QuasiNewtonSolver::SetMaxRestart(int max_num_restart)
174 : {
175 0 : max_restart = max_num_restart;
176 0 : }
177 :
178 0 : void QuasiNewtonSolver::SetExtraSystemMatrix(
179 : std::function<std::unique_ptr<ComplexOperator>(double)> A2)
180 : {
181 0 : funcA2 = A2;
182 0 : }
183 :
184 0 : void QuasiNewtonSolver::SetPreconditionerUpdate(
185 : std::function<std::unique_ptr<ComplexOperator>(
186 : std::complex<double>, std::complex<double>, std::complex<double>, double)>
187 : P)
188 : {
189 0 : funcP = P;
190 0 : }
191 :
192 0 : void QuasiNewtonSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &M,
193 : EigenvalueSolver::ScaleType type)
194 : {
195 0 : MFEM_VERIFY(!opK || K.Height() == n, "Invalid modification of eigenvalue problem size!");
196 0 : bool first = (opK == nullptr);
197 0 : opK = &K;
198 0 : opM = &M;
199 :
200 0 : if (first && type != ScaleType::NONE)
201 : {
202 0 : normK = linalg::SpectralNorm(comm, *opK, opK->IsReal());
203 0 : normM = linalg::SpectralNorm(comm, *opM, opM->IsReal());
204 0 : MFEM_VERIFY(normK >= 0.0 && normM >= 0.0,
205 : "Invalid matrix norms for Quasi-Newton scaling!");
206 0 : if (normK > 0 && normM > 0.0)
207 : {
208 0 : gamma = std::sqrt(normK / normM);
209 0 : delta = 2.0 / (normK);
210 : }
211 : }
212 :
213 0 : n = opK->Height();
214 :
215 : // Set up workspace.
216 0 : x1.SetSize(opK->Height());
217 0 : y1.SetSize(opK->Height());
218 0 : x1.UseDevice(true);
219 0 : y1.UseDevice(true);
220 0 : }
221 :
222 0 : void QuasiNewtonSolver::SetOperators(const ComplexOperator &K, const ComplexOperator &C,
223 : const ComplexOperator &M,
224 : EigenvalueSolver::ScaleType type)
225 : {
226 0 : MFEM_VERIFY(!opK || K.Height() == n, "Invalid modification of eigenvalue problem size!");
227 0 : bool first = (opK == nullptr);
228 0 : opK = &K;
229 0 : opC = &C;
230 0 : opM = &M;
231 :
232 0 : if (first && type != ScaleType::NONE)
233 : {
234 0 : normK = linalg::SpectralNorm(comm, *opK, opK->IsReal());
235 0 : normC = linalg::SpectralNorm(comm, *opC, opC->IsReal());
236 0 : normM = linalg::SpectralNorm(comm, *opM, opM->IsReal());
237 0 : MFEM_VERIFY(normK >= 0.0 && normC >= 0.0 && normM >= 0.0,
238 : "Invalid matrix norms for Quasi-Newton scaling!");
239 0 : if (normK > 0 && normC >= 0.0 && normM > 0.0)
240 : {
241 0 : gamma = std::sqrt(normK / normM);
242 0 : delta = 2.0 / (normK + gamma * normC);
243 : }
244 : }
245 :
246 0 : n = opK->Height();
247 :
248 : // Set up workspace.
249 0 : x1.SetSize(opK->Height());
250 0 : y1.SetSize(opK->Height());
251 0 : x1.UseDevice(true);
252 0 : y1.UseDevice(true);
253 0 : }
254 :
255 0 : void QuasiNewtonSolver::SetInitialGuess()
256 : {
257 :
258 0 : MFEM_VERIFY(n > 0, "Must call SetOperators before using SetInitialguess for nonlinear "
259 : "eigenvalue solver!");
260 0 : MFEM_VERIFY(nev > 0, "Must call SetNumModes before using SetInitialguess for nonlinear "
261 : "eigenvalue solver!");
262 :
263 : // Get eigenmodes initial guesses from linear eigensolver
264 0 : eigenvalues.resize(nev_linear);
265 0 : eigenvectors.resize(nev_linear);
266 0 : for (int i = 0; i < nev_linear; i++)
267 : {
268 0 : eigenvalues[i] = linear_eigensolver_->GetEigenvalue(i);
269 0 : linear_eigensolver_->GetEigenvector(i, x1);
270 : eigenvectors[i] = x1;
271 : }
272 :
273 : // Compute errors.
274 0 : RescaleEigenvectors(nev_linear);
275 :
276 : // Initialize eigenpairs ordering.
277 0 : perm = std::make_unique<int[]>(nev_linear);
278 0 : std::iota(perm.get(), perm.get() + nev_linear, 0);
279 :
280 : // Early return if nonlinear Newton won't be used.
281 0 : if (!refine_nonlinear)
282 0 : return;
283 :
284 : // If the number of initial guesses is greater than the number of requested modes
285 : // de-prioritize the initial guesses that have larger errors.
286 0 : std::vector<size_t> indices(nev_linear);
287 : std::iota(indices.begin(), indices.end(), 0);
288 0 : if (nev_linear > nev)
289 : {
290 0 : double min_error = res.get()[0];
291 0 : for (int i = 0; i < nev_linear; i++)
292 : {
293 0 : min_error = std::min(min_error, res.get()[i]);
294 : }
295 0 : const double threshold = 100.0 * min_error;
296 0 : std::sort(indices.begin(), indices.end(),
297 0 : [&](const auto i, const auto j)
298 : {
299 0 : if (res.get()[i] < threshold && res.get()[j] > threshold)
300 : {
301 : return true;
302 : }
303 0 : else if (res.get()[i] > threshold && res.get()[j] < threshold)
304 : {
305 : return false;
306 : }
307 : else
308 : {
309 0 : return eigenvalues[i].imag() < eigenvalues[j].imag();
310 : }
311 : });
312 : }
313 0 : for (int i = 0; i < nev_linear; i++)
314 : {
315 0 : eigenvalues[i] = linear_eigensolver_->GetEigenvalue(indices[i]);
316 0 : linear_eigensolver_->GetEigenvector(indices[i], x1);
317 0 : linalg::NormalizePhase(comm, x1);
318 : eigenvectors[i] = x1;
319 : }
320 :
321 : // Get ordering of the eigenpairs.
322 0 : std::sort(perm.get(), perm.get() + nev_linear, [&eig = this->eigenvalues](auto l, auto r)
323 0 : { return eig[l].imag() < eig[r].imag(); });
324 : }
325 :
326 : namespace
327 : {
328 : // Multiply an (n x k) matrix (vector of size k of ComplexVectors of size n) by a vector of
329 : // size k, returning a ComplexVector of size n.
330 0 : ComplexVector MatVecMult(const std::vector<ComplexVector> &X, const Eigen::VectorXcd &y)
331 : {
332 : MFEM_ASSERT(X.size() == y.size(), "Mismatch in dimension of input vectors!");
333 : const size_t k = X.size();
334 : const size_t n = X[0].Size();
335 : const bool use_dev = X[0].UseDevice();
336 0 : ComplexVector z;
337 0 : z.SetSize(n);
338 0 : z.UseDevice(use_dev);
339 : z = 0.0;
340 0 : for (int j = 0; j < k; j++)
341 : {
342 0 : linalg::AXPBYPCZ(y(j).real(), X[j].Real(), -y(j).imag(), X[j].Imag(), 1.0, z.Real());
343 0 : linalg::AXPBYPCZ(y(j).imag(), X[j].Real(), y(j).real(), X[j].Imag(), 1.0, z.Imag());
344 : }
345 0 : return z;
346 : }
347 :
348 : } // namespace
349 :
350 0 : int QuasiNewtonSolver::Solve()
351 : {
352 : // Quasi-Newton method for nonlinear eigenvalue problems.
353 : // Reference: Jarlebring, Koskela, Mele, Disguised and new quasi-Newton methods for
354 : // nonlinear eigenvalue problems, Numerical Algorithms (2018).
355 : // Using the deflation scheme used by SLEPc's NEP solver with minimality index set to 1.
356 : // Reference: Effenberger, Robust successive computation of eigenpairs for nonlinear
357 : // eigenvalue problems, SIAM J. Matrix Anal. Appl. (2013).
358 : // The deflation scheme solves an extended problem of size n + k, where n is the original
359 : // problem size and k is the number of converged eigenpairs. The extended operators are
360 : // never explicitly constructed and two separate vectors of length n and k are used to
361 : // store the extended solution: [v, v2] where v is a ComplexVector distributed across all
362 : // processes and v2 is an Eigen::VectorXcd stored redundantly on all processes.
363 :
364 : // Set initial guess from linear eigensolver.
365 0 : SetInitialGuess();
366 :
367 : // Return early if not refining the eigenmodes with Newton.
368 0 : if (!refine_nonlinear)
369 : {
370 0 : nev = eigenvalues.size();
371 0 : return nev;
372 : }
373 :
374 : // Palace ComplexVectors of size n.
375 0 : ComplexVector v, u, w, c, w0, z;
376 0 : v.SetSize(n);
377 0 : u.SetSize(n);
378 0 : w.SetSize(n);
379 0 : c.SetSize(n);
380 0 : w0.SetSize(n);
381 0 : z.SetSize(n);
382 0 : v.UseDevice(true);
383 0 : u.UseDevice(true);
384 0 : w.UseDevice(true);
385 0 : c.UseDevice(true);
386 0 : w0.UseDevice(true);
387 0 : z.UseDevice(true);
388 :
389 : // Eigen Matrix/Vectors for extended operator of size k.
390 : Eigen::MatrixXcd H;
391 : Eigen::VectorXcd u2, z2, c2, w2, v2;
392 :
393 : // Storage for eigenpairs.
394 : std::vector<ComplexVector> X;
395 : std::vector<std::complex<double>> eigs;
396 0 : X.reserve(nev);
397 :
398 : // Set defaults.
399 0 : if (nleps_it <= 0)
400 : {
401 0 : nleps_it = 100;
402 : }
403 :
404 : // Delta used in to compute divided difference Jacobian.
405 : const auto delta = std::sqrt(std::numeric_limits<double>::epsilon());
406 :
407 : // Set a seed and distribution for random Eigen vectors to ensure the same values on all
408 : // ranks.
409 : unsigned int seed = 111;
410 : std::mt19937 gen(seed);
411 : std::uniform_real_distribution<> dist_real(-1.0, 1.0);
412 :
413 0 : const int num_init_guess = eigenvalues.size();
414 0 : std::uniform_int_distribution<int> dist_int(0, num_init_guess - 1);
415 0 : int k = 0, restart = 0, guess_idx = 0;
416 0 : while (k < nev)
417 : {
418 : // If > max_restart with the same initial guess, skip to next initial guess.
419 : // If we tried all initial guesses and the random guess, end search even if k < nev.
420 0 : if (restart > max_restart)
421 : {
422 0 : if (guess_idx < num_init_guess)
423 : {
424 0 : guess_idx++;
425 0 : restart = 0;
426 : }
427 : else
428 : {
429 : break;
430 : }
431 : }
432 :
433 : // Set the eigenpair estimate to the initial guess.
434 0 : std::complex<double> eig, eig_opInv;
435 0 : if (guess_idx < num_init_guess)
436 : {
437 0 : eig = eigenvalues[guess_idx];
438 : v = eigenvectors[guess_idx];
439 : }
440 : else
441 : {
442 0 : eig = sigma;
443 0 : if (num_init_guess < 3)
444 : {
445 : // Set purely random vector.
446 0 : linalg::SetRandom(GetComm(), v);
447 0 : linalg::SetSubVector(
448 0 : v, *dynamic_cast<const ComplexParOperator *>(opK)->GetEssentialTrueDofs(), 0.0);
449 : }
450 : else
451 : {
452 : // Set random vector as the average of two distinct randomly-chosen initial guesses.
453 : int i1 = dist_int(gen);
454 : int i2;
455 : do
456 : {
457 : i2 = dist_int(gen);
458 0 : } while (i2 == i1);
459 0 : v.AXPBYPCZ(0.5, eigenvectors[i1], 0.5, eigenvectors[i2], 0.0);
460 : }
461 : }
462 0 : eig_opInv = eig; // eigenvalue estimate used in the (lagged) preconditioner
463 :
464 : // Set the "random" c vector and the deflation component of the eigenpair initial guess.
465 0 : linalg::SetRandom(GetComm(), c, seed); // set seed for deterministic behavior
466 0 : c2.conservativeResize(k);
467 0 : v2.conservativeResize(k);
468 0 : for (int i = 0; i < k; i++)
469 : {
470 0 : c2(i) = std::complex<double>(dist_real(gen), dist_real(gen));
471 0 : v2(i) = std::complex<double>(dist_real(gen), dist_real(gen));
472 : }
473 :
474 : // Normalize random c vector.
475 0 : double norm_c = std::sqrt(std::abs(linalg::Dot(GetComm(), c, c)) + c2.squaredNorm());
476 0 : c *= 1.0 / norm_c;
477 0 : c2 *= 1.0 / norm_c;
478 :
479 : // Normalize eigenvector estimate.
480 0 : double norm_v = std::sqrt(std::abs(linalg::Dot(GetComm(), v, v)) + v2.squaredNorm());
481 0 : v *= 1.0 / norm_v;
482 0 : v2 *= 1.0 / norm_v;
483 :
484 : // Set the linear solver operators.
485 0 : opA2 = (*funcA2)(std::abs(eig.imag()));
486 0 : opA = BuildParSumOperator({1.0 + 0.0i, eig, eig * eig, 1.0 + 0.0i},
487 0 : {opK, opC, opM, opA2.get()}, true);
488 0 : opP = (*funcP)(1.0 + 0.0i, eig, eig * eig, eig.imag());
489 0 : opInv->SetOperators(*opA, *opP);
490 :
491 : // Linear solve with the extended operator of the deflated problem.
492 0 : auto deflated_solve = [&](const ComplexVector &b1, const Eigen::VectorXcd &b2,
493 : ComplexVector &x1, Eigen::VectorXcd &x2)
494 : {
495 : // Solve the block linear system
496 : // |T(σ) U(σ)| |x1| = |b1|
497 : // |A(σ) B(σ)| |x2| |b2|
498 : // x1 = T^-1 b1
499 : // x2 = SS^-1 (b2 - A x1) where SS = (B - A T^-1 U) = - X^* X S^-1
500 : // x1 = x1 - X S x2
501 0 : opInv->Mult(b1, x1);
502 0 : if (k == 0) // no deflation
503 : {
504 0 : return;
505 : }
506 0 : x2.conservativeResize(k);
507 0 : for (int j = 0; j < k; j++)
508 : {
509 0 : x2(j) = b2(j) - linalg::Dot(GetComm(), x1, X[j]);
510 : }
511 0 : Eigen::MatrixXcd SS(k, k);
512 0 : for (int i = 0; i < k; i++)
513 : {
514 0 : for (int j = 0; j < k; j++)
515 : {
516 0 : SS(i, j) = linalg::Dot(GetComm(), X[i], X[j]);
517 : }
518 : }
519 0 : const Eigen::MatrixXcd S = eig_opInv * Eigen::MatrixXcd::Identity(k, k) - H;
520 0 : SS = -S.fullPivLu().solve(SS);
521 0 : x2 = SS.fullPivLu().solve(x2);
522 0 : const ComplexVector XSx2 = MatVecMult(X, S.fullPivLu().solve(x2));
523 0 : linalg::AXPY(-1.0, XSx2, x1);
524 0 : };
525 :
526 : // Compute w0 = T^-1 c and normalize it.
527 0 : deflated_solve(c, c2, w0, w2);
528 0 : double norm_w0 = std::sqrt(std::abs(linalg::Dot(GetComm(), w0, w0)) + w2.squaredNorm());
529 0 : w0 *= 1.0 / norm_w0;
530 0 : w2 *= 1.0 / norm_w0;
531 :
532 : // Newton iterations.
533 0 : double res = mfem::infinity();
534 0 : int it = 0, diverged_it = 0;
535 0 : while (it < nleps_it)
536 : {
537 : // Compute u = A * v.
538 0 : auto A2n = (*funcA2)(std::abs(eig.imag()));
539 0 : auto A = BuildParSumOperator({1.0 + 0.0i, eig, eig * eig, 1.0 + 0.0i},
540 0 : {opK, opC, opM, A2n.get()}, true);
541 0 : A->Mult(v, u);
542 0 : if (k > 0) // Deflation
543 : {
544 : // u1 = T(l) v1 + U(l) v2 = T(l) v1 + T(l)X(lI - H)^-1 v2.
545 0 : const Eigen::MatrixXcd S = eig * Eigen::MatrixXcd::Identity(k, k) - H;
546 0 : const ComplexVector XSv2 = MatVecMult(X, S.fullPivLu().solve(v2));
547 0 : A->AddMult(XSv2, u, 1.0);
548 : // u2 = X^* v1.
549 0 : u2.conservativeResize(k);
550 0 : for (int j = 0; j < k; j++)
551 : {
552 0 : u2(j) = linalg::Dot(GetComm(), v, X[j]);
553 : }
554 : }
555 :
556 : // Compute residual.
557 0 : res = std::sqrt(std::abs(linalg::Dot(GetComm(), u, u)) + u2.squaredNorm());
558 0 : if (print > 0)
559 : {
560 0 : Mpi::Print(GetComm(),
561 : "{:d} NLEPS (nconv={:d}, restart={:d}) residual norm {:.6e}\n", it, k,
562 : restart, res);
563 : }
564 :
565 : // End if residual below tolerance and eigenvalue above the target.
566 0 : if (res < rtol)
567 : {
568 0 : if (print > 0)
569 : {
570 0 : Mpi::Print(GetComm(),
571 : "Eigenvalue {:d}, Quasi-Newton converged in {:d} iterations "
572 : "({:.3e}{:+.3e}i).\n",
573 0 : k, it, eig.real(), eig.imag());
574 : }
575 : // Update the invariant pair with normalization.
576 0 : const auto scale = linalg::Norml2(GetComm(), v);
577 0 : v *= 1.0 / scale;
578 0 : eigs.resize(k + 1);
579 0 : eigs[k] = eig;
580 0 : X.resize(k + 1);
581 0 : X[k] = v;
582 0 : H.conservativeResizeLike(Eigen::MatrixXd::Zero(k + 1, k + 1));
583 0 : H.col(k).head(k) = v2 / scale;
584 0 : H(k, k) = eig;
585 0 : k++;
586 : // If the eigenvalue is inside the desired range, increment initial guess index
587 : // Otherwise, use the same initial guess again and increment number of desired
588 : // eigenvalues.
589 0 : if (eig.imag() > sigma.imag())
590 : {
591 0 : guess_idx++;
592 : }
593 : else
594 : {
595 0 : nev++;
596 : }
597 0 : restart = 0; // reset restart counter
598 : break;
599 : }
600 : // Stop if large residual for 10 consecutive iterations.
601 0 : diverged_it = (res > 0.9) ? diverged_it + 1 : 0;
602 0 : if (diverged_it > 10)
603 : {
604 0 : if (print > 0)
605 : {
606 0 : Mpi::Print(GetComm(),
607 : "Eigenvalue {:d}, Quasi-Newton not converging after {:d} iterations, "
608 : "restarting.\n",
609 : k, it);
610 : }
611 0 : restart++;
612 0 : break;
613 : }
614 :
615 : // Compute w = J * v.
616 0 : auto opA2p = (*funcA2)(std::abs(eig.imag()) * (1.0 + delta));
617 : const std::complex<double> denom =
618 0 : std::complex<double>(0.0, delta * std::abs(eig.imag()));
619 : std::unique_ptr<ComplexOperator> opAJ =
620 0 : BuildParSumOperator({1.0 / denom, -1.0 / denom}, {opA2p.get(), A2n.get()}, true);
621 0 : auto opJ = BuildParSumOperator({0.0 + 0.0i, 1.0 + 0.0i, 2.0 * eig, 1.0 + 0.0i},
622 0 : {opK, opC, opM, opAJ.get()}, true);
623 0 : opJ->Mult(v, w);
624 0 : if (k > 0) // Deflation
625 : {
626 : // w1 = T'(l) v1 + U'(l) v2 = T'(l) v1 + T'(l)XS v2 - T(l)XS^2 v2.
627 0 : const Eigen::MatrixXcd S = eig * Eigen::MatrixXcd::Identity(k, k) - H;
628 0 : const Eigen::VectorXcd Sv2 = S.fullPivLu().solve(v2);
629 0 : const ComplexVector XSv2 = MatVecMult(X, Sv2);
630 0 : const ComplexVector XSSv2 = MatVecMult(X, S.fullPivLu().solve(Sv2));
631 0 : opJ->AddMult(XSv2, w, 1.0);
632 0 : A->AddMult(XSSv2, w, -1.0);
633 : }
634 :
635 : // Compute delta = - dot(w0, u) / dot(w0, w).
636 0 : const std::complex<double> u2_w0 = std::complex<double>(w2.adjoint() * u2);
637 : const std::complex<double> delta =
638 0 : -(linalg::Dot(GetComm(), u, w0) + u2_w0) / linalg::Dot(GetComm(), w, w0);
639 :
640 : // Update eigenvalue.
641 : eig += delta;
642 :
643 : // Compute z = -(delta * w + u).
644 0 : z.AXPBYPCZ(-delta, w, -1.0, u, 0.0);
645 0 : z2 = -u2;
646 :
647 : // Update preconditioner if needed. Updating the preconditioner as infrequently as
648 : // possible gives the best performance and robustness. Updating the preconditioner
649 : // close to an eigenvalue can lead to numerical instability.
650 0 : if (it > 0 && it % preconditioner_lag == 0 && res > preconditioner_tol)
651 : {
652 0 : eig_opInv = eig;
653 0 : opA2 = (*funcA2)(std::abs(eig.imag()));
654 0 : opA = BuildParSumOperator({1.0 + 0.0i, eig, eig * eig, 1.0 + 0.0i},
655 0 : {opK, opC, opM, opA2.get()}, true);
656 0 : opP = (*funcP)(1.0 + 0.0i, eig, eig * eig, eig.imag());
657 0 : opInv->SetOperators(*opA, *opP);
658 : // Recompute w0 and normalize.
659 0 : deflated_solve(c, c2, w0, w2);
660 : double norm_w0 =
661 0 : std::sqrt(std::abs(linalg::Dot(GetComm(), w0, w0)) + w2.squaredNorm());
662 0 : w0 *= 1.0 / norm_w0;
663 0 : w2 *= 1.0 / norm_w0;
664 : }
665 :
666 : // Solve M (v_k+1 - v_k) = z.
667 0 : deflated_solve(z, z2, u, u2);
668 :
669 : // Update and normalize eigenvector estimate.
670 : v += u;
671 : v2 += u2;
672 0 : norm_v = std::sqrt(std::abs(linalg::Dot(GetComm(), v, v)) + v2.squaredNorm());
673 0 : v *= 1.0 / norm_v;
674 0 : v2 *= 1.0 / norm_v;
675 :
676 0 : it++;
677 0 : if (it == nleps_it)
678 : {
679 0 : if (print > 0)
680 : {
681 0 : Mpi::Print(GetComm(),
682 : "Eigenvalue {:d}, Quasi-Newton did not converge in {:d} iterations, "
683 : "restarting.\n",
684 0 : k, nleps_it);
685 : }
686 0 : restart++;
687 : }
688 : }
689 : }
690 0 : nev = k; // in case some guesses did not converge
691 :
692 : // Eigenpair extraction from the invariant pair (X, H).
693 0 : Eigen::ComplexEigenSolver<Eigen::MatrixXcd> eps;
694 0 : eps.compute(H);
695 : // H eigenvectors are ordered arbitrarily, need to match them to order of X.
696 0 : std::vector<int> order(nev), order_eigen(nev), order2(nev);
697 : std::iota(order.begin(), order.end(), 0);
698 : std::iota(order_eigen.begin(), order_eigen.end(), 0);
699 : std::iota(order2.begin(), order2.end(), 0);
700 0 : std::sort(order.begin(), order.end(),
701 0 : [&](auto l, auto r) { return eigs[l].imag() < eigs[r].imag(); });
702 0 : std::sort(order_eigen.begin(), order_eigen.end(),
703 : [&epseig = eps.eigenvalues()](auto l, auto r)
704 0 : { return epseig(l).imag() < epseig(r).imag(); });
705 0 : std::sort(order2.begin(), order2.end(),
706 0 : [&](auto l, auto r) { return order[l] < order[r]; });
707 :
708 : // Sort Eigen eigenvectors.
709 : std::vector<Eigen::VectorXcd> Xeig;
710 0 : for (int i = 0; i < nev; i++)
711 : {
712 0 : Xeig.push_back(eps.eigenvectors().col(order_eigen[i]));
713 : }
714 :
715 : // Recover the eigenvectors in the target range.
716 : eigenvalues.clear();
717 : eigenvectors.clear();
718 0 : for (int i = 0; i < nev; i++)
719 : {
720 0 : if (eigs[i].imag() > sigma.imag())
721 : {
722 0 : ComplexVector eigv = MatVecMult(X, Xeig[order2[i]]);
723 0 : eigenvalues.push_back(eigs[i]);
724 0 : eigenvectors.push_back(eigv);
725 : }
726 : }
727 0 : nev = eigenvalues.size();
728 :
729 : // Get ordering of the reported eigenpairs.
730 0 : perm = std::make_unique<int[]>(nev);
731 0 : std::iota(perm.get(), perm.get() + nev, 0);
732 0 : std::sort(perm.get(), perm.get() + nev, [&eig = this->eigenvalues](auto l, auto r)
733 0 : { return eig[l].imag() < eig[r].imag(); });
734 :
735 : // Compute the eigenpair residuals for eigenvalue λ.
736 0 : RescaleEigenvectors(nev);
737 :
738 0 : return nev;
739 0 : }
740 :
741 0 : double QuasiNewtonSolver::GetResidualNorm(std::complex<double> l, const ComplexVector &x,
742 : ComplexVector &r) const
743 : {
744 : // Compute the i-th eigenpair residual: || P(λ) x ||₂ = || (K + λ C + λ² M + A2(λ)) x ||₂
745 : // for eigenvalue λ.
746 0 : opK->Mult(x, r);
747 0 : if (opC)
748 : {
749 0 : opC->AddMult(x, r, l);
750 : }
751 0 : opM->AddMult(x, r, l * l);
752 0 : auto A2 = (*funcA2)(std::abs(l.imag()));
753 0 : A2->AddMult(x, r, 1.0);
754 0 : return linalg::Norml2(comm, r);
755 : }
756 :
757 0 : double QuasiNewtonSolver::GetBackwardScaling(std::complex<double> l) const
758 : {
759 : // Make sure not to use norms from scaling as this can be confusing if they are different.
760 : // Note that SLEPc uses ||.||∞, not the 2-norm.
761 0 : if (normK <= 0.0)
762 : {
763 0 : normK = linalg::SpectralNorm(comm, *opK, opK->IsReal());
764 : }
765 0 : if (normC <= 0.0 && opC)
766 : {
767 0 : normC = linalg::SpectralNorm(comm, *opC, opC->IsReal());
768 : }
769 0 : if (normM <= 0.0)
770 : {
771 0 : normM = linalg::SpectralNorm(comm, *opM, opM->IsReal());
772 : }
773 : double t = std::abs(l);
774 0 : return normK + t * normC + t * t * normM;
775 : }
776 :
777 0 : NewtonInterpolationOperator::NewtonInterpolationOperator(
778 0 : std::function<std::unique_ptr<ComplexOperator>(double)> funcA2, int size)
779 0 : : funcA2(funcA2)
780 : {
781 0 : rhs.SetSize(size);
782 0 : rhs.UseDevice(true);
783 0 : }
784 :
785 : // Compute the elementary symmetric polynomial. Used to convert from Newton to monomial
786 : // basis.
787 : template <typename ScalarType>
788 0 : ScalarType elementarySymmetric(const std::vector<ScalarType> &points, int k, int n)
789 : {
790 0 : if (k == 0)
791 : {
792 0 : return 1.0;
793 : }
794 0 : if (k > n || k < 0 || n == 0)
795 : {
796 0 : return 0.0;
797 : }
798 0 : return elementarySymmetric(points, k, n - 1) +
799 0 : points[n - 1] * elementarySymmetric(points, k - 1, n - 1);
800 : }
801 :
802 0 : void NewtonInterpolationOperator::Interpolate(const std::complex<double> sigma_min,
803 : const std::complex<double> sigma_max)
804 : {
805 : // Reset operators and sample points each time Interpolate is called.
806 : ops.clear();
807 0 : ops.resize(num_points);
808 : points.clear();
809 0 : points.resize(num_points);
810 :
811 : // Linearly spaced sample points.
812 0 : for (int j = 0; j < num_points; j++)
813 : {
814 0 : points[j] = sigma_min + (double)j * (sigma_max - sigma_min) / (double)(num_points - 1);
815 : }
816 :
817 : // Build divided difference matrices.
818 0 : for (int k = 0; k < num_points; k++)
819 : {
820 0 : for (int j = 0; j < num_points - k; j++)
821 : {
822 0 : if (k == 0)
823 : {
824 0 : auto A2j = (funcA2)(points[j].imag());
825 : ops[k].push_back(std::move(A2j));
826 : }
827 : else
828 : {
829 0 : std::complex<double> denom = points[j + k] - points[j];
830 : auto A2dd =
831 0 : BuildParSumOperator({1.0 / denom, -1.0 / denom},
832 0 : {ops[k - 1][j + 1].get(), ops[k - 1][j].get()}, true);
833 0 : ops[k].push_back(std::move(A2dd));
834 : }
835 : }
836 : }
837 :
838 : // Compute monomial coefficients as a function of the Newton polynomial coefficients.
839 : coeffs.clear();
840 0 : coeffs.assign(num_points, std::vector<std::complex<double>>(num_points, 0.0));
841 0 : for (int k = 0; k < num_points; k++)
842 : {
843 0 : for (int j = k; j < num_points; j++)
844 : {
845 0 : double sign = ((j - k) % 2 == 0) ? 1 : -1;
846 0 : coeffs[k][j] = sign * elementarySymmetric(points, j - k, j);
847 : }
848 : }
849 0 : }
850 :
851 : std::unique_ptr<ComplexOperator>
852 0 : NewtonInterpolationOperator::GetInterpolationOperator(int order) const
853 : {
854 0 : MFEM_VERIFY(order >= 0 && order < num_points,
855 : "Order must be greater than or equal to 0 and smaller than the number of "
856 : "interpolation points!");
857 0 : return BuildParSumOperator({coeffs[order][0], coeffs[order][1], coeffs[order][2]},
858 0 : {ops[0][0].get(), ops[1][0].get(), ops[2][0].get()}, true);
859 : }
860 :
861 0 : void NewtonInterpolationOperator::Mult(int order, const ComplexVector &x,
862 : ComplexVector &y) const
863 : {
864 0 : MFEM_VERIFY(order >= 0 && order < num_points,
865 : "Order must be greater than or equal to 0 and smaller than the number of "
866 : "interpolation points!");
867 :
868 : y = 0.0;
869 0 : for (int j = 0; j < num_points; j++)
870 : {
871 0 : if (coeffs[order][j] != 0.0)
872 : {
873 0 : ops[j][0]->AddMult(x, y, coeffs[order][j]);
874 : }
875 : }
876 0 : }
877 :
878 0 : void NewtonInterpolationOperator::AddMult(int order, const ComplexVector &x,
879 : ComplexVector &y, std::complex<double> a) const
880 : {
881 0 : this->Mult(order, x, rhs);
882 0 : rhs *= a;
883 : y += rhs;
884 0 : }
885 :
886 : } // namespace palace
|
