Spectrum extremes (real parts) of a small dense matrix: signed wave-speed bounds supplied by a model (exact HLL via flux Jacobian). More...

#include <cmath>
#include <limits>
#include <pops/core/foundation/types.hpp>

Include dependency graph for dense_eig.hpp:

Go to the source code of this file.

Classes
struct	pops::EigBounds
	Result of real_eig_minmax: real-part extremes + diagnostic. More...

Namespaces
namespace	pops

namespace	pops::detail

Enumerations
enum class	pops::Spectrum : int { pops::kReal = 0 , pops::kComplexPair = 1 , pops::kUnknown = 2 }
	Tri-state classification of a small block's spectrum (ADC-276), returned by pops::real_spectrum. More...

Functions
template<int N>
POPS_HD void	pops::detail::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 (disks centered at a_ii of radius the sum of the \|off-diagonal\| terms of the row).

template<int N>
POPS_HD void	pops::detail::hessenberg_reduce (Real(&H)[N][N])
	Upper Hessenberg reduction by Householder reflections, IN PLACE, without accumulating the transformations (eigenvalues only).

POPS_HD Real	pops::detail::hqr_copysign (Real mag, Real sgn)

POPS_HD void	pops::detail::record_eig (Real re, Real im, Real &lmin, Real &lmax, Real &max_im, bool &first)
	Accumulate an eigenvalue (re, im) into the current extremes.

template<int N>
POPS_HD bool	pops::detail::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, eigenvalues only, blocks processed bottom-up with deflation).

template<int N>
POPS_HD EigBounds	pops::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\| encountered and a convergence indicator (see the file header for the full contract: Gershgorin fallback on non-convergence, max_im as a hyperbolicity-loss detector).

template<int N>
POPS_HD Spectrum	pops::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 GENERIC, device-safe predicate over the SAME Francis-QR path as real_eig_minmax (no second algorithm to keep in sync).

template<int N>
POPS_HD bool	pops::detail::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, into `inv`.

template<int N>
POPS_HD Real	pops::detail::mat_norm_inf (const Real(&A)[N][N])
	Max absolute row sum (infinity norm) of a dense N x N matrix.

template<int N>
POPS_HD bool	pops::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 A * sign(A).

Detailed Description

Spectrum extremes (real parts) of a small dense matrix: signed wave-speed bounds supplied by a model (exact HLL via flux Jacobian).

GENERIC utility: the core knows nothing of the calling model – it receives a dense block Real[N][N] and returns min/max of the real parts of its spectrum. Intended consumer: the DSL codegen wave_speeds_from_jacobian (the model supplies dF/dU, possibly by diagonal blocks; each block passes through here). Any other "small eigenvalues on the stack" use is legitimate.

Algorithm (eigenvalues ONLY, never the vectors): N == 1, 2: closed form (trace/determinant); N >= 3: Hessenberg reduction (Householder, without accumulation) then QR iteration with implicit Francis DOUBLE SHIFT and deflation (EISPACK/hqr formulation) – complex conjugate pairs stay in real arithmetic (2x2 blocks), exceptional shifts after 10 and 20 iterations on a single block.

ROBUSTNESS CONTRACT: if a block does not converge under the iteration cap, the result is the Gershgorin FALLBACK over the WHOLE matrix (converged = false): an EXTERNAL bound always valid for all real parts (sL <= all speeds <= sR, hence a stable HLL flux, just more diffusive). If the fallback triggers, lmin/lmax are NOT the eigenvalues and max_im is 0 by CONVENTION (nothing was computed, this is NOT a real-spectrum signal): any bit-match against an eig reference is then void. converged (or the fallback output parameter) decides: never read lmin/lmax/max_im without consulting it. The default cap (100) is sized so that the usual near-degenerate companion blocks converge; the fallback remains the safety net for out-of-envelope cases.

HYPERBOLICITY: max_im returns the largest |Im(lambda)| encountered. A hyperbolic system has a real spectrum (max_im ~ 0); a model that loses hyperbolicity does not silently receive a plausible-but-wrong speed – the caller decides (assert, warning, clamp). A caller wanting only a real/complex verdict (not the extremes) uses pops::real_spectrum / EigBounds::all_real (ADC-276), which couple this max_im check with convergence so a non-converged block is never reported real.

PRECISION: simple, separated eigenvalues -> machine precision (tested at rtol 1e-10). CLUSTERED eigenvalues of a non-symmetric matrix: conditioning ~ eps^(1/m) for a near-multiplicity m (a limit of the problem, not of the algorithm) – do not expect 1e-10 at a triple point. No exactness claim is made under the iteration cap: the cap bounds the COST, the fallback bounds the RESULT.

Device: POPS_HD, stack buffers (O(N^2)), zero allocation, zero std:: container, bounded loops, no recursion – only std::sqrt / std::fabs (resolved on device under nvcc/Kokkos, as on the flux path).

include/pops/numerics/linalg/dense_eig.hpp File Reference¶

Classes

Namespaces

Enumerations

Functions

Detailed Description