Skip to content

num::mixed_lu

Defined in header <linear/factorization/mixed_lu.hpp>

struct mixed_lu;

\(A\) factored in single precision, solved to double-precision accuracy.

The \(O(n^3)\) factorization runs in float, where a SIMD register holds twice as many values, so it moves half the memory and does about twice the arithmetic per instruction. Each solve then refines in double: with \(x_0 = \hat{A}^{-1}b\) from the float factors, repeat \(r = b - Ax_k\) in double (\(O(n^2)\)) and \(x_{k+1} = x_k + \hat{A}^{-1}r\). Each step shrinks the error by about \(\kappa(A)\,u_{\text{float}}\), so it reaches double accuracy in a few steps whenever \(\kappa(A)\) is well below \(1/u_{\text{float}} \approx 10^7\).

Construction checks that with inverse_norm1_estimate on the float factors. When \(A\) is too ill-conditioned it also factors in double, and solves use that instead: refined() says which.

Members

member description
static constexpr real max_contraction = 0.1 \(\kappa(A)u_{\text{float}}\) above which refinement is not trusted to converge.
static constexpr int max_steps = 30 Refinement steps before a solve settles for its current iterate, as LAPACK dsgesv.
mat<real> A
lu_result<float> low
std::optional<lu_result<real>> fallback
real condition_estimate = 0.0 \(\Vert A\Vert _1\) times the estimate of \(\Vert A^{-1}\Vert _1\).
real norm_1 = 0.0 \(\Vert A\Vert _1 = \Vert A^T\Vert _\infty\), for the transposed stopping test.
real norm_inf = 0.0 \(\Vert A\Vert _\infty\), for the stopping test.
idx size() const
bool refined() const True when solves refine the float factors; false when they use fallback.