The SPIKE algorithm is a hybrid parallel solver for banded linear systems developed by Eric Polizzi and Ahmed Sameh[#endnote_] ;templatestyles src="Citation/styles.css"/ ;^ [#endnote_]
Overview.
The SPIKE algorithm deals with a linear system AX 
 F, where A is a banded formula_1 matrix of bandwidth much less than formula_2, and F is an formula_3 matrix containing formula_4 right-hand sides. It is divided into a preprocessing stage and a postprocessing stage.
Preprocessing stage.
In the preprocessing stage, the linear system AX 
 F is partitioned into a block tridiagonal form
formula_5
Assume, for the time being, that the diagonal blocks Aj (j 
 1...,p with p ≥ 2) are nonsingular. Define a block diagonal matrix
D 
 diag(A1...,Ap),
then D is also nonsingular. Left-multiplying D−1 to both sides of the system gives
formula_6
which is to be solved in the postprocessing stage. Left-multiplication by D−1 is equivalent to solving formula_7 systems of the form
Aj[Vj Wj Gj] 
 [Bj Cj Fj]
(omitting W1 and C1 for formula_8, and Vp and Bp for formula_9), which can be carried out in parallel.
Due to the banded nature of A, only a few leftmost columns of each Vj and a few rightmost columns of each Wj can be nonzero. These columns are called the "spikes".
Postprocessing stage.
Without loss of generality, assume that each spike contains exactly formula_10 columns (formula_10 is much less than formula_2) (pad the spike with columns of zeroes if necessary). Partition the spikes in all Vj and Wj into
formula_13 and formula_14
where , , and are of dimensions formula_15. Partition similarly all Xj and Gj into
formula_16 and formula_17
Notice that the system produced by the preprocessing stage can be reduced to a block pentadiagonal system of much smaller size (recall that formula_10 is much less than formula_2)
formula_20
which we call the "reduced system" and denote by S̃X̃ 
 G̃.
Once all and are found, all X′j can be recovered with perfect parallelism via
formula_21
SPIKE as a polyalgorithmic banded linear system solver.
Despite being logically divided into two stages, computationally, the SPIKE algorithm comprises three stages:
Each of these stages can be accomplished in several ways, allowing a multitude of variants. Two notable variants are the "recursive SPIKE" algorithm for non-diagonally-dominant cases and the "truncated SPIKE" algorithm for diagonally-dominant cases. Depending on the variant, a system can be solved either exactly or approximately. In the latter case, SPIKE is used as a preconditioner for iterative schemes like Krylov subspace methods and iterative refinement.
Recursive SPIKE.
Preprocessing stage.
The first step of the preprocessing stage is to factorize the diagonal blocks Aj. For numerical stability, one can use LAPACK's codice_1 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‖1.
If the pivot does not satisfy the condition, it is then boosted by
formula_22
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 codice_2 routines. After the diagonal blocks are factorized, the spikes are computed and passed on to the postprocessing stage.
Postprocessing stage.
The two-partition case.
In the two-partition case, i.e., when p 
 2, the reduced system S̃X̃ 
 G̃ has the form
formula_23
An even smaller system can be extracted from the center:
formula_24
which can be solved using the block LU factorization
formula_25
Once and are found, and can be computed via
 − ,
 − .
The multiple-partition case.
Assume that p is a power of two, i.e., p 
 2d. Consider a block diagonal matrix
D̃1 
 diag(...,)
where
formula_26
for k 
 1...,p/2. Notice that D̃1 essentially consists of diagonal blocks of order 4m extracted from S̃. Now we factorize S̃ as
S̃ 
 D̃1S̃2.
The new matrix S̃2 has the form
formula_27
Its structure is very similar to that of S̃2, 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̃2 to produce
S̃2 
 D̃2S̃3
and
S̃ 
 D̃1D̃2S̃3.
Such factorization steps can be performed recursively. After d − 1 steps, we obtain the factorization
S̃ 
 D̃1⋯D̃d−1S̃d,
where S̃d has only two spikes. The reduced system will then be solved via
X̃ 
 ⋯G̃.
The block LU factorization technique in the two-partition case can be used to handle the solving steps involving D̃1, ..., 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.
Truncated SPIKE.
When A is diagonally-dominant, in the reduced system
formula_20
the blocks and are often negligible. With them omitted, the reduced system becomes block diagonal
formula_29
and can be easily solved in parallel [#endnote_].
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.
SPIKE for tridiagonal systems.
The first SPIKE partitioning and algorithm was presented in [#endnote_] and was designed as the means to improve the stability properties of a parallel Givens rotations-based solver for tridiagonal systems. A version of the algorithm, termed g-Spike, that is based on serial Givens rotations applied independently on each block was designed for the NVIDIA GPU [#endnote_]. A SPIKE-based algorithm for the GPU that is based on a special block diagonal pivoting strategy is described in [#endnote_].
SPIKE as a preconditioner.
The SPIKE algorithm can also function as a preconditioner for iterative methods for solving linear systems. To solve a linear system Ax 
 b using a SPIKE-preconditioned iterative solver, one extracts center bands from A to form a banded preconditioner M and solves linear systems involving M in each iteration with the SPIKE algorithm.
In order for the preconditioner to be effective, row and/or column permutation is usually necessary to move "heavy" elements of A close to the diagonal so that they are covered by the preconditioner. This can be accomplished by computing the weighted spectral reordering of A.
The SPIKE algorithm can be generalized by not restricting the preconditioner to be strictly banded. In particular, the diagonal block in each partition can be a general matrix and thus handled by a direct general linear system solver rather than a banded solver. This enhances the preconditioner, and hence allows better chance of convergence and reduces the number of iterations.
Implementations.
Intel offers an implementation of the SPIKE algorithm under the name "Intel Adaptive Spike-Based Solver" [#endnote_]. Tridiagonal solvers have also been developed for the NVIDIA GPU 
[#endnote_][#endnote_] and the Xeon Phi co-processors. The method in [#endnote_] is the basis for a tridiagonal solver in the cuSPARSE library. The Givens rotations based solver was also implemented for the 
GPU and the Intel Xeon Phi.
References.
 ;templatestyles src="Reflist/styles.css" / ;