Finite-budget computed Francis Schur reduction with explicit statuses

Input a real or complex square matrix $A\in {\mathbb{F}}^{n\times n}$ with $n\ge 1$, a finite nonnegative deflation tolerance $\tau$, a nonnegative integer sweep budget $N$, a positive integer stagnation period $h$, a finite family of exceptional-shift tables $T_{a,b}$ indexed by original contiguous index intervals $[a,b]\subseteq [1,n]$ with $a<b$, optional exact rational/dyadic encodings of its stored real values, optionally a terminating exact evaluator of Exact sign predicate for a real 2x2 characteristic discriminant, and optionally a terminating outward evaluator for nonnegative expressions and Frobenius norms. The outward evaluator's upadd, upmul, and upsqrt calls return finite nonnegative represented upper bounds for the exact sum, product, and square root, while upnormF(X) returns a finite represented upper bound for $\mathrm{\parallel }X{\mathrm{\parallel }}_F$; any call may instead return unavailable. In the declared scalar-cost model, each upadd, upmul, and upsqrt call has unit cost and upnormF(X) makes one pass over the entries of $X$ with a fixed number of such calls per real component. The evaluator acts only on certificate bookkeeping. A formation breakdown is propagated with the operation index and partial factors. Every scalar operation in the main procedure is evaluated in the operationwise model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. Any invalid, overflow, or underflow operation in the main procedure returns range-failure with the operation index and partial factors.

Each complex table entry is a finite represented scalar $\mu \in \mathbb{C}$; each real table entry is a pair of finite represented real coefficients $(s,t)$ for $z^2-sz+t$. Each table is an ordered finite list, possibly empty; an omitted interval has the empty list. These descriptors use the field of $A$. Tables support direct indexed lookup by interval and cursor; unused entries are not scanned. All interval-state updates use a single ordered pass over the current partition. Counters and table indices are exact integers. Maintain state $(c_{a,b},q_{a,b})$, initially $(0,0)$, for each current maximal diagonal interval separated by exact zero subdiagonals: $c$ counts completed ordinary block sweeps since the last reset and $q$ counts exceptional descriptors already consumed. After a split, discard the parent state and initialize every new child interval to $(0,0)$, using its own table. An unchanged interval retains its state. Completed one-by-one blocks, and real two-by-two blocks with verified negative discriminant, require no shift.

Start with $H=A$ and the empty factor product $Q=I$. Set the completed global sweep count to $k=0$ after initialization. Before attempting any sweep, and after every completed sweep, perform steps 2 and 5, in that order. If neither returns, snapshot the remaining active intervals in increasing left endpoint, and perform steps 3 and 4 once for each snapshotted interval in that order. No deflation scan or repartition occurs during this pass. After all its intervals finish, increment $k$ by one and repeat the scan and stopping test. All intervals in a pass use their current diagonal block when reached; previously applied similarities on disjoint blocks do not reschedule a block.

Before every real or complex formation in steps 1 and 4, test the target vector componentwise for exact zero, including an external initial polynomial vector. If it is zero, mark the stage identity, copy all affected matrix entries unchanged, skip formation and both left and right updates, append the embedded identity $P_j=I$ to the product representation, and record $E_j=0$ and arithmetic coefficient ${\alpha }_j={\bar{\alpha }}_j=0$. Do not evaluate $P(v)$ or a stored-vector norm in this branch. For a matrix target its column is already zero, so its target zero writes change nothing; for an external vector no matrix column is overwritten. Continue with the invariant checks, sweep counters and stopping tests specified below. Record and process every deflation write separately exactly as specified below; this branch does not discard a deflation perturbation. These identity instructions take precedence over the formation/update instructions in steps 1 and 4.

