Scalar multifrontal Cholesky and solves on the elimination forest

Input a Hermitian positive-definite matrix $A\in {\mathbb{C}}^{n\times n}$ with $n\ge 1$, its Hermitian structural graph from Undirected structural graph of a Hermitian matrix, and an ordering $\pi =(v_0,\dots ,v_{n-1})$. Let $P$ be the permutation matrix whose row $r$ is row $v_r$ of $I_n$, and put $B=PA{P}^{\ast }$. Use the filled graph and parent forest of Structural graphs and ordered symmetric elimination and the column patterns of Elimination-tree column patterns and reach. For each original vertex $v_r$, let its front index set be $F_r=\{v_r\}\cup N_r^{+},$ where $N_r^{+}$ is its later-neighbor set in the filled graph. Assign each nonzero off-diagonal matrix entry $A_{pq}$ and each diagonal entry $A_{pp}$ once to the front of the endpoint eliminated first (the diagonal goes to its own front); zero off-diagonal entries need no assignment. Treat front arrays as indexed by original-vertex labels, but store the factor in ordered coordinates of $B$: for $v_q\in F_r\setminus \{v_r\}$, its factor entry is $L_{qr}$, where $q>r$ is the ordering position. Initialize $L$ to zero on ordered diagonal and lower-triangular positions, allocating every position permitted by the filled column pattern, including structural zeros. Use the operation model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add for all finite arithmetic, including contribution assembly and solves. A range failure during assembly, factorization or solving returns range-failure with the current position.

Traverse fronts in increasing pivot position $r$ (every child has smaller position than its parent). At front $r$, initialize a zero dense Hermitian array ${\mathrm{\Phi }}_r$ indexed by $F_r$ and add the entries assigned to $v_r$. Process the children in increasing pivot position $s$; for each child, add its stored contribution array to the entries whose indices lie in $F_s\setminus \{v_s\}$, using the common original-vertex labels. The inclusion needed for this addition is supplied by Every later filled neighbor lies on the elimination-tree parent path. For the pivot $v_r$, let $d=({\mathrm{\Phi }}_r{)}_{v_r{v}_r}$. If it is not finite and positive, return invalid-pivot$(r,d)$; otherwise compute $l_{rr}=\sqrt{d}$ and, for each later position $q$ with $v_q\in F_r$, compute $l_{qr}=({\mathrm{\Phi }}_r{)}_{v_q{v}_r}\mathrm{/}{l}_{rr}$ in increasing $q$. Apply Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add to each operation and return range-failure$(r)$ at the first range error. Store the computed diagonal and all entries $L_{qr}$ in the preallocated slots. For lower-triangular trailing positions with later positions $q\ge p>r$ and $v_q,v_p\in F_r$, in increasing $(q,p)$ order, update $({\mathrm{\Phi }}_r{)}_{v_q{v}_p}\leftarrow \mathrm{fl}(({\mathrm{\Phi }}_r{)}_{v_q{v}_p}-l_{qr}\overline{l_{pr}})$; set the Hermitian counterpart by conjugation and keep the diagonal real. Send the trailing contribution array to the parent; at a root discard only the empty final remainder. A failed pivot or range operation stops the factorization with its indexed status. These scalar steps have the same exact pivot and Schur update as Ordered sparse Cholesky with explicit pivot and range statuses.

After a successful factorization, input a right-hand side $b\in {\mathbb{C}}^n$ and use these solve instructions:

1. Set $x=Pb$.
2. For $j=0,\dots ,n-1$, divide $x_j$ by $l_{jj}$, then for every stored $l_{ij}$ with $i>j$ in increasing $i$ set $x_i\leftarrow x_i-l_{ij}{x}_j$.
3. For $j=n-1,\dots ,0$, set $x_j\leftarrow (x_j-{\sum }_{i>j}\overline{l_{ij}}{x}_i)\mathrm{/}\overline{l_{jj}}$, accumulating the sum from left to right in increasing $i$, using the backward-substitution convention of Forward and backward substitution for triangular linear systems.
4. Return $z=P^{\ast }x$.

The factorization returns either $(P,L)$ or a pivot/range status. The solve returns the computed vector $z$ or range-failure with the current position. The solve arithmetic uses the declared finite-precision model; no exact-solution or finite-precision error guarantee is asserted here.