Reusable Householder QR solve for a shifted matrix with residual status

Input: a stored real square matrix $M\in {\mathbb{R}}^{d\times d}$ with $d\ge 1$; a unit roundoff $u>0$ with $du<1$; and a terminating certificate evaluator evaluator. For each requested Euclidean or Frobenius norm of stored data, evaluator returns a finite stored nonnegative upper bound for the exact norm or unavailable. For each requested nonnegative expression formed from such upper bounds, $u$, and $d$, it returns a finite stored nonnegative upper bound for the exact expression or unavailable. Every returned bound is sound. All reflector, substitution, dot-product, and residual operations below use Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add; certificate-evaluator calls return certified bounds rather than relying on ordinary rounded evaluation of a bound expression.

Denote the top-level factor-only operation by factor(M,u,evaluator). It returns either a reusable factor object or a factorization failure status.

1. Set $R:=M$. For each column $j=0,\dots ,d-1$, let $x$ be the active part of column $j$ of $R$ and $B$ the active trailing block to its right. Call Scaled Householder reflector formation and application on $(x,B)$. If it returns range-failure, return that status and no factor object. If it returns a marked identity or reflector result, replace the active block by its returned $C_{+}$, store its marker and (when present) $(v,\widehat{\tau })$, and continue. For any other return, stop and return factorization-failure with no factor object. The resulting $R$ is upper triangular with its lower entries set to exact zeros by the called procedure; this is the full/reduced computational QR convention of Full, reduced, and column-pivoted computational QR factorisations.
2. If any diagonal entry $R_{jj}$ is zero, return singular-factor. Otherwise return a factor object containing the original stored $M$, $R$, the ordered stored reflectors, $u$, and the certificate evaluator. The factor object's operation solve(b,epsilon), whose input is any stored real $b\in {\mathbb{R}}^d$ and finite $ϵ\ge 0$, is:
   1. Apply the stored reflector updates in the same order, each using the successful rounded dot, scaling, product, and subtraction operations specified by Scaled Householder reflector formation and application, to obtain $c$. Return the first range-failure if an operation fails.
   2. Compute $x$ by descending rounded back substitution for $Rx=c$. For $i=d-1,d-2,\dots ,0$, initialize $s:=0$; for $j=i+1,i+2,\dots ,d-1$ in increasing order, compute $p:=\mathrm{fl}(R_{ij}{x}_j)$ and then $s:=\mathrm{fl}(s+p)$; next compute $t:=\mathrm{fl}(c_i-s)$ and $x_i:=\mathrm{fl}(t\mathrm{/}{R}_{ii})$. If any product, addition, subtraction, or division fails, return the first range-failure. This is the descending recurrence of Forward and backward substitution for triangular linear systems with its rounded operation order made explicit.
   3. Compute ${\widehat{v}}_i$ for $i=0,\dots ,d-1$ as the rounded left-to-right dot product of row $i$ of the original $M$ with $x$: initialize the accumulator with $\mathrm{fl}(M_{i0}{x}_0)$ and for $j=1,\dots ,d-1$ compute $p:=\mathrm{fl}(M_{ij}{x}_j)$ and then replace the accumulator $s$ by $\mathrm{fl}(s+p)$. Set the resulting accumulator to ${\widehat{v}}_i$. Then compute each ${\widehat{r}}_i:=\mathrm{fl}({\widehat{v}}_i-b_i)$. Return the first range-failure if any operation fails.
   4. Request evaluator bounds $\bar{m}\ge \mathrm{\parallel }M{\mathrm{\parallel }}_F$, $\bar{x}\ge \mathrm{\parallel }x{\mathrm{\parallel }}_2$, $\bar{b}\ge \mathrm{\parallel }b{\mathrm{\parallel }}_2$, $\bar{v}\ge \mathrm{\parallel }\widehat{v}{\mathrm{\parallel }}_2$, and $\bar{r}\ge \mathrm{\parallel }\widehat{r}{\mathrm{\parallel }}_2$. If any request is unavailable, return (residual-unverified, x, null). Request a certified upper bound $\bar{W}$ for ${\gamma }_d\bar{m}\bar{x}+u(\bar{v}+\bar{b})$, where ${\gamma }_d:=du\mathrm{/}(1-du)$. If this request is unavailable, return (residual-unverified, x, null). Request a certified upper bound $\rho$ for $\bar{r}+\bar{W}$; if unavailable, return (residual-unverified, x, null). By A computed shifted-system residual has a certified upper bound, $\rho \ge \mathrm{\parallel }Mx-b{\mathrm{\parallel }}_2$. Return (success, x, rho) if $\rho \le ϵ$; otherwise return (residual-not-certified, x, rho).

This is a factor-only operation followed by an independently callable solve operation: the factor object accepts each later stored right-hand side and tolerance without receiving a precomputed list. Each solve performs one triangular solve and one residual check using the same stored factorization. A returned residual bound is an equation-residual certificate, not a forward-error estimate; converting it to forward error requires a separate inverse bound. No rank-revealing guarantee, sparse cost bound, or success guarantee for an ill-conditioned matrix is part of this procedure.

Factorization terminates with singular-factor, range-failure, factorization-failure, or a reusable factor object. Each call to the object's solve operation terminates with success, range-failure, residual-unverified, or residual-not-certified, together with any completed approximate solution and residual bound.