For each nonidentity stored reflector with active indices $S=(s_1,\dots ,s_m)$, nonzero vector $v\in {\mathbb{F}}^m$, and real computed coefficient $\widehat{\tau }$, the real or complex formation routine first performs the stated finite left application. Then perform the right application in this fixed order: for each matrix row $i=1,\dots ,n$, compute $d_i=\mathrm{fl}({\sum }_{a=1}^m{H}_{i,s_a}^L{v}_a)$ by left-to-right rounded products and additions; for $a=1,\dots ,m$, compute $c_a=\mathrm{fl}(\widehat{\tau }\overline{v_a})$, $g_{ia}=\mathrm{fl}(d_i{c}_a)$, and $H_{i,s_a}^R=\mathrm{fl}(H_{i,s_a}^L-g_{ia})$. Each complex product uses the six real operations of the arithmetic model. The exact factor appended to $Q$ is $P(v)=I-2v{v}^{\ast }\mathrm{/}(v^{\ast }v)$; the computed right update uses $\widehat{\tau }$ as displayed and is not asserted equal to multiplication by $P(v)$. For a matrix-target column, the formation routine's explicit zero writes are retained. In the initial polynomial-vector step there is no matrix-target column: use the entire active row panel for the left update and then apply the right update to the active columns. Record deflation writes separately, including a pure deflation with no following reflector as an identity-transform stage; a local error certificate must account for all other rounding and target-column writes.

