include/pops/numerics/linalg/dense_eig.hpp Source FileΒΆ

adc_cpp: include/pops/numerics/linalg/dense_eig.hpp Source File
adc_cpp 0.3.0
Model-free C++23 core for coupled hyperbolic-elliptic systems on adaptive (AMR) meshes, with MPI and GPU (Kokkos) backends
dense_eig.hpp
Go to the documentation of this file.
1#pragma once
42
43#include <cmath>
44#include <limits>
45
46#include <pops/core/foundation/types.hpp>
47
48namespace pops {
49
53struct EigBounds {
54 Real lmin;
55 Real lmax;
56 Real max_im;
59 bool converged;
60
76 POPS_HD bool all_real(Real im_tol = Real(1e-5)) const {
77 const Real rho = std::fabs(lmin) > std::fabs(lmax) ? std::fabs(lmin) : std::fabs(lmax);
78 const Real scale = rho > Real(1) ? rho : Real(1);
79 return converged && max_im <= im_tol * scale;
80 }
84 POPS_HD bool has_complex_pair(Real im_tol = Real(1e-5)) const {
85 return converged && !all_real(im_tol);
86 }
87};
88
93enum class Spectrum : int { kReal = 0, kComplexPair = 1, kUnknown = 2 };
94
95namespace detail {
96
100template <int N>
101POPS_HD inline void gershgorin_bounds(const Real (&A)[N][N], Real& lo, Real& hi) {
102 for (int i = 0; i < N; ++i) {
103 Real r = Real(0);
104 for (int j = 0; j < N; ++j)
105 if (j != i)
106 r += std::fabs(A[i][j]);
107 const Real l = A[i][i] - r, h = A[i][i] + r;
108 if (i == 0 || l < lo)
109 lo = l;
110 if (i == 0 || h > hi)
111 hi = h;
112 }
113}
114
117template <int N>
118POPS_HD inline void hessenberg_reduce(Real (&H)[N][N]) {
119 Real v[N]; // Householder vector of the current step (components k..N-1)
120 for (int k = 1; k <= N - 2; ++k) {
121 Real scale = Real(0);
122 for (int i = k; i < N; ++i)
123 scale += std::fabs(H[i][k - 1]);
124 if (scale == Real(0))
125 continue; // column already zero below the subdiagonal
126 Real h = Real(0);
127 for (int i = k; i < N; ++i) {
128 v[i] = H[i][k - 1] / scale;
129 h += v[i] * v[i];
130 }
131 Real g = std::sqrt(h);
132 if (v[k] > Real(0))
133 g = -g;
134 h -= v[k] * g; // h = v.v / 2 after updating v[k]
135 v[k] -= g;
136 if (h == Real(0))
137 continue;
138 // P = I - v v^T / h; H <- P H P (column k-1 set explicitly, exact zeros below)
139 for (int j = k; j < N; ++j) { // P * H on rows k..N-1
140 Real f = Real(0);
141 for (int i = k; i < N; ++i)
142 f += v[i] * H[i][j];
143 f /= h;
144 for (int i = k; i < N; ++i)
145 H[i][j] -= f * v[i];
146 }
147 for (int i = 0; i < N; ++i) { // H * P on columns k..N-1
148 Real f = Real(0);
149 for (int j = k; j < N; ++j)
150 f += H[i][j] * v[j];
151 f /= h;
152 for (int j = k; j < N; ++j)
153 H[i][j] -= f * v[j];
154 }
155 H[k][k - 1] = scale * g;
156 for (int i = k + 1; i < N; ++i)
157 H[i][k - 1] = Real(0);
158 }
159}
160
161POPS_HD inline Real hqr_copysign(Real mag, Real sgn) {
162 return sgn >= Real(0) ? std::fabs(mag) : -std::fabs(mag);
163}
164
167POPS_HD inline void record_eig(Real re, Real im, Real& lmin, Real& lmax, Real& max_im, bool& first) {
168 if (first || re < lmin)
169 lmin = re;
170 if (first || re > lmax)
171 lmax = re;
172 const Real ai = std::fabs(im);
173 if (first || ai > max_im)
174 max_im = ai;
175 first = false;
176}
177
182template <int N>
183POPS_HD inline bool hqr_minmax(Real (&H)[N][N], Real& lmin, Real& lmax, Real& max_im,
184 int max_iter_per_eig) {
185 constexpr Real kEps = std::numeric_limits<Real>::epsilon(); // follows the Real type
186 Real anorm = Real(0); // norm of the Hessenberg part (deflation criterion for the s == 0 cases)
187 for (int i = 0; i < N; ++i)
188 for (int j = (i > 0 ? i - 1 : 0); j < N; ++j)
189 anorm += std::fabs(H[i][j]);
190 if (anorm == Real(0)) { // null matrix: spectrum {0}
191 lmin = lmax = max_im = Real(0);
192 return true;
193 }
194
195 bool first = true;
196
197 int nn = N - 1; // top index of the active block
198 Real t = Real(0); // cumulative shift (exceptional shifts)
199 Real p = Real(0), q = Real(0), r = Real(0), x, y, z, w, s;
200 while (nn >= 0) {
201 int its = 0;
202 int l;
203 do {
204 // deflation: smallest l such that H[l][l-1] is negligible
205 for (l = nn; l >= 1; --l) {
206 s = std::fabs(H[l - 1][l - 1]) + std::fabs(H[l][l]);
207 if (s == Real(0))
208 s = anorm;
209 if (std::fabs(H[l][l - 1]) <= kEps * s) {
210 H[l][l - 1] = Real(0);
211 break;
212 }
213 }
214 x = H[nn][nn];
215 if (l == nn) { // real 1x1 eigenvalue
216 record_eig(x + t, Real(0), lmin, lmax, max_im, first);
217 --nn;
218 } else {
219 y = H[nn - 1][nn - 1];
220 w = H[nn][nn - 1] * H[nn - 1][nn];
221 if (l == nn - 1) { // 2x2 block: real pair or complex conjugate pair
222 p = Real(0.5) * (y - x);
223 q = p * p + w;
224 z = std::sqrt(std::fabs(q));
225 x += t;
226 if (q >= Real(0)) { // two real values
227 z = p + hqr_copysign(z, p);
228 record_eig(x + z, Real(0), lmin, lmax, max_im, first);
229 record_eig(z != Real(0) ? x - w / z : x + z, Real(0), lmin, lmax, max_im, first);
230 } else { // complex pair: Re = x + p, |Im| = z
231 record_eig(x + p, z, lmin, lmax, max_im, first);
232 record_eig(x + p, -z, lmin, lmax, max_im, first);
233 }
234 nn -= 2;
235 } else { // block > 2: one Francis double-shift iteration
236 if (its == max_iter_per_eig)
237 return false; // cap reached -> caller fallback
238 if (its == 10 || its == 20) { // exceptional shift (slow cycles)
239 t += x;
240 for (int i = 0; i <= nn; ++i)
241 H[i][i] -= x;
242 s = std::fabs(H[nn][nn - 1]) + std::fabs(H[nn - 1][nn - 2]);
243 y = x = Real(0.75) * s;
244 w = Real(-0.4375) * s * s;
245 }
246 ++its;
247 int m;
248 for (m = nn - 2; m >= l; --m) { // two consecutive small subdiagonals
249 z = H[m][m];
250 r = x - z;
251 s = y - z;
252 p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
253 q = H[m + 1][m + 1] - z - r - s;
254 r = H[m + 2][m + 1];
255 s = std::fabs(p) + std::fabs(q) + std::fabs(r);
256 p /= s;
257 q /= s;
258 r /= s;
259 if (m == l)
260 break;
261 const Real u = std::fabs(H[m][m - 1]) * (std::fabs(q) + std::fabs(r));
262 const Real v = std::fabs(p) *
263 (std::fabs(H[m - 1][m - 1]) + std::fabs(z) + std::fabs(H[m + 1][m + 1]));
264 if (u <= kEps * v)
265 break;
266 }
267 for (int i = m + 2; i <= nn; ++i) {
268 H[i][i - 2] = Real(0);
269 if (i > m + 2)
270 H[i][i - 3] = Real(0);
271 }
272 for (int k = m; k <= nn - 1; ++k) { // QR double-shift sweep on columns m..nn-1
273 if (k != m) {
274 p = H[k][k - 1];
275 q = H[k + 1][k - 1];
276 r = (k != nn - 1) ? H[k + 2][k - 1] : Real(0);
277 x = std::fabs(p) + std::fabs(q) + std::fabs(r);
278 if (x != Real(0)) {
279 p /= x;
280 q /= x;
281 r /= x;
282 }
283 }
284 s = hqr_copysign(std::sqrt(p * p + q * q + r * r), p);
285 if (s == Real(0))
286 continue;
287 if (k == m) {
288 if (l != m)
289 H[k][k - 1] = -H[k][k - 1];
290 } else {
291 H[k][k - 1] = -s * x;
292 }
293 p += s;
294 x = p / s;
295 y = q / s;
296 z = r / s;
297 q /= p;
298 r /= p;
299 for (int j = k; j <= nn; ++j) { // transform rows k..k+2
300 p = H[k][j] + q * H[k + 1][j];
301 if (k != nn - 1) {
302 p += r * H[k + 2][j];
303 H[k + 2][j] -= p * z;
304 }
305 H[k + 1][j] -= p * y;
306 H[k][j] -= p * x;
307 }
308 const int mmin = (nn < k + 3) ? nn : k + 3;
309 for (int i = l; i <= mmin; ++i) { // transform columns k..k+2
310 p = x * H[i][k] + y * H[i][k + 1];
311 if (k != nn - 1) {
312 p += z * H[i][k + 2];
313 H[i][k + 2] -= p * r;
314 }
315 H[i][k + 1] -= p * q;
316 H[i][k] -= p;
317 }
318 }
319 }
320 }
321 } while (l < nn - 1);
322 }
323 return true;
324}
325
326} // namespace detail
327
342template <int N>
343POPS_HD inline EigBounds real_eig_minmax(const Real (&A)[N][N], int max_iter_per_eig = 100,
344 bool* fallback = nullptr) {
345 static_assert(N >= 1, "real_eig_minmax: N >= 1");
346 static_assert(N <= 16,
347 "real_eig_minmax: block limited to 16x16 (stack buffer O(N^2) per device "
348 "thread, ~2 KB; beyond that, a dense solver with allocation is more "
349 "appropriate than this path)");
350 EigBounds b{Real(0), Real(0), Real(0), true};
351 if constexpr (N == 1) {
352 b.lmin = b.lmax = A[0][0];
353 } else if constexpr (N == 2) { // closed form: trace / determinant
354 const Real tr2 = Real(0.5) * (A[0][0] + A[1][1]);
355 const Real disc = Real(0.25) * (A[0][0] - A[1][1]) * (A[0][0] - A[1][1]) + A[0][1] * A[1][0];
356 if (disc >= Real(0)) {
357 const Real z = std::sqrt(disc);
358 b.lmin = tr2 - z;
359 b.lmax = tr2 + z;
360 } else {
361 b.lmin = b.lmax = tr2;
362 b.max_im = std::sqrt(-disc);
363 }
364 } else {
365 Real H[N][N]; // working copy (A is not modified)
366 for (int i = 0; i < N; ++i)
367 for (int j = 0; j < N; ++j)
368 H[i][j] = A[i][j];
369 detail::hessenberg_reduce(H);
370 if (!detail::hqr_minmax(H, b.lmin, b.lmax, b.max_im, max_iter_per_eig)) {
371 // non-convergence: external Gershgorin bound on the ORIGINAL matrix (the working copy is in
372 // an intermediate state) -- safe bound, not the spectrum. max_im forced to 0 by CONVENTION
373 // (nothing was computed), never to be interpreted as a real spectrum.
374 detail::gershgorin_bounds(A, b.lmin, b.lmax);
375 b.max_im = Real(0);
376 b.converged = false;
377 }
378 }
379 if (fallback)
380 *fallback = !b.converged;
381 return b;
382}
383
395template <int N>
396POPS_HD inline Spectrum real_spectrum(const Real (&A)[N][N], Real im_tol = Real(1e-5),
397 int max_iter_per_eig = 100) {
398 const EigBounds b = real_eig_minmax(A, max_iter_per_eig);
399 if (!b.converged)
400 return Spectrum::kUnknown; // non-convergence is EXPLICIT, never kReal
401 return b.all_real(im_tol) ? Spectrum::kReal : Spectrum::kComplexPair;
402}
403
404namespace detail {
405
409template <int N>
410POPS_HD inline bool mat_inverse(const Real (&A)[N][N], Real (&inv)[N][N],
411 Real pivot_tol = Real(1e-300)) {
412 Real M[N][N];
413 for (int i = 0; i < N; ++i)
414 for (int j = 0; j < N; ++j) {
415 M[i][j] = A[i][j];
416 inv[i][j] = (i == j) ? Real(1) : Real(0);
417 }
418 for (int col = 0; col < N; ++col) {
419 int piv = col;
420 Real best = std::fabs(M[col][col]);
421 for (int r = col + 1; r < N; ++r) {
422 const Real v = std::fabs(M[r][col]);
423 if (v > best) {
424 best = v;
425 piv = r;
426 }
427 }
428 if (best < pivot_tol)
429 return false;
430 if (piv != col)
431 for (int j = 0; j < N; ++j) {
432 Real t = M[col][j];
433 M[col][j] = M[piv][j];
434 M[piv][j] = t;
435 t = inv[col][j];
436 inv[col][j] = inv[piv][j];
437 inv[piv][j] = t;
438 }
439 const Real invd = Real(1) / M[col][col];
440 for (int j = 0; j < N; ++j) {
441 M[col][j] *= invd;
442 inv[col][j] *= invd;
443 }
444 for (int r = 0; r < N; ++r) {
445 if (r == col)
446 continue;
447 const Real f = M[r][col];
448 if (f != Real(0))
449 for (int j = 0; j < N; ++j) {
450 M[r][j] -= f * M[col][j];
451 inv[r][j] -= f * inv[col][j];
452 }
453 }
454 }
455 return true;
456}
457
459template <int N>
460POPS_HD inline Real mat_norm_inf(const Real (&A)[N][N]) {
461 Real m = Real(0);
462 for (int i = 0; i < N; ++i) {
463 Real r = Real(0);
464 for (int j = 0; j < N; ++j)
465 r += std::fabs(A[i][j]);
466 if (r > m)
467 m = r;
468 }
469 return m;
470}
471
472} // namespace detail
473
489template <int N>
490POPS_HD inline bool roe_abs_apply(const Real (&A)[N][N], const Real (&dU)[N], Real (&out)[N],
491 int max_iter = 80, Real tol = Real(1e-13)) {
492 static_assert(N >= 1 && N <= 16, "roe_abs_apply: 1 <= N <= 16");
493 if (real_spectrum(A) != Spectrum::kReal)
494 return false; // complex/unknown -> caller falls back
495 Real S[N][N], Sinv[N][N];
496 for (int i = 0; i < N; ++i)
497 for (int j = 0; j < N; ++j)
498 S[i][j] = A[i][j];
499 bool converged = false;
500 for (int it = 0; it < max_iter; ++it) {
501 if (!detail::mat_inverse(S, Sinv))
502 return false; // singular iterate (zero eigenvalue)
503 const Real ns = detail::mat_norm_inf(S), nsi = detail::mat_norm_inf(Sinv);
504 Real mu = Real(1);
505 if (ns > Real(0) && nsi > Real(0))
506 mu = std::sqrt(nsi / ns);
507 const Real a = Real(0.5) * mu, b = Real(0.5) / mu;
508 Real diff = Real(0), nrm = Real(0);
509 for (int i = 0; i < N; ++i)
510 for (int j = 0; j < N; ++j) {
511 const Real snext = a * S[i][j] + b * Sinv[i][j];
512 const Real d = snext - S[i][j];
513 diff += d * d;
514 nrm += snext * snext;
515 S[i][j] = snext;
516 }
517 if (nrm > Real(0) && diff <= tol * tol * nrm) {
518 converged = true;
519 break;
520 }
521 }
522 if (!converged)
523 return false;
524 // |A| dU = (A sign(A)) dU = A (S dU) (S = sign(A) commutes with A)
525 Real SdU[N];
526 for (int i = 0; i < N; ++i) {
527 Real s = Real(0);
528 for (int j = 0; j < N; ++j)
529 s += S[i][j] * dU[j];
530 SdU[i] = s;
531 }
532 for (int i = 0; i < N; ++i) {
533 Real s = Real(0);
534 for (int j = 0; j < N; ++j)
535 s += A[i][j] * SdU[j];
536 out[i] = s;
537 }
538 return true;
539}
540
541} // namespace pops
pops::detail::mat_inverse
POPS_HD bool mat_inverse(const Real(&A)[N][N], Real(&inv)[N][N], Real pivot_tol=Real(1e-300))
Inverse of a dense N x N matrix by Gauss-Jordan elimination with partial pivoting,...
Definition dense_eig.hpp:410
pops::detail::hqr_copysign
POPS_HD Real hqr_copysign(Real mag, Real sgn)
Definition dense_eig.hpp:161
pops::detail::mat_norm_inf
POPS_HD Real mat_norm_inf(const Real(&A)[N][N])
Max absolute row sum (infinity norm) of a dense N x N matrix.
Definition dense_eig.hpp:460
pops::detail::hqr_minmax
POPS_HD bool hqr_minmax(Real(&H)[N][N], Real &lmin, Real &lmax, Real &max_im, int max_iter_per_eig)
QR iteration with implicit Francis double shift on a Hessenberg matrix (EISPACK/hqr formulation,...
Definition dense_eig.hpp:183
pops::detail::hessenberg_reduce
POPS_HD void hessenberg_reduce(Real(&H)[N][N])
Upper Hessenberg reduction by Householder reflections, IN PLACE, without accumulating the transformat...
Definition dense_eig.hpp:118
pops::detail::record_eig
POPS_HD void record_eig(Real re, Real im, Real &lmin, Real &lmax, Real &max_im, bool &first)
Accumulate an eigenvalue (re, im) into the current extremes.
Definition dense_eig.hpp:167
pops::detail::gershgorin_bounds
POPS_HD void gershgorin_bounds(const Real(&A)[N][N], Real &lo, Real &hi)
Gershgorin bound on the REAL PARTS: every lambda of the spectrum satisfies lo <= Re(lambda) <= hi (di...
Definition dense_eig.hpp:101
pops
Definition amr_hierarchy.hpp:29
pops::Spectrum
Spectrum
Tri-state classification of a small block's spectrum (ADC-276), returned by pops::real_spectrum.
Definition dense_eig.hpp:93
pops::roe_abs_apply
POPS_HD bool roe_abs_apply(const Real(&A)[N][N], const Real(&dU)[N], Real(&out)[N], int max_iter=80, Real tol=Real(1e-13))
Roe matrix-absolute-value applied to a state jump: out = |A| dU, with |A| the SPECTRAL absolute value...
Definition dense_eig.hpp:490
pops::Real
double Real
Definition types.hpp:30
pops::real_spectrum
POPS_HD Spectrum real_spectrum(const Real(&A)[N][N], Real im_tol=Real(1e-5), int max_iter_per_eig=100)
Classify the spectrum of a small dense block A as kReal / kComplexPair / kUnknown (ADC-276): a GENERI...
Definition dense_eig.hpp:396
pops::real_eig_minmax
POPS_HD EigBounds real_eig_minmax(const Real(&A)[N][N], int max_iter_per_eig=100, bool *fallback=nullptr)
Extremes of the REAL PARTS of the spectrum of a small dense block A, plus the largest |Im| encountere...
Definition dense_eig.hpp:343
pops::EigBounds
Result of real_eig_minmax: real-part extremes + diagnostic.
Definition dense_eig.hpp:53
pops::EigBounds::lmax
Real lmax
largest real part (or Gershgorin upper bound if !converged)
Definition dense_eig.hpp:55
pops::EigBounds::all_real
POPS_HD bool all_real(Real im_tol=Real(1e-5)) const
REAL-SPECTRUM PREDICATE (ADC-276).
Definition dense_eig.hpp:76
pops::EigBounds::lmin
Real lmin
smallest real part (or Gershgorin lower bound if !converged)
Definition dense_eig.hpp:54
pops::EigBounds::converged
bool converged
false -> Gershgorin fallback (valid external bound, NOT the spectrum)
Definition dense_eig.hpp:59
pops::EigBounds::has_complex_pair
POPS_HD bool has_complex_pair(Real im_tol=Real(1e-5)) const
Complement of all_real RESTRICTED to converged blocks: true iff the spectrum was computed AND carries...
Definition dense_eig.hpp:84
pops::EigBounds::max_im
Real max_im
largest |Im(lambda)| encountered (0 = real spectrum: hyperbolic).
Definition dense_eig.hpp:56
types.hpp
Base scalar types and the POPS_HD macro (host+device portability).
POPS_HD
#define POPS_HD
Definition types.hpp:25
Generated by doxygen 1.9.8