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 "romoperator.hpp"
5 :
6 : #include <Eigen/SVD>
7 : #include <mfem.hpp>
8 : #include "linalg/orthog.hpp"
9 : #include "models/spaceoperator.hpp"
10 : #include "utils/communication.hpp"
11 : #include "utils/iodata.hpp"
12 : #include "utils/timer.hpp"
13 :
14 : // Eigen does not provide a complex-valued generalized eigenvalue solver, so we use LAPACK
15 : // for this.
16 : extern "C"
17 : {
18 : void zggev_(char *, char *, int *, std::complex<double> *, int *, std::complex<double> *,
19 : int *, std::complex<double> *, std::complex<double> *, std::complex<double> *,
20 : int *, std::complex<double> *, int *, std::complex<double> *, int *, double *,
21 : int *);
22 : }
23 :
24 : namespace palace
25 : {
26 :
27 : using namespace std::complex_literals;
28 :
29 : namespace
30 : {
31 :
32 : constexpr auto ORTHOG_TOL = 1.0e-12;
33 :
34 : template <typename VecType, typename ScalarType>
35 0 : inline void OrthogonalizeColumn(Orthogonalization type, MPI_Comm comm,
36 : const std::vector<VecType> &V, VecType &w, ScalarType *Rj,
37 : int j)
38 : {
39 : // Orthogonalize w against the leading j columns of V.
40 0 : switch (type)
41 : {
42 0 : case Orthogonalization::MGS:
43 0 : linalg::OrthogonalizeColumnMGS(comm, V, w, Rj, j);
44 0 : break;
45 0 : case Orthogonalization::CGS:
46 0 : linalg::OrthogonalizeColumnCGS(comm, V, w, Rj, j);
47 0 : break;
48 0 : case Orthogonalization::CGS2:
49 0 : linalg::OrthogonalizeColumnCGS(comm, V, w, Rj, j, true);
50 0 : break;
51 : }
52 0 : }
53 :
54 0 : inline void ProjectMatInternal(MPI_Comm comm, const std::vector<Vector> &V,
55 : const ComplexOperator &A, Eigen::MatrixXcd &Ar,
56 : ComplexVector &r, int n0)
57 : {
58 : // Update Ar = Vᴴ A V for the new basis dimension n0 -> n. V is real and thus the result
59 : // is complex symmetric if A is symmetric (which we assume is the case). Ar is replicated
60 : // across all processes as a sequential n x n matrix.
61 : const auto n = Ar.rows();
62 0 : MFEM_VERIFY(n0 < n, "Invalid dimensions in PROM matrix projection!");
63 0 : for (int j = n0; j < n; j++)
64 : {
65 : // Fill block of Vᴴ A V = [ | Vᴴ A vj ] . We can optimize the matrix-vector product
66 : // since the columns of V are real.
67 0 : MFEM_VERIFY(A.Real() || A.Imag(),
68 : "Invalid zero ComplexOperator for PROM matrix projection!");
69 0 : if (A.Real())
70 : {
71 0 : A.Real()->Mult(V[j], r.Real());
72 : }
73 0 : if (A.Imag())
74 : {
75 0 : A.Imag()->Mult(V[j], r.Imag());
76 : }
77 0 : for (int i = 0; i < n; i++)
78 : {
79 0 : Ar(i, j).real(A.Real() ? V[i] * r.Real() : 0.0); // Local inner product
80 0 : Ar(i, j).imag(A.Imag() ? V[i] * r.Imag() : 0.0);
81 : }
82 : }
83 0 : Mpi::GlobalSum((n - n0) * n, Ar.data() + n0 * n, comm);
84 :
85 : // Fill lower block of Vᴴ A V = [ ____________ | ]
86 : // [ vjᴴ A V[1:n0] | ] .
87 0 : for (int j = 0; j < n0; j++)
88 : {
89 0 : for (int i = n0; i < n; i++)
90 : {
91 0 : Ar(i, j) = Ar(j, i);
92 : }
93 : }
94 0 : }
95 :
96 0 : inline void ProjectVecInternal(MPI_Comm comm, const std::vector<Vector> &V,
97 : const ComplexVector &b, Eigen::VectorXcd &br, int n0)
98 : {
99 : // Update br = Vᴴ b for the new basis dimension n0 -> n. br is replicated across all
100 : // processes as a sequential n-dimensional vector.
101 : const auto n = br.size();
102 0 : MFEM_VERIFY(n0 < n, "Invalid dimensions in PROM vector projection!");
103 0 : for (int i = n0; i < n; i++)
104 : {
105 0 : br(i).real(V[i] * b.Real()); // Local inner product
106 0 : br(i).imag(V[i] * b.Imag());
107 : }
108 0 : Mpi::GlobalSum(n - n0, br.data() + n0, comm);
109 0 : }
110 :
111 0 : inline void ComputeMRI(const Eigen::MatrixXcd &R, Eigen::VectorXcd &q)
112 : {
113 : // Compute the coefficients of the minimal rational interpolation (MRI):
114 : // u = [sum_i u_i q_i / (z - z_i)] / [sum_i q_i / (z - z_i)]. The coefficients are given
115 : // by the right singular vector of R corresponding to the minimum singular value.
116 : const auto S = R.rows();
117 : MFEM_ASSERT(S > 0 && R.cols() == S, "Invalid dimension mismatch when computing MRI!");
118 : // For Eigen = v3.4.0 (latest tagged release as of 10/2023)
119 0 : Eigen::JacobiSVD<Eigen::MatrixXcd> svd;
120 0 : svd.compute(R, Eigen::ComputeFullV);
121 : // For Eigen > v3.4.0 (GitLab repo is at v3.4.90 as of 10/2023)
122 : // Eigen::JacobiSVD<Eigen::MatrixXcd, Eigen::ComputeFullV> svd;
123 : // svd.compute(R);
124 : const auto &sigma = svd.singularValues();
125 0 : auto m = S - 1;
126 0 : while (m > 0 && sigma[m] < ORTHOG_TOL * sigma[0])
127 : {
128 0 : Mpi::Warning("Minimal rational interpolation encountered rank-deficient matrix: "
129 : "σ[{:d}] = {:.3e} (σ[0] = {:.3e})!\n",
130 : m, sigma[m], sigma[0]);
131 0 : m--;
132 : }
133 0 : q = svd.matrixV().col(m);
134 0 : }
135 :
136 : template <typename MatType, typename VecType>
137 : inline void ZGGEV(MatType &A, MatType &B, VecType &D, MatType &VR)
138 : {
139 : // Wrapper for LAPACK's (z)ggev. A and B are overwritten by their Schur decompositions.
140 : MFEM_VERIFY(A.rows() == A.cols() && B.rows() == B.cols() && A.rows() == B.rows(),
141 : "Generalized eigenvalue problem expects A, B matrices to be square and have "
142 : "same dimensions!");
143 : char jobvl = 'N', jobvr = 'V';
144 : int n = static_cast<int>(A.rows()), lwork = 2 * n;
145 : std::vector<std::complex<double>> alpha(n), beta(n), work(lwork);
146 : std::vector<double> rwork(8 * n);
147 : MatType VL(0, 0);
148 : VR.resize(n, n);
149 : int info = 0;
150 :
151 : zggev_(&jobvl, &jobvr, &n, A.data(), &n, B.data(), &n, alpha.data(), beta.data(),
152 : VL.data(), &n, VR.data(), &n, work.data(), &lwork, rwork.data(), &info);
153 : MFEM_VERIFY(info == 0, "ZGGEV failed with info = " << info << "!");
154 :
155 : // Postprocess the eigenvalues and eigenvectors (return unit 2-norm eigenvectors).
156 : D.resize(n);
157 : for (int i = 0; i < n; i++)
158 : {
159 : D(i) = (beta[i] == 0.0)
160 : ? ((alpha[i] == 0.0) ? std::numeric_limits<std::complex<double>>::quiet_NaN()
161 : : mfem::infinity())
162 : : alpha[i] / beta[i];
163 : VR.col(i) /= VR.col(i).norm();
164 : }
165 : }
166 :
167 : template <typename VecType>
168 0 : inline void ProlongatePROMSolution(std::size_t n, const std::vector<Vector> &V,
169 : const VecType &y, ComplexVector &u)
170 : {
171 : u = 0.0;
172 0 : for (std::size_t j = 0; j < n; j += 2)
173 : {
174 0 : if (j + 1 < n)
175 : {
176 0 : linalg::AXPBYPCZ(y(j).real(), V[j], y(j + 1).real(), V[j + 1], 1.0, u.Real());
177 0 : linalg::AXPBYPCZ(y(j).imag(), V[j], y(j + 1).imag(), V[j + 1], 1.0, u.Imag());
178 : }
179 : else
180 : {
181 0 : linalg::AXPY(y(j).real(), V[j], u.Real());
182 0 : linalg::AXPY(y(j).imag(), V[j], u.Imag());
183 : }
184 : }
185 0 : }
186 :
187 : } // namespace
188 :
189 0 : MinimalRationalInterpolation::MinimalRationalInterpolation(int max_size)
190 : {
191 0 : Q.resize(max_size, ComplexVector());
192 0 : }
193 :
194 0 : void MinimalRationalInterpolation::AddSolutionSample(double omega, const ComplexVector &u,
195 : const SpaceOperator &space_op,
196 : Orthogonalization orthog_type)
197 : {
198 : MPI_Comm comm = space_op.GetComm();
199 :
200 : // Compute the coefficients for the minimal rational interpolation of the state u used
201 : // as an error indicator. The complex-valued snapshot matrix U = [{u_i, (iω) u_i}] is
202 : // stored by its QR decomposition.
203 0 : MFEM_VERIFY(dim_Q + 1 <= Q.size(),
204 : "Unable to increase basis storage size, increase maximum number of vectors!");
205 0 : R.conservativeResizeLike(Eigen::MatrixXd::Zero(dim_Q + 1, dim_Q + 1));
206 : {
207 0 : std::vector<const ComplexVector *> blocks = {&u, &u};
208 0 : std::vector<std::complex<double>> s = {1.0, 1i * omega};
209 0 : Q[dim_Q].SetSize(2 * u.Size());
210 0 : Q[dim_Q].UseDevice(true);
211 0 : Q[dim_Q].SetBlocks(blocks, s);
212 : }
213 0 : OrthogonalizeColumn(orthog_type, comm, Q, Q[dim_Q], R.col(dim_Q).data(), dim_Q);
214 0 : R(dim_Q, dim_Q) = linalg::Norml2(comm, Q[dim_Q]);
215 0 : Q[dim_Q] *= 1.0 / R(dim_Q, dim_Q);
216 0 : dim_Q++;
217 0 : ComputeMRI(R, q);
218 : if constexpr (false)
219 : {
220 : Mpi::Print("MRI (S = {}):\nR = {}\nq = {}", dim_Q, R, q);
221 : }
222 0 : z.push_back(omega);
223 0 : }
224 :
225 0 : std::vector<double> MinimalRationalInterpolation::FindMaxError(int N) const
226 : {
227 : // Return an estimate for argmax_z ||u(z) - V y(z)|| as argmin_z |Q(z)| with Q(z) =
228 : // sum_i q_z / (z - z_i) (denominator of the barycentric interpolation of u). The roots of
229 : // Q are given analytically as the solution to an S + 1 dimensional eigenvalue problem.
230 0 : BlockTimer bt(Timer::CONSTRUCT_PROM);
231 0 : const auto S = dim_Q;
232 0 : MFEM_VERIFY(S >= 2, "Maximum error can only be found once two sample points have been "
233 : "added to the PROM to define the parameter domain!");
234 0 : double start = *std::min_element(z.begin(), z.end());
235 0 : double end = *std::max_element(z.begin(), z.end());
236 0 : Eigen::Map<const Eigen::VectorXd> z_map(z.data(), S);
237 0 : std::vector<std::complex<double>> z_star(N, 0.0);
238 :
239 : // XX TODO: For now, we explicitly minimize Q on the real line since we don't allow
240 : // samples at complex-valued points (yet).
241 :
242 : // Eigen::MatrixXcd A = Eigen::MatrixXcd::Zero(S + 1, S + 1);
243 : // A.diagonal().head(S) = z_map.array();
244 : // A.row(S).head(S) = q;
245 : // A.col(S).head(S) = Eigen::VectorXcd::Ones(S);
246 :
247 : // Eigen::MatrixXcd B = Eigen::MatrixXcd::Identity(S + 1, S + 1);
248 : // B(S, S) = 0.0;
249 :
250 : // Eigen::VectorXcd D;
251 : // Eigen::MatrixXcd X;
252 : // ZGGEV(A, B, D, X);
253 :
254 : // // If there are multiple roots in [start, end], pick the ones furthest from the
255 : // // existing set of samples.
256 : // {
257 : // std::vector<double> dist_star(N, 0.0);
258 : // for (auto d : D)
259 : // {
260 : // if (std::real(d) < start || std::real(d) > end)
261 : // {
262 : // continue;
263 : // }
264 : // const double dist = (z_map.array() - std::real(d)).abs().maxCoeff();
265 : // for (int i = 0; i < N; i++)
266 : // {
267 : // if (dist > dist_star[i])
268 : // {
269 : // for (int j = i + 1; j < N; j++)
270 : // {
271 : // z_star[j] = z_star[j - 1];
272 : // dist_star[j] = dist_star[j - 1];
273 : // }
274 : // z_star[i] = start;
275 : // dist_star[i] = dist;
276 : // }
277 : // }
278 : // }
279 : // }
280 :
281 : // Fall back to sampling Q on discrete points if no roots exist in [start, end].
282 0 : if (std::abs(z_star[0]) == 0.0)
283 : {
284 0 : const auto delta = (end - start) / 1.0e6;
285 0 : std::vector<double> Q_star(N, mfem::infinity());
286 0 : while (start <= end)
287 : {
288 0 : const double Q = std::abs((q.array() / (z_map.array() - start)).sum());
289 0 : for (int i = 0; i < N; i++)
290 : {
291 0 : if (Q < Q_star[i])
292 : {
293 0 : for (int j = i + 1; j < N; j++)
294 : {
295 0 : z_star[j] = z_star[j - 1];
296 0 : Q_star[j] = Q_star[j - 1];
297 : }
298 : z_star[i] = start;
299 0 : Q_star[i] = Q;
300 : }
301 : }
302 0 : start += delta;
303 : }
304 0 : MFEM_VERIFY(
305 : N == 0 || std::abs(z_star[0]) > 0.0,
306 : fmt::format("Could not locate a maximum error in the range [{}, {}]!", start, end));
307 : }
308 0 : std::vector<double> vals(z_star.size());
309 : std::transform(z_star.begin(), z_star.end(), vals.begin(),
310 : [](std::complex<double> z) { return std::real(z); });
311 0 : return vals;
312 0 : }
313 :
314 0 : RomOperator::RomOperator(const IoData &iodata, SpaceOperator &space_op,
315 0 : int max_size_per_excitation)
316 0 : : space_op(space_op), orthog_type(iodata.solver.linear.gs_orthog)
317 : {
318 : // Construct the system matrices defining the linear operator. PEC boundaries are handled
319 : // simply by setting diagonal entries of the system matrix for the corresponding dofs.
320 : // Because the Dirichlet BC is always homogeneous, no special elimination is required on
321 : // the RHS. The damping matrix may be nullptr.
322 0 : K = space_op.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
323 0 : C = space_op.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
324 0 : M = space_op.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
325 0 : MFEM_VERIFY(K && M, "Invalid empty HDM matrices when constructing PROM!");
326 :
327 : // Initialize working vector storage.
328 0 : r.SetSize(K->Height());
329 0 : r.UseDevice(true);
330 :
331 : // Set up the linear solver and set operators but don't set the operators yet (this will
332 : // be done during an HDM solve at a given parameter point). The preconditioner for the
333 : // complex linear system is constructed from a real approximation to the complex system
334 : // matrix.
335 0 : ksp = std::make_unique<ComplexKspSolver>(iodata, space_op.GetNDSpaces(),
336 0 : &space_op.GetH1Spaces());
337 :
338 : auto excitation_helper = space_op.GetPortExcitations();
339 :
340 : // The initial PROM basis is empty. The provided maximum dimension is the number of sample
341 : // points (2 basis vectors per point). Basis orthogonalization method is configured using
342 : // GMRES/FGMRES settings.
343 0 : MFEM_VERIFY(max_size_per_excitation * excitation_helper.Size() > 0,
344 : "Reduced order basis storage must have > 0 columns!");
345 0 : V.resize(2 * max_size_per_excitation * excitation_helper.Size(), Vector());
346 :
347 : // Set up MRI.
348 0 : for (const auto &[excitation_idx, data] : excitation_helper)
349 : {
350 0 : mri.emplace(excitation_idx, MinimalRationalInterpolation(max_size_per_excitation));
351 : }
352 0 : }
353 :
354 0 : void RomOperator::SetExcitationIndex(int excitation_idx)
355 : {
356 : // Set up RHS vector (linear in frequency part) for the incident field at port boundaries,
357 : // and the vector for the solution, which satisfies the Dirichlet (PEC) BC.
358 0 : excitation_idx_cache = excitation_idx;
359 0 : has_RHS1 = space_op.GetExcitationVector1(excitation_idx_cache, RHS1);
360 0 : if (!has_RHS1)
361 : {
362 0 : RHS1.SetSize(0);
363 : }
364 : else
365 : {
366 : // Project RHS1 to RHS1r with current PROM.
367 0 : if (dim_V > 0)
368 : {
369 0 : auto comm = space_op.GetComm();
370 0 : RHS1r.conservativeResize(dim_V);
371 0 : ProjectVecInternal(comm, V, RHS1, RHS1r, 0);
372 : }
373 : }
374 0 : }
375 :
376 0 : void RomOperator::SolveHDM(int excitation_idx, double omega, ComplexVector &u)
377 : {
378 0 : if (excitation_idx_cache != excitation_idx)
379 : {
380 0 : SetExcitationIndex(excitation_idx);
381 : }
382 : // Compute HDM solution at the given frequency. The system matrix, A = K + iω C - ω² M +
383 : // A2(ω) is built by summing the underlying operator contributions.
384 0 : A2 = space_op.GetExtraSystemMatrix<ComplexOperator>(omega, Operator::DIAG_ZERO);
385 0 : has_A2 = (A2 != nullptr);
386 0 : auto A = space_op.GetSystemMatrix(std::complex<double>(1.0, 0.0), 1i * omega,
387 0 : std::complex<double>(-omega * omega, 0.0), K.get(),
388 0 : C.get(), M.get(), A2.get());
389 0 : auto P = space_op.GetPreconditionerMatrix<ComplexOperator>(1.0 + 0.0i, 1i * omega,
390 0 : -omega * omega + 0.0i, omega);
391 0 : ksp->SetOperators(*A, *P);
392 :
393 : // The HDM excitation vector is computed as RHS = iω RHS1 + RHS2(ω).
394 0 : Mpi::Print("\n");
395 0 : if (has_RHS2)
396 : {
397 0 : has_RHS2 = space_op.GetExcitationVector2(excitation_idx, omega, r);
398 : }
399 : else
400 : {
401 0 : r = 0.0;
402 : }
403 0 : if (has_RHS1)
404 : {
405 0 : r.Add(1i * omega, RHS1);
406 : }
407 :
408 : // Solve the linear system.
409 0 : ksp->Mult(r, u);
410 0 : }
411 :
412 0 : void RomOperator::UpdatePROM(const ComplexVector &u)
413 : {
414 :
415 : // Update V. The basis is always real (each complex solution adds two basis vectors if it
416 : // has a nonzero real and imaginary parts).
417 0 : BlockTimer bt(Timer::CONSTRUCT_PROM);
418 0 : MPI_Comm comm = space_op.GetComm();
419 0 : const double normr = linalg::Norml2(comm, u.Real());
420 0 : const double normi = linalg::Norml2(comm, u.Imag());
421 0 : const bool has_real = (normr > ORTHOG_TOL * std::sqrt(normr * normr + normi * normi));
422 0 : const bool has_imag = (normi > ORTHOG_TOL * std::sqrt(normr * normr + normi * normi));
423 0 : MFEM_VERIFY(dim_V + has_real + has_imag <= V.size(),
424 : "Unable to increase basis storage size, increase maximum number of vectors!");
425 : const std::size_t dim_V0 = dim_V;
426 : std::vector<double> H(dim_V + static_cast<std::size_t>(has_real) +
427 0 : static_cast<std::size_t>(has_imag));
428 0 : if (has_real)
429 : {
430 0 : V[dim_V] = u.Real();
431 0 : OrthogonalizeColumn(orthog_type, comm, V, V[dim_V], H.data(), dim_V);
432 0 : H[dim_V] = linalg::Norml2(comm, V[dim_V]);
433 0 : V[dim_V] *= 1.0 / H[dim_V];
434 0 : dim_V++;
435 : }
436 0 : if (has_imag)
437 : {
438 0 : V[dim_V] = u.Imag();
439 0 : OrthogonalizeColumn(orthog_type, comm, V, V[dim_V], H.data(), dim_V);
440 0 : H[dim_V] = linalg::Norml2(comm, V[dim_V]);
441 0 : V[dim_V] *= 1.0 / H[dim_V];
442 0 : dim_V++;
443 : }
444 :
445 : // Update reduced-order operators. Resize preserves the upper dim0 x dim0 block of each
446 : // matrix and first dim0 entries of each vector and the projection uses the values
447 : // computed for the unchanged basis vectors.
448 0 : Kr.conservativeResize(dim_V, dim_V);
449 0 : ProjectMatInternal(comm, V, *K, Kr, r, dim_V0);
450 0 : if (C)
451 : {
452 0 : Cr.conservativeResize(dim_V, dim_V);
453 0 : ProjectMatInternal(comm, V, *C, Cr, r, dim_V0);
454 : }
455 0 : Mr.conservativeResize(dim_V, dim_V);
456 0 : ProjectMatInternal(comm, V, *M, Mr, r, dim_V0);
457 0 : Ar.resize(dim_V, dim_V);
458 0 : if (RHS1.Size())
459 : {
460 0 : RHS1r.conservativeResize(dim_V);
461 0 : ProjectVecInternal(comm, V, RHS1, RHS1r, dim_V0);
462 : }
463 0 : RHSr.resize(dim_V);
464 0 : }
465 :
466 0 : void RomOperator::UpdateMRI(int excitation_idx, double omega, const ComplexVector &u)
467 : {
468 0 : BlockTimer bt(Timer::CONSTRUCT_PROM);
469 0 : mri.at(excitation_idx).AddSolutionSample(omega, u, space_op, orthog_type);
470 0 : }
471 :
472 0 : void RomOperator::SolvePROM(int excitation_idx, double omega, ComplexVector &u)
473 : {
474 0 : if (excitation_idx_cache != excitation_idx)
475 : {
476 0 : SetExcitationIndex(excitation_idx);
477 : }
478 :
479 : // Assemble the PROM linear system at the given frequency. The PROM system is defined by
480 : // the matrix Aᵣ(ω) = Kᵣ + iω Cᵣ - ω² Mᵣ + Vᴴ A2 V(ω) and source vector RHSᵣ(ω) =
481 : // iω RHS1ᵣ + Vᴴ RHS2(ω). A2(ω) and RHS2(ω) are constructed only if required and are
482 : // only nonzero on boundaries, will be empty if not needed.
483 0 : if (has_A2 && Ar.rows() > 0)
484 : {
485 0 : A2 = space_op.GetExtraSystemMatrix<ComplexOperator>(omega, Operator::DIAG_ZERO);
486 0 : ProjectMatInternal(space_op.GetComm(), V, *A2, Ar, r, 0);
487 : }
488 : else
489 : {
490 : Ar.setZero();
491 : }
492 0 : Ar += Kr;
493 0 : if (C)
494 : {
495 0 : Ar += (1i * omega) * Cr;
496 : }
497 0 : Ar += (-omega * omega) * Mr;
498 :
499 0 : if (has_RHS2 && RHSr.size() > 0)
500 : {
501 0 : space_op.GetExcitationVector2(excitation_idx, omega, RHS2);
502 0 : ProjectVecInternal(space_op.GetComm(), V, RHS2, RHSr, 0);
503 : }
504 : else
505 : {
506 : RHSr.setZero();
507 : }
508 0 : if (has_RHS1)
509 : {
510 0 : RHSr += (1i * omega) * RHS1r;
511 : }
512 :
513 : // Compute PROM solution at the given frequency and expand into high-dimensional space.
514 : // The PROM is solved on every process so the matrix-vector product for vector expansion
515 : // does not require communication.
516 0 : BlockTimer bt(Timer::SOLVE_PROM);
517 : if constexpr (false)
518 : {
519 : // LDLT solve.
520 : RHSr = Ar.ldlt().solve(RHSr);
521 : RHSr = Ar.selfadjointView<Eigen::Lower>().ldlt().solve(RHSr);
522 : }
523 : else
524 : {
525 : // LU solve.
526 0 : RHSr = Ar.partialPivLu().solve(RHSr);
527 : }
528 0 : ProlongatePROMSolution(dim_V, V, RHSr, u);
529 0 : }
530 :
531 0 : std::vector<std::complex<double>> RomOperator::ComputeEigenvalueEstimates() const
532 : {
533 : // XX TODO: Not yet implemented
534 0 : MFEM_ABORT("Eigenvalue estimates for PROM operators are not yet implemented!");
535 : return {};
536 : }
537 :
538 : } // namespace palace
|
