Preconditioner algorithms and residual metrics

Let $F$ be a field, let $n\ge 1$, and let the input be a valid $n\times n$ CSR record as in Compressed sparse row matrices with stored-entry semantics, representing $A\in M_n(F)$. Matrix positions in this definition use one-based notation: $(i,j)$ refers to CSR row $i-1$ and column index $j-1$. This row/column order is fixed before the procedure begins; vector component label $i$ refers to coordinate $i-1$. The procedure performs no pivoting or reordering. All field operations below are exact.

Let $P=\{(i,j):a_{ij}\mathrm{\ne }0\}$. Scan each row once to check the canonical and increasing-index conditions and to record whether its diagonal position occurs. If the input fails either pattern condition, return invalid-pattern. If a row has no diagonal position, return missing-diagonal(i) for the least such $i$.

Copy the value in each input slot to a work array $W$, and regard the unit diagonal of $L$ as implicit. Process rows $i=1,\dots ,n$ in increasing order. For each $k<i$ with $(i,k)\in P$, in increasing $k$ order, do the following:

1. If $W_{kk}=0$, return zero-pivot(k).
2. Set $W_{ik}\leftarrow W_{ik}\mathrm{/}{W}_{kk}$; this value is the stored factor $L_{ik}$.
3. Let $R$ be the increasing list of one-based matrix column labels $j>k$ represented by row $i$'s CSR segment, and let $C$ be the corresponding increasing list represented by row $k$'s CSR segment. Set each pointer to the first entry of its list, or to that list's end marker when it is empty. While both pointers are within their lists, compare their current indices. If they are equal, call their common value $j$ and update $W_{ij}\leftarrow W_{ij}-W_{ik}{W}_{kj}$ and advance both pointers once. If the current index in $R$ is smaller, advance only its pointer once; if the current index in $C$ is smaller, advance only its pointer once, making no arithmetic update in either unequal-index branch. Stop when either pointer passes the end of its list. Do not skip a matching slot when a current work value is zero.

After all such $k$ have been processed in row $i$, if $W_{ii}=0$, return zero-pivot(i). Otherwise the upper row entries $U_{ij}$ are the current $W_{ij}$ for $(i,j)\in P$ and $j\ge i$. Continue to the next row.

If every row completes, return ok with the unit-lower factor $L$ (stored strict-lower entries $L_{ik}=W_{ik}$ for $(i,k)\in P$) and the upper factor $U$ (stored entries $U_{ij}=W_{ij}$ for $(i,j)\in P$, $j\ge i$). Both factors use only positions in $P$; the diagonal of $L$ is implicit, and all other unstored off-diagonal entries are zero. The numeric values remain stored at their structural slots if an update makes a value zero.

For a successful factor pair and $r\in F^n$, apply the preconditioner as follows. For $i=1,\dots ,n$ in order, initialize $y_i:=r_i$ and, for each $k<i$ with $(i,k)\in P$ in increasing order, update $y_i:=y_i-L_{ik}{y}_k$. Then for $i=n,n-1,\dots ,1$, initialize $x_i:=y_i$; for each $k>i$ with $(i,k)\in P$ in decreasing order, update $x_i:=x_i-U_{ik}{x}_k$; finally set $x_i:=x_i\mathrm{/}{U}_{ii}$. Visit every such stored factor slot exactly once, including slots whose current factor value is zero. These are the forward and backward substitution recurrences, with their sums traversed only over stored strict-triangular slots. The pivot checks ensure every denominator in this application is nonzero. This procedure defines setup and application only; it asserts no factorization identity $A=LU$, Krylov convergence, conditioning improvement, or accuracy guarantee.

Let $F$ be a field, $n\ge 1$, and $A\in M_n(F)$ have a valid canonical CSR record with strictly increasing column indices in each row. Let $P=\{(i,j):a_{ij}\mathrm{\ne }0\}$ and $s=\mathrm{nnz}(A)$. Run Fixed-pattern ILU(0) factorization and application. The procedure terminates with a missing-diagonal or zero-pivot status, or returns ok.

On an ok return, define $P_L=\{(i,k)\in P:i>k\},P_U=\{(k,j)\in P:k\le j\},$ and put $\mathrm{\ell }=\mathrm{\mid }{P}_L\mathrm{\mid }$, $q=\mathrm{\mid }{P}_U\mathrm{\mid }$, and $u=q-n$. The stored strict-lower slots of $L$ are exactly $P_L$, and the stored slots of $U$ are exactly $P_U$; the unit diagonal of $L$ is implicit. Thus the factors use exactly $\mathrm{\ell }+q=s$ numeric slots. Here a factor's stored-slot count retains a slot even if an update makes its value zero.

