Threshold partial-pivot sparse LU with explicit fill and failure statuses

Input a nonsingular real matrix $A\in {\mathbb{R}}^{n\times n}$ with $n\ge 1$, an initial simultaneous row-and-column permutation matrix $P$ under the row-action convention of Permutation matrices, partial pivoting, and the pivot-growth factor, and a threshold $\theta \in (0,1]$. Choose the exact-arithmetic branch or the finite-precision branch. In the latter, additionally supply a format from Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add; the entries of $A$ must already be stored values in that format, with no input conversion. The exact-arithmetic branch permits arbitrary real entries and evaluates every scalar operation exactly, without finite-format range restrictions. Initialize the active matrix $S^{(0)}=PA{P}^T$ by permuted copies, the accumulated row permutation $Q=I_n$, the unit lower-triangular factor $L=I_n$, and the upper-triangular factor $U=0$. Treat an unstored active entry as zero. In finite precision every real scalar operation uses the supplied model; if an exact intermediate is nonzero and subnormal, overflows, or the operation reports range failure, return range-failure with its stage. Pivot tests and tie-breaking comparisons are exact guards on current values. For $k=0,\dots ,n-1$, ties for a maximum are resolved by choosing the smallest current row index.

At stage $k$, let $m_k={\max }_{i\ge k}\mathrm{\mid }{S}_{ik}^{(k)}\mathrm{\mid }$. If $m_k=0$, return singular-pivot$(k)$. If $\mathrm{\mid }{S}_{kk}^{(k)}\mathrm{\mid }\ge \theta m_k$, select row $p=k$; otherwise select the smallest row attaining $m_k$. If $p\mathrm{\ne }k$, swap rows $p$ and $k$ of the active matrix and the already stored entries in columns $0,\dots ,k-1$ of those rows of $L$. Also replace $Q$ by $E_{kp}Q$, where $E_{kp}$ is the permutation matrix swapping rows $k$ and $p$. Do not change the previously stored rows of $U$. If the selected pivot is zero, return singular-pivot$(k)$.

Save the active pivot row by setting $U_{kj}=S_{kj}^{(k)}$ for every $j\ge k$; entries $U_{kj}$ with $j<k$ remain zero. Set $L_{ik}=S_{ik}^{(k)}\mathrm{/}{S}_{kk}^{(k)}$ for every $i>k$, interpreting a missing $S_{ik}^{(k)}$ as zero. For every $i,j>k$, update $S_{ij}^{(k+1)}=S_{ij}^{(k)}-L_{ik}{U}_{kj}\mathrm{.}$ In a sparse representation it suffices to visit positions where $S_{ij}^{(k)}\mathrm{\ne }0$ or $L_{ik}{U}_{kj}\mathrm{\ne }0$; missing entries in this formula are zero. Allocate any position whose rounded update is nonzero; this is numerical fill. If an arithmetic operation is not representable or reports a range error, return range-failure with its stage. Any failure return is a status only and asserts no factorization. On success return $(Q,P,L,U)$. In exact arithmetic the updates give $QPA{P}^T=LU$; the procedure makes no claim of a universal fill or growth bound.