SPIKE algorithm

In the preprocessing stage, the linear system AX = F is partitioned into a block tridiagonal form

${\displaystyle {\begin{bmatrix}{\boldsymbol {A}}_{1}&{\boldsymbol {B}}_{1}\\{\boldsymbol {C}}_{2}&{\boldsymbol {A}}_{2}&{\boldsymbol {B}}_{2}\\&\ddots &\ddots &\ddots \\&&{\boldsymbol {C}}_{p-1}&{\boldsymbol {A}}_{p-1}&{\boldsymbol {B}}_{p-1}\\&&&{\boldsymbol {C}}_{p}&{\boldsymbol {A}}_{p}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}\\{\boldsymbol {X}}_{2}\\\vdots \\{\boldsymbol {X}}_{p-1}\\{\boldsymbol {X}}_{p}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {F}}_{1}\\{\boldsymbol {F}}_{2}\\\vdots \\{\boldsymbol {F}}_{p-1}\\{\boldsymbol {F}}_{p}\end{bmatrix}}.}$

Assume, for the time being, that the diagonal blocks A_j (j = 1,...,p with p ≥ 2) are nonsingular. Define a block diagonal matrix

D = diag(A₁,...,A_p),

then D is also nonsingular. Left-multiplying D⁻¹ to both sides of the system gives

${\displaystyle {\begin{bmatrix}{\boldsymbol {I}}&{\boldsymbol {V}}_{1}\\{\boldsymbol {W}}_{2}&{\boldsymbol {I}}&{\boldsymbol {V}}_{2}\\&\ddots &\ddots &\ddots \\&&{\boldsymbol {W}}_{p-1}&{\boldsymbol {I}}&{\boldsymbol {V}}_{p-1}\\&&&{\boldsymbol {W}}_{p}&{\boldsymbol {I}}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}\\{\boldsymbol {X}}_{2}\\\vdots \\{\boldsymbol {X}}_{p-1}\\{\boldsymbol {X}}_{p}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}\\{\boldsymbol {G}}_{2}\\\vdots \\{\boldsymbol {G}}_{p-1}\\{\boldsymbol {G}}_{p}\end{bmatrix}},}$

which is to be solved in the postprocessing stage. Left-multiplication by D⁻¹ is equivalent to solving ${\displaystyle p}$ systems of the form

A_j[V_j W_j G_j] = [B_j C_j F_j]

(omitting W₁ and C₁ for ${\displaystyle j=1}$, and V_p and B_p for ${\displaystyle j=p}$), which can be carried out in parallel.

Due to the banded nature of A, only a few leftmost columns of each V_j and a few rightmost columns of each W_j can be nonzero. These columns are called the spikes.

Without loss of generality, assume that each spike contains exactly ${\displaystyle m}$ columns (${\displaystyle m}$ is much less than ${\displaystyle n}$) (pad the spike with columns of zeroes if necessary). Partition the spikes in all V_j and W_j into

${\displaystyle {\begin{bmatrix}{\boldsymbol {V}}_{j}^{(t)}\\{\boldsymbol {V}}_{j}'\\{\boldsymbol {V}}_{j}^{(b)}\end{bmatrix}}}$ and ${\displaystyle {\begin{bmatrix}{\boldsymbol {W}}_{j}^{(t)}\\{\boldsymbol {W}}_{j}'\\{\boldsymbol {W}}_{j}^{(b)}\\\end{bmatrix}}}$

where V (t)
j , V (b)
j , W (t)
j  and W (b)
j  are of dimensions ${\displaystyle m\times m}$. Partition similarly all X_j and G_j into

${\displaystyle {\begin{bmatrix}{\boldsymbol {X}}_{j}^{(t)}\\{\boldsymbol {X}}_{j}'\\{\boldsymbol {X}}_{j}^{(b)}\end{bmatrix}}}$ and ${\displaystyle {\begin{bmatrix}{\boldsymbol {G}}_{j}^{(t)}\\{\boldsymbol {G}}_{j}'\\{\boldsymbol {G}}_{j}^{(b)}\\\end{bmatrix}}.}$

Notice that the system produced by the preprocessing stage can be reduced to a block pentadiagonal system of much smaller size (recall that ${\displaystyle m}$ is much less than ${\displaystyle n}$)