1. Reduce $A$ to upper Hessenberg form with Householder similarities. For $k=1,\dots ,\max (n-2,0)$ in increasing order, use active indices $S=\{k+1,\dots ,n\}$, target $x=H_{S,k}$, and trailing left panel $B=H_{S,\{k+1,\dots ,n\}}$ with columns in increasing order. In the real branch use Scaled Householder reflector formation and application for the finite left application and overwrite; in the complex branch use Scaled complex Householder reflector step with explicit range status. Complete each similarity by the right application above on columns $S$. After the entire initialization loop, including when it is empty, check the resulting matrix's exact upper-Hessenberg zero pattern; on failure return bulge-chase-invariant-failure with the partial factors. Store the exact unitary reflector factors as a product representation $Q_0$.
2. In increasing global row index $i=2,\dots ,n$, test each nonzero first-subdiagonal entry by $\mathrm{\mid }{h}_{i,i-1}\mathrm{\mid }\le \tau (\mathrm{\mid }{h}_{i-1,i-1}\mathrm{\mid }+\mathrm{\mid }{h}_{i,i}\mathrm{\mid })\mathrm{.}$ For a complex entry $z=a+ib$, evaluate $\mathrm{\mid }z\mathrm{\mid }$ by $s=\max (\mathrm{\mid }a\mathrm{\mid },\mathrm{\mid }b\mathrm{\mid })$ and $s\sqrt{(a\mathrm{/}s{)}^2+(b\mathrm{/}s{)}^2}$ when $z\mathrm{\ne }0$, with divisions, squares, sum, square root and final scaling rounded in the stated order; use $\mathrm{\mid }0\mathrm{\mid }=0$. Return range-failure if this evaluation fails. Compute the diagonal-magnitude sum first, then multiply by $\tau$ in the operationwise model; compare the resulting finite values exactly. When the test holds, set that entry to zero and record its original value as a deflation perturbation. After the scan, form the maximal intervals separated by exact zero subdiagonals and update interval states by the rule above. The equality case deflates. Do not use the Hermitian-tridiagonal Wilkinson shift rule for a general Hessenberg matrix.
3. For the current snapshotted interval $[a,b]$, let $L=H_{a:b,a:b}$ and $\mathrm{\ell }=b-a+1$. If $c_{a,b}=h$, read entry $q_{a,b}+1$ of $T_{a,b}$. If absent, return stagnation-exhausted before this block's chase. Otherwise consume it by incrementing $q_{a,b}$ once and use its complex scalar shift or real coefficient pair. After this exceptional block sweep completes, set $c_{a,b}=0$. If $c_{a,b}<h$, use the ordinary shift: in the complex branch $\mu =L_{\mathrm{\ell },\mathrm{\ell }}$; in the real branch, for the trailing block $\left(\begin{array}{cc}{a}^{\mathrm{\prime }}& b^{\mathrm{\prime }}\\ c^{\mathrm{\prime }}& d^{\mathrm{\prime }}\end{array}\right)$, compute $s=\mathrm{fl}(a^{\mathrm{\prime }}+d^{\mathrm{\prime }})$ then $t=\mathrm{fl}(\mathrm{fl}(a^{\mathrm{\prime }}{d}^{\mathrm{\prime }})-\mathrm{fl}(b^{\mathrm{\prime }}{c}^{\mathrm{\prime }}))$ in that order. Return range-failure on a failed operation. After this ordinary block sweep completes, increment $c_{a,b}$ by one. A split at the next scan resets each new child as specified above; unchanged intervals retain these updated states.
4. On the diagonal block $L$ selected in step 3, apply the index pattern of Exact finite Francis single- or double-shift Hessenberg QR step, with its dimension $\mathrm{\ell }$ and local indices $1,\dots ,\mathrm{\ell }$. In the vector formulas in this step, $H$ means that current diagonal block. Translate every local index $i$ to global index $a+i-1$; embed each factor in the full $n\times n$ identity. The left panel extends through global column $n$, and the right application runs through all global rows, so upper off-diagonal blocks are updated as well. For the initial single-shift vector compute $x_1=\mathrm{fl}(h_{11}-\mu )$ and $x_i=h_{i1}$ for $i>1$. For the initial real double-shift vector, compute $u=H{e}_1$, then $w=Hu$ by left-to-right rounded matrix-vector products, then $x=\mathrm{fl}(\mathrm{fl}(w-su)+t{e}_1)$ componentwise. Compute only the first $r+1$ components from the Hessenberg support pattern; on range failure return the operation index and partial factors. For a chase at local column $j$, use the globally translated active indices $S=I_j$, the corresponding global target column $a+j-1$, and the trailing left panel through global columns $a+j,\dots ,n$ in increasing order. Call the real or complex formation routine, write the target tail zeros exactly, complete the two-sided update by the right-kernel order above on columns $S$, and append the exact reflector representation to $Q$. In the initial polynomial-vector step there is no matrix-target column; use the full active row panel for the left update and then apply the right update to its active columns. After the block chase, perform the following exact upper-Hessenberg zero-pattern check, including when all its stages were identities. For $i=1,\dots ,n$ in increasing order, if $a\le i\le b$, test $h_{ij}=0$ for $j=1,\dots ,i-2$ in increasing order; otherwise test $h_{ij}=0$ for $j=a,\dots ,\min (b,i-2)$ in increasing order. Skip every empty index range. On the first nonzero tested entry return bulge-chase-invariant-failure with the partial matrix and factors. These are precisely the below-first-subdiagonal positions whose row or column belongs to $[a,b]$. Enumerate only these ranges; do not scan all $n^2$ positions and filter them. The unchanged complementary positions retain the zero values established by initialization or the previous accepted check. Deflation writes affect only the first subdiagonal and do not change those positions. Use Left and right scaled Householder errors give a concrete similarity bound to store the proved two-sided coefficient ${\alpha }_{m_j}=a_{m_j}+b_{m_j}+a_{m_j}{b}_{m_j}$, with $a$ the target-overwrite coefficient or separate-block coefficient as actually used and $b$ the right separate-block coefficient. This gives $\mathrm{\parallel }{E}_j{\mathrm{\parallel }}_F\le {\alpha }_{m_j}\mathrm{\parallel }{\widehat{H}}_{j-1}{\mathrm{\parallel }}_F$ when its $m_{j}u<1$ and range hypotheses hold. Apply the same certificate to every Hessenberg-initialization similarity. If any step lacks that certificate, retain the residual diagnostic and set the backward-error-certificate status to inconclusive.
5. Inspect maximal intervals in increasing left endpoint. A size-one interval is completed. A complex interval of size at least two, or a real interval of size at least three, is active. For each real size-two interval determine its exact discriminant sign as below: a negative sign completes that interval, while zero or a positive sign leaves it active. Stop with Schur-computed when all intervals are completed. Use the supplied evaluator, or the finite integer recipe in Rational or dyadic stored real blocks have a finite exact discriminant-sign evaluator when exact rational/dyadic encodings are available, for each real block's $\mathrm{discSign}$ as defined in Exact sign predicate for a real 2x2 characteristic discriminant. If no exact sign is available, return schur-type-undecided and do not claim Schur-computed. If any active interval remains and $k=N$, return budget-exhausted; in particular, with $N=0$ return Schur-computed only when the initialized matrix already has the target block structure. Return the partial matrix and factors on every failure or nonconvergence status.

