Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #ifndef PALACE_LINALG_NLEPS_HPP
5 : #define PALACE_LINALG_NLEPS_HPP
6 :
7 : #include <complex>
8 : #include <memory>
9 : #include <optional>
10 : #include <mpi.h>
11 : #include "linalg/eps.hpp"
12 : #include "linalg/ksp.hpp"
13 : #include "linalg/operator.hpp"
14 : #include "linalg/vector.hpp"
15 : #include "models/spaceoperator.hpp"
16 :
17 : namespace palace
18 : {
19 :
20 : //
21 : // Abstract base class for nonlinear eigensolvers.
22 : // Currently only implemented for complex scalar interface.
23 : //
24 : class NonLinearEigenvalueSolver : public EigenvalueSolver
25 : {
26 : protected:
27 : // MPI communicator.
28 : MPI_Comm comm;
29 :
30 : // Control print level for debugging.
31 : int print;
32 :
33 : // Number of eigenvalues to compute and problem size.
34 : int nev, n;
35 :
36 : // Relative eigenvalue error convergence tolerance for the solver.
37 : double rtol;
38 :
39 : // Maximum number of iterations.
40 : int nleps_it;
41 :
42 : // Specifies which part of the spectrum to search for.
43 : EigenvalueSolver::WhichType which_type;
44 :
45 : // Variables for scaling, from Higham et al., IJNME 2008.
46 : double gamma, delta;
47 :
48 : // Parameters defining the spectral transformation.
49 : std::complex<double> sigma;
50 : bool sinvert;
51 :
52 : // Storage for computed eigenvalues and eigenvectors.
53 : std::vector<std::complex<double>> eigenvalues;
54 : std::vector<ComplexVector> eigenvectors;
55 : std::unique_ptr<int[]> perm;
56 :
57 : // Storage for computed residual norms and eigenvector scalings.
58 : std::unique_ptr<double[]> res, xscale;
59 :
60 : // Reference to linear solver used for operator action of P(σ)⁻¹ (not owned).
61 : ComplexKspSolver *opInv;
62 :
63 : // Reference to solver for projecting an intermediate vector onto a divergence-free space
64 : // (not owned).
65 : const DivFreeSolver<ComplexVector> *opProj;
66 :
67 : // Reference to matrix used for weighted inner products (not owned). May be nullptr, in
68 : // which case identity is used.
69 : const Operator *opB;
70 :
71 : // Workspace vector for operator applications.
72 : mutable ComplexVector x1, y1;
73 :
74 : // Helper routine for computing the eigenvector normalization.
75 : double GetEigenvectorNorm(const ComplexVector &x, ComplexVector &Bx) const;
76 :
77 : // Helper routine for computing the eigenpair residual.
78 : virtual double GetResidualNorm(std::complex<double> l, const ComplexVector &x,
79 : ComplexVector &r) const = 0;
80 :
81 : // Helper routine for computing the backward error.
82 : virtual double GetBackwardScaling(std::complex<double> l) const = 0;
83 :
84 : // Get the associated MPI communicator.
85 : virtual MPI_Comm GetComm() const = 0;
86 :
87 : // Return problem type name.
88 : virtual const char *GetName() const = 0;
89 :
90 : public:
91 : NonLinearEigenvalueSolver(MPI_Comm comm, int print);
92 :
93 : // The linear solver will be configured to compute the action of T(σ)⁻¹
94 : // where σ is the current eigenvalue estimate.
95 : void SetLinearSolver(ComplexKspSolver &ksp) override;
96 :
97 : // Set the projection operator for enforcing the divergence-free constraint.
98 : void SetDivFreeProjector(const DivFreeSolver<ComplexVector> &divfree) override;
99 :
100 : // Set optional B matrix used for weighted inner products. This must be set explicitly
101 : // even for generalized problems, otherwise the identity will be used.
102 : void SetBMat(const Operator &B) override;
103 :
104 : // Get scaling factors used by the solver.
105 0 : double GetScalingGamma() const override { return gamma; }
106 0 : double GetScalingDelta() const override { return delta; }
107 :
108 : // Set the number of required eigenmodes.
109 : void SetNumModes(int num_eig, int num_vec = 0) override;
110 :
111 : // Set solver tolerance.
112 : void SetTol(double tol) override;
113 :
114 : // Set maximum number of Arnoldi update iterations.
115 : void SetMaxIter(int max_it) override;
116 :
117 : // Set target spectrum for the eigensolver. When a spectral transformation is used, this
118 : // applies to the spectrum of the shifted operator.
119 : void SetWhichEigenpairs(WhichType type) override;
120 :
121 : // Set shift-and-invert spectral transformation.
122 : void SetShiftInvert(std::complex<double> s, bool precond = false) override;
123 :
124 : // Set an initial vector for the solution subspace.
125 : void SetInitialSpace(const ComplexVector &v) override;
126 :
127 : // Solve the eigenvalue problem. Returns the number of converged eigenvalues.
128 : int Solve() override = 0;
129 :
130 : // Get the corresponding eigenvalue.
131 : std::complex<double> GetEigenvalue(int i) const override;
132 :
133 : // Get the corresponding eigenvector. Eigenvectors are normalized such that ||x||₂ = 1,
134 : // unless the B-matrix is set for weighted inner products.
135 : void GetEigenvector(int i, ComplexVector &x) const override;
136 :
137 : // Get the corresponding eigenpair error.
138 : double GetError(int i, ErrorType type) const override;
139 :
140 : // Re-normalize the given number of eigenvectors, for example if the matrix B for weighted
141 : // inner products has changed. This does not perform re-orthogonalization with respect to
142 : // the new matrix, only normalization.
143 : void RescaleEigenvectors(int num_eig) override;
144 : };
145 :
146 : // Quasi-Newton nonlinear eigenvalue solver for (K + λ C + λ² M + A2(λ)) x = 0.
147 : class QuasiNewtonSolver : public NonLinearEigenvalueSolver
148 : {
149 : private:
150 : // References to matrices defining the nonlinear eigenvalue problem
151 : // (not owned).
152 : const ComplexOperator *opK, *opC, *opM;
153 :
154 : // Operators used in the iterative linear solver.
155 : std::unique_ptr<ComplexOperator> opA2, opA, opP;
156 :
157 : // Function to compute the A2 operator.
158 : std::optional<std::function<std::unique_ptr<ComplexOperator>(double)>> funcA2;
159 :
160 : // Function to compute the preconditioner matrix.
161 : std::optional<std::function<std::unique_ptr<ComplexOperator>(
162 : std::complex<double>, std::complex<double>, std::complex<double>, double)>>
163 : funcP;
164 :
165 : // Linear eigenvalue solver used to set initial guess.
166 : std::unique_ptr<EigenvalueSolver> linear_eigensolver_;
167 :
168 : // Number of eigenmode initial guesses.
169 : int nev_linear;
170 :
171 : // Operator norms for scaling.
172 : mutable double normK, normC, normM;
173 :
174 : // Update frequency of the preconditioner during Newton iterations.
175 : int preconditioner_lag;
176 :
177 : // Update tolerance of the preconditioner (no update below tol).
178 : double preconditioner_tol;
179 :
180 : // Maximum number of Newton attempts with the same initial guess.
181 : int max_restart;
182 :
183 : // Refine linear eigenvalues with nonlinear Newton eigenvalue solver.
184 : bool refine_nonlinear;
185 :
186 : // Set the initial guesses from the linear eigenvalue solver results.
187 : void SetInitialGuess();
188 :
189 : protected:
190 : double GetResidualNorm(std::complex<double> l, const ComplexVector &x,
191 : ComplexVector &r) const override;
192 :
193 : double GetBackwardScaling(std::complex<double> l) const override;
194 :
195 0 : const char *GetName() const override { return "QuasiNewton"; }
196 :
197 : public:
198 : QuasiNewtonSolver(MPI_Comm comm, std::unique_ptr<EigenvalueSolver> linear_eigensolver,
199 : int num_conv, int print, bool refine);
200 :
201 : using NonLinearEigenvalueSolver::SetOperators;
202 : void SetOperators(const ComplexOperator &K, const ComplexOperator &M,
203 : ScaleType type) override;
204 : void SetOperators(const ComplexOperator &K, const ComplexOperator &C,
205 : const ComplexOperator &M, ScaleType type) override;
206 :
207 : // Set the frequency-dependent A2 matrix function.
208 : void SetExtraSystemMatrix(
209 : std::function<std::unique_ptr<ComplexOperator>(double)>) override;
210 :
211 : // Set the preconditioner update function.
212 : void SetPreconditionerUpdate(std::function<std::unique_ptr<ComplexOperator>(
213 : std::complex<double>, std::complex<double>,
214 : std::complex<double>, double)>) override;
215 :
216 : // Set the update frequency of the preconditioner.
217 : void SetPreconditionerLag(int preconditioner_update_freq,
218 : double preconditioner_update_tol);
219 :
220 : // Set the maximum number of restarts with the same initial guess.
221 : void SetMaxRestart(int max_num_restart);
222 :
223 : // Solve the nonlinear eigenvalue problem.
224 : int Solve() override;
225 :
226 0 : MPI_Comm GetComm() const override { return comm; }
227 : };
228 :
229 : //
230 : // Interpolation operators to approximate the nonlinear A2 operator.
231 : //
232 : class Interpolation
233 : {
234 : public:
235 : Interpolation() = default;
236 : virtual ~Interpolation() = default;
237 : virtual void Interpolate(const std::complex<double> sigma_min,
238 : const std::complex<double> sigma_max) = 0;
239 : virtual std::unique_ptr<ComplexOperator> GetInterpolationOperator(int order) const = 0;
240 : virtual void Mult(int order, const ComplexVector &x, ComplexVector &y) const = 0;
241 : virtual void AddMult(int order, const ComplexVector &x, ComplexVector &y,
242 : std::complex<double> a = 1.0) const = 0;
243 : };
244 :
245 : // Newton polynomial interpolation to approximate the nonlinear A2 operator.
246 : class NewtonInterpolationOperator : public Interpolation
247 : {
248 : private:
249 : // Function to compute the A2 operator.
250 : std::function<std::unique_ptr<ComplexOperator>(double)> funcA2;
251 :
252 : // Number of points used in the interpolation (currently always second order).
253 : int num_points = 3;
254 :
255 : // Interpolation points.
256 : std::vector<std::complex<double>> points;
257 :
258 : // Monomial basis coefficients.
259 : std::vector<std::vector<std::complex<double>>> coeffs;
260 :
261 : // Divided difference operators.
262 : std::vector<std::vector<std::unique_ptr<ComplexOperator>>> ops;
263 :
264 : // Workspace objects for solver application.
265 : mutable ComplexVector rhs;
266 :
267 : public:
268 : NewtonInterpolationOperator(
269 : std::function<std::unique_ptr<ComplexOperator>(double)> funcA2, const int size);
270 :
271 : // Interpolate the A2 matrix between sigma_min and sigma_max with a Newton polynomial.
272 : void Interpolate(const std::complex<double> sigma_min,
273 : const std::complex<double> sigma_max);
274 :
275 : // Get the interpolation operator of specified order.
276 : std::unique_ptr<ComplexOperator> GetInterpolationOperator(int order) const;
277 :
278 : // Perform multiplication with interpolation operator of specified order.
279 : void Mult(int order, const ComplexVector &x, ComplexVector &y) const;
280 : void AddMult(int order, const ComplexVector &x, ComplexVector &y,
281 : std::complex<double> a = 1.0) const;
282 : };
283 :
284 : } // namespace palace
285 :
286 : #endif // PALACE_LINALG_NLEPS_HPP
|