${\displaystyle {\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{1}^{(t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{2}^{(t)}\\&{\boldsymbol {W}}_{2}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{2}^{(b)}&{\boldsymbol {0}}\\&&\ddots &\ddots &\ddots &\ddots &\ddots \\&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p-1}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{p-1}^{(t)}\\&&&&{\boldsymbol {W}}_{p-1}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{p-1}^{(b)}&{\boldsymbol {0}}\\&&&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&&&&&&{\boldsymbol {W}}_{p}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(t)}\\{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\\{\boldsymbol {X}}_{2}^{(b)}\\\vdots \\{\boldsymbol {X}}_{p-1}^{(t)}\\{\boldsymbol {X}}_{p-1}^{(b)}\\{\boldsymbol {X}}_{p}^{(t)}\\{\boldsymbol {X}}_{p}^{(b)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(t)}\\{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\\{\boldsymbol {G}}_{2}^{(b)}\\\vdots \\{\boldsymbol {G}}_{p-1}^{(t)}\\{\boldsymbol {G}}_{p-1}^{(b)}\\{\boldsymbol {G}}_{p}^{(t)}\\{\boldsymbol {G}}_{p}^{(b)}\end{bmatrix}}{\text{,}}}$

which we call the reduced system and denote by S̃X̃ = G̃.

Once all X (t)
j  and X (b)
j  are found, all X′_j can be recovered with perfect parallelism via

${\displaystyle {\begin{cases}{\boldsymbol {X}}_{1}'={\boldsymbol {G}}_{1}'-{\boldsymbol {V}}_{1}'{\boldsymbol {X}}_{2}^{(t)}{\text{,}}\\{\boldsymbol {X}}_{j}'={\boldsymbol {G}}_{j}'-{\boldsymbol {V}}_{j}'{\boldsymbol {X}}_{j+1}^{(t)}-{\boldsymbol {W}}_{j}'{\boldsymbol {X}}_{j-1}^{(b)}{\text{,}}&j=2,\ldots ,p-1{\text{,}}\\{\boldsymbol {X}}_{p}'={\boldsymbol {G}}_{p}'-{\boldsymbol {W}}_{p}{\boldsymbol {X}}_{p-1}^{(b)}{\text{.}}\end{cases}}}$

The first step of the preprocessing stage is to factorize the diagonal blocks A_j. For numerical stability, one can use LAPACK's XGBTRF routines to LU factorize them with partial pivoting. Alternatively, one can also factorize them without partial pivoting but with a "diagonal boosting" strategy. The latter method tackles the issue of singular diagonal blocks.

In concrete terms, the diagonal boosting strategy is as follows. Let 0_ε denote a configurable "machine zero". In each step of LU factorization, we require that the pivot satisfy the condition

|pivot| > 0_ε‖A‖₁.

If the pivot does not satisfy the condition, it is then boosted by

${\displaystyle \mathrm {pivot} ={\begin{cases}\mathrm {pivot} +\epsilon \lVert {\boldsymbol {A}}_{j}\rVert _{1}&{\text{if }}\mathrm {pivot} \geq 0{\text{,}}\\\mathrm {pivot} -\epsilon \lVert {\boldsymbol {A}}_{j}\rVert _{1}&{\text{if }}\mathrm {pivot} <0\end{cases}}}$

where ε is a positive parameter depending on the machine's unit roundoff, and the factorization continues with the boosted pivot. This can be achieved by modified versions of ScaLAPACK's XDBTRF routines. After the diagonal blocks are factorized, the spikes are computed and passed on to the postprocessing stage.

In the two-partition case, i.e., when p = 2, the reduced system S̃X̃ = G̃ has the form

${\displaystyle {\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{1}^{(t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&{\boldsymbol {W}}_{2}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(t)}\\{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\\{\boldsymbol {X}}_{2}^{(b)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(t)}\\{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\\{\boldsymbol {G}}_{2}^{(b)}\end{bmatrix}}{\text{.}}}$

An even smaller system can be extracted from the center:

${\displaystyle {\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}\\{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\end{bmatrix}}{\text{,}}}$

which can be solved using the block LU factorization

${\displaystyle {\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}\\{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {I}}_{m}\\{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}\\&{\boldsymbol {I}}_{m}-{\boldsymbol {W}}_{2}^{(t)}{\boldsymbol {V}}_{1}^{(b)}\end{bmatrix}}{\text{.}}}$

Once X (b)
1  and X (t)
2  are found, X (t)
1  and X (b)
2  can be computed via

X (t)
1  = G (t)
1  − V (t)
1 X (t)
2 ,

X (b)
2  = G (b)
2  − W (b)
2 X (b)
1 .

Assume that p is a power of two, i.e., p = 2^d. Consider a block diagonal matrix

D̃₁ = diag(D̃ [1]
1 ,...,D̃ [1]
p/2 )

where

${\displaystyle {\boldsymbol {\tilde {D}}}_{k}^{[1]}={\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{2k-1}^{(t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{2k-1}^{(b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2k}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&{\boldsymbol {W}}_{2k}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}\end{bmatrix}}}$

for k = 1,...,p/2. Notice that D̃₁ essentially consists of diagonal blocks of order 4m extracted from S̃. Now we factorize S̃ as

S̃ = D̃₁S̃₂.

The new matrix S̃₂ has the form

${\displaystyle {\begin{bmatrix}{\boldsymbol {I}}_{3m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{1}^{[2](t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{[2](b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2}^{[2](t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{2}^{[2](t)}\\&{\boldsymbol {W}}_{2}^{[2](b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{3m}&{\boldsymbol {V}}_{2}^{[2](b)}&{\boldsymbol {0}}\\&&\ddots &\ddots &\ddots &\ddots &\ddots \\&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p/2-1}^{[2](t)}&{\boldsymbol {I}}_{3m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{p/2-1}^{[2](t)}\\&&&&{\boldsymbol {W}}_{p/2-1}^{[2](b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{p/2-1}^{[2](b)}&{\boldsymbol {0}}\\&&&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p/2}^{[2](t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&&&&&&{\boldsymbol {W}}_{p/2}^{[2](b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{3m}\end{bmatrix}}{\text{.}}}$

Its structure is very similar to that of S̃₂, only differing in the number of spikes and their height (their width stays the same at m). Thus, a similar factorization step can be performed on S̃₂ to produce

S̃₂ = D̃₂S̃₃

and

S̃ = D̃₁D̃₂S̃₃.

Such factorization steps can be performed recursively. After d − 1 steps, we obtain the factorization

S̃ = D̃₁⋯D̃_d−1S̃_d,

where S̃_d has only two spikes. The reduced system will then be solved via

X̃ = S̃ −1
d D̃ −1
d−1 ⋯D̃ −1
1 G̃.

The block LU factorization technique in the two-partition case can be used to handle the solving steps involving D̃₁, ..., D̃_d−1 and S̃_d for they essentially solve multiple independent systems of generalized two-partition forms.

Generalization to cases where p is not a power of two is almost trivial.

When A is diagonally-dominant, in the reduced system

${\displaystyle {\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{1}^{(t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{2}^{(t)}\\&{\boldsymbol {W}}_{2}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{2}^{(b)}&{\boldsymbol {0}}\\&&\ddots &\ddots &\ddots &\ddots &\ddots \\&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p-1}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{p-1}^{(t)}\\&&&&{\boldsymbol {W}}_{p-1}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{p-1}^{(b)}&{\boldsymbol {0}}\\&&&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&&&&&&{\boldsymbol {W}}_{p}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(t)}\\{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\\{\boldsymbol {X}}_{2}^{(b)}\\\vdots \\{\boldsymbol {X}}_{p-1}^{(t)}\\{\boldsymbol {X}}_{p-1}^{(b)}\\{\boldsymbol {X}}_{p}^{(t)}\\{\boldsymbol {X}}_{p}^{(b)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(t)}\\{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\\{\boldsymbol {G}}_{2}^{(b)}\\\vdots \\{\boldsymbol {G}}_{p-1}^{(t)}\\{\boldsymbol {G}}_{p-1}^{(b)}\\{\boldsymbol {G}}_{p}^{(t)}\\{\boldsymbol {G}}_{p}^{(b)}\end{bmatrix}}{\text{,}}}$

the blocks V (t)
j  and W (b)
j  are often negligible. With them omitted, the reduced system becomes block diagonal

${\displaystyle {\begin{bmatrix}{\boldsymbol {I}}_{m}\\&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}\\&{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}\\&&&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{2}^{(b)}\\&&&\ddots &\ddots &\ddots \\&&&&{\boldsymbol {W}}_{p-1}^{(t)}&{\boldsymbol {I}}_{m}\\&&&&&&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{p-1}^{(b)}\\&&&&&&{\boldsymbol {W}}_{p}^{(t)}&{\boldsymbol {I}}_{m}\\&&&&&&&&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(t)}\\{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\\{\boldsymbol {X}}_{2}^{(b)}\\\vdots \\{\boldsymbol {X}}_{p-1}^{(t)}\\{\boldsymbol {X}}_{p-1}^{(b)}\\{\boldsymbol {X}}_{p}^{(t)}\\{\boldsymbol {X}}_{p}^{(b)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(t)}\\{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\\{\boldsymbol {G}}_{2}^{(b)}\\\vdots \\{\boldsymbol {G}}_{p-1}^{(t)}\\{\boldsymbol {G}}_{p-1}^{(b)}\\{\boldsymbol {G}}_{p}^{(t)}\\{\boldsymbol {G}}_{p}^{(b)}\end{bmatrix}}}$

and can be easily solved in parallel .

The truncated SPIKE algorithm can be wrapped inside some outer iterative scheme (e.g., BiCGSTAB or iterative refinement) to improve the accuracy of the solution.