For certificate bookkeeping, start before Hessenberg initialization with ${\rho }_0=\mathrm{upnormF}(A)$ and $C_0=0$. If the evaluator is absent or this call is unavailable, mark the backward-error certificate inconclusive. Process each local similarity update and each deflation-only write as a separate finite stage $j$. For a local similarity, take the actual coefficient ${\alpha }_j=a_{m_j}+b_{m_j}+a_{m_j}{b}_{m_j}$ from Left and right scaled Householder errors give a concrete similarity bound and compute ${\bar{\alpha }}_j=\mathrm{upadd}(\mathrm{upadd}(a_{m_j},b_{m_j}),\mathrm{upmul}(a_{m_j},b_{m_j}))\ge {\alpha }_j$; for an identity update put ${\bar{\alpha }}_j=0$. A recorded deflation entry stores its old value, while its additive matrix $D_j$ has entry equal to zero minus that old value at the written position. For a deflation-only stage use $P_j=I$, $E_j=0$, ${\bar{\alpha }}_j=0$, and obtain ${\bar{\beta }}_j\ge \mathrm{\parallel }{D}_j{\mathrm{\parallel }}_F$ by calling upnormF on the one-entry list containing the signed recorded value, whose Frobenius norm equals $\mathrm{\parallel }{D}_j{\mathrm{\parallel }}_F$; do not form or scan a dense matrix $D_j$; for a similarity stage set ${\bar{\beta }}_j=0$. With the outward evaluator compute

$g_j=\mathrm{upadd}(\mathrm{upmul}({\bar{\alpha }}_j,{\rho }_{j-1}),{\bar{\beta }}_j),{\rho }_j=\mathrm{upadd}({\rho }_{j-1},g_j),C_j=\mathrm{upadd}(C_{j-1},g_j)\mathrm{.}$

If any local coefficient certificate or evaluator call is unavailable, keep the residual diagnostic and mark the backward-error-certificate status inconclusive; this ancillary failure does not change the main algorithm status. Otherwise return $C_J$ as the accumulated backward-error upper bound and mark that status certified. The number of recurrence calls is finite: one initialization norm, a fixed number of scalar upper-bound calls per similarity stage, and one norm call per deflation-only stage.

Let $\widehat{T}$ be the final returned matrix. The output includes the product representation $Q$, residual matrix $R=Q^{\ast }AQ-\widehat{T}$ and its Frobenius norm as diagnostics, the exact orthogonality bound $0$ for that reflector product, the recorded deflation entries, and either the accumulated backward-error bound when all local certificates exist or certificate status inconclusive. Bound expressions and any claimed upper bounds are evaluated exactly or with rigorous outward allowances; uncorrected rounded residual norms are diagnostics only. A rounded dense matrix formed from $Q$ has its own orthogonality residual and is not asserted exactly unitary. The residual diagnostic does not itself give a forward eigenvalue bound. The procedure asserts no unconditional convergence of any shift schedule.