For each $(i,k)\in P_L$, let $R_{ik}=\{j>k:(i,j)\in P\},C_k=\{j>k:(k,j)\in P\},$ and set $I={\sum }_{(i,k)\in P_L}\mathrm{\mid }{R}_{ik}\cap C_k\mathrm{\mid }$ and $J={\sum }_{(i,k)\in P_L}(\mathrm{\mid }{R}_{ik}\mathrm{\mid }+\mathrm{\mid }{C}_k\mathrm{\mid })$. On a successful setup, the factorization performs exactly $\mathrm{\ell }$ divisions, $I$ multiplications, and $I$ subtractions. The two-pointer merges perform at most $J$ column-index comparisons, in addition to one validation visit for each of the $n$ row intervals and the $s$ stored entries.

Applying a successful factor pair to one vector by the specified two triangular solves performs exactly $\mathrm{\ell }+u=s-n$ multiplications, $s-n$ subtractions, and $n$ divisions, and returns the unique vector produced by forward and backward substitution. These counts concern setup and one application separately. They assert no convergence, conditioning improvement, or accuracy guarantee. The stored-pattern conclusion is about $L$ and $U$; it does not assert that the product $LU$ has support contained in $P$.

Let $n,m\in \mathbb{N}$ with $n,m\ge 1$ (The natural numbers $\mathbb{N}$ (von Neumann)), let $A\in M_n(\mathbb{C})$ (The vector space $M_{m\times n}(F):=F^{m\times n}$ of $m$ by $n$ matrices over a field, with entrywise operations), and let $b,x_0\in {\mathbb{C}}^n$. The complex scalars form a field ($\mathbb{C}=\mathbb{R}[x]\mathrm{/}(x^2+1)$ is a field, every element is uniquely $a+bi$, and every nonzero element has inverse $(a-bi)\mathrm{/}(a^2+b^2)$). As additional input, take maps ${\mathcal{P}}_1,\dots ,{\mathcal{P}}_m:{\mathbb{C}}^n\to {\mathbb{C}}^n$. The maps may vary with $j$; no linearity or fixedness is assumed. In the usual matrix case, ${\mathcal{P}}_j(v)=M_j^{-1}v$ for an invertible $M_j$ (Invertible matrices and the general linear group ${\mathrm{GL}}_n(F)$). Equip each coordinate space ${\mathbb{C}}^q$ used below with the standard Hermitian inner product $⟨z,w⟩={\sum }_{a<q}{z}_a\overline{w_a}$ and its induced Euclidean norm $\mathrm{\parallel }z{\mathrm{\parallel }}_2=\sqrt{⟨z,z⟩}$ (The standard formulas $⟨x,y⟩={\sum }_{k<n}{x}_k{y}_k$ on ${\mathbb{R}}^n$ and ${\sum }_{k<n}{x}_k\overline{y_k}$ on ${\mathbb{C}}^n$ are inner products, The norm $\parallel v\parallel =\sqrt{⟨v,v⟩}$ induced by a real or complex inner product). All arithmetic and orthogonalization are exact. Iteration and auxiliary Hessenberg labels below start at $1$; coordinates of $A$ and ambient vectors retain the library's zero-based convention.

The algorithm uses and stores each actual $z_j$ and makes no fixed-matrix assumption.

Set $r_0=b-A{x}_0$ and $\beta =\mathrm{\parallel }{r}_0{\mathrm{\parallel }}_2$. If $\beta =0$, return initial-solution($x_0$). Otherwise set $v_1=r_0\mathrm{/}\beta$. For $j=1,\dots ,m$, perform:

1. Obtain $z_j={\mathcal{P}}_j(v_j)$ and set $w:=A{z}_j$.
2. For $i=1,\dots ,j$ in order, set $h_{ij}:=⟨w,v_i⟩,w:=w-h_{ij}{v}_i\mathrm{.}$
3. Set $h_{j+1,j}:=\mathrm{\parallel }w{\mathrm{\parallel }}_2$. If $h_{j+1,j}=0$, return breakdown($j,x_{j-1}$) with the last fully formed iterate (for $j=1$ it is $x_0$); do not normalize or claim a new iterate at this step.
4. If $h_{j+1,j}>0$, set $v_{j+1}:=w\mathrm{/}{h}_{j+1,j}$. Let $e_0\in {\mathbb{C}}^{j+1}$ be the standard unit vector at coordinate $0$ (coordinate $0$ is $1$ and all other coordinates are $0$; The standard list $e:n\to F^n$ with $e_i(i)=1_F$ and $e_i(j)=0_F$ for $j\mathrm{\ne }i$ is an ordered basis of $F^n$; hence ${\dim }_F{F}^n=n$, and $F^0$ is the zero space with basis $\mathrm{\varnothing }$ and dimension $0$). The iteration labels $h_{ik}$ remain one-based, while matrix slots are zero-based: put $Z_j\in M_{n\times j}(\mathbb{C})$ with $(Z_j{)}_{a,k-1}=(z_k{)}_a$ for $0\le a<n$, $1\le k\le j$, and let ${\bar{H}}_j\in M_{(j+1)\times j}(\mathbb{C})$ satisfy $({\bar{H}}_j{)}_{i-1,k-1}=h_{ik}$ for $1\le k\le j$, $1\le i\le k+1$, and $({\bar{H}}_j{)}_{a,b}=0$ when $a>b+1$. Choose any $y_j\in {\mathrm{argmin}}_{y\in {\mathbb{C}}^j}\mathrm{\parallel }\beta e_0-{\bar{H}}_{j}y{\mathrm{\parallel }}_2$ and set $x_j:=x_0+Z_j{y}_j$.

If all $m$ steps complete without breakdown, return budget-complete($m,x_m,Z_m$). The finite-dimensional least-squares minimizer exists because each coordinate space ${\mathbb{C}}^q$ is finite-dimensional by The standard list $e:n\to F^n$ with $e_i(i)=1_F$ and $e_i(j)=0_F$ for $j\mathrm{\ne }i$ is an ordered basis of $F^n$; hence ${\dim }_F{F}^n=n$, and $F^0$ is the zero space with basis $\mathrm{\varnothing }$ and dimension $0$, and the least-squares theorem For a linear map $T:V\to W$ between finite-dimensional inner-product spaces, $x$ minimises $\parallel Tx-b\parallel$ if and only if $T^{\ast }(Tx-b)=0$, equivalently $T^{\ast }Tx=T^{\ast }b$; minimisers exist and any two differ by an element of $\ker T$ applies; the output relation permits any minimizer. The algorithm has no convergence test: budget-complete and breakdown assert neither convergence nor that a breakdown iterate solves $Ax=b$. The residual-minimization property over the actual span of the $z_j$ is proved separately in Flexible GMRES minimizes over the actual preconditioned-vector span. Every division is by a positive real norm: the initial $\beta =0$ case returns before normalizing $r_0$, and each later $h_{j+1,j}=0$ case returns before normalizing $w$.

Let $A\in M_n(\mathbb{C})$ and $b,x_0\in {\mathbb{C}}^n$, with the standard Hermitian inner product and induced Euclidean norm. Suppose the algorithm Flexible GMRES with variable right preconditioning completes $m\ge 1$ steps without breakdown, producing $Z_m=[z_1,\dots ,z_m]$ and ${\bar{H}}_m$. Define $\beta =\mathrm{\parallel }b-A{x}_0{\mathrm{\parallel }}_2$, let $e_0\in {\mathbb{C}}^{m+1}$ be the standard unit vector at coordinate $0$, and define $V_{m+1}\in M_{n\times (m+1)}(\mathbb{C})$ by $(V_{m+1}{)}_{a,i-1}=(v_i{)}_a$ for $0\le a<n$, $1\le i\le m+1$. The one-based column label $i$ therefore occupies zero-based matrix column $i-1$. The matrices $Z_m$ and ${\bar{H}}_m$ use the algorithm's one-based-to-zero-based slot mapping. Then $A{Z}_m=V_{m+1}{\bar{H}}_m,$ and every $y\in {\mathbb{C}}^m$ satisfies $b-A(x_0+Z_{m}y)=V_{m+1}(\beta e_0-{\bar{H}}_{m}y),\mathrm{\parallel }b-A(x_0+Z_{m}y){\mathrm{\parallel }}_2=\mathrm{\parallel }\beta e_0-{\bar{H}}_{m}y{\mathrm{\parallel }}_2\mathrm{.}$ Consequently every least-squares minimizer $y_m$ in the algorithm produces a true-residual minimizer $x_m=x_0+Z_m{y}_m$ over $x_0+\mathrm{span}\{z_1,\dots ,z_m\}\mathrm{.}$ Moreover, ${\bar{H}}_m$ has full column rank, so its least-squares minimizer is unique. The result makes no fixed-preconditioner spectral or convergence claim.