SPIKE algorithm

The SPIKE algorithm is a hybrid parallel solver for banded linear systems developed by Eric Polizzi and Ahmed SamehTemplate:RefTemplate:NoteTemplate:Ref

Overview

The SPIKE algorithm deals with a linear system Template:Math, where Template:Math is a banded $n\times n$ matrix of bandwidth much less than $n$, and Template:Math is an $n\times s$ matrix containing $s$ right-hand sides. It is divided into a preprocessing stage and a postprocessing stage.

Preprocessing stage

In the preprocessing stage, the linear system Template:Math is partitioned into a block tridiagonal form

$\left[\begin{array}{cc}{𝑨}_1& {𝑩}_1\\ {𝑪}_2& {𝑨}_2& {𝑩}_2\\ & \ddots & \ddots & \ddots \\ & & {𝑪}_{p-1}& {𝑨}_{p-1}& {𝑩}_{p-1}\\ & & & {𝑪}_p& {𝑨}_p\end{array}\right]\left[\begin{array}{c}{𝑿}_1\\ {𝑿}_2\\ ⋮\\ {𝑿}_{p-1}\\ {𝑿}_p\end{array}\right]=\left[\begin{array}{c}{𝑭}_1\\ {𝑭}_2\\ ⋮\\ {𝑭}_{p-1}\\ {𝑭}_p\end{array}\right].$

Assume, for the time being, that the diagonal blocks Template:Math (Template:Math with Template:Math) are nonsingular. Define a block diagonal matrix

Template:Math,

then Template:Math is also nonsingular. Left-multiplying Template:Math to both sides of the system gives

$\left[\begin{array}{cc}𝑰& {𝑽}_1\\ {𝑾}_2& 𝑰& {𝑽}_2\\ & \ddots & \ddots & \ddots \\ & & {𝑾}_{p-1}& 𝑰& {𝑽}_{p-1}\\ & & & {𝑾}_p& 𝑰\end{array}\right]\left[\begin{array}{c}{𝑿}_1\\ {𝑿}_2\\ ⋮\\ {𝑿}_{p-1}\\ {𝑿}_p\end{array}\right]=\left[\begin{array}{c}{𝑮}_1\\ {𝑮}_2\\ ⋮\\ {𝑮}_{p-1}\\ {𝑮}_p\end{array}\right],$

which is to be solved in the postprocessing stage. Left-multiplication by Template:Math is equivalent to solving $p$ systems of the form

Template:Math

(omitting Template:Math and Template:Math for $j=1$, and Template:Math and Template:Math for $j=p$), which can be carried out in parallel.

Due to the banded nature of Template:Math, only a few leftmost columns of each Template:Math and a few rightmost columns of each Template:Math can be nonzero. These columns are called the spikes.

Postprocessing stage

Without loss of generality, assume that each spike contains exactly $m$ columns ($m$ is much less than $n$) (pad the spike with columns of zeroes if necessary). Partition the spikes in all Template:Math and Template:Math into

$\left[\begin{array}{c}{𝑽}_j^{(t)}\\ {{𝑽}_j}^{\prime }\\ {𝑽}_j^{(b)}\end{array}\right]$ and $\left[\begin{array}{c}{𝑾}_j^{(t)}\\ {{𝑾}_j}^{\prime }\\ {𝑾}_j^{(b)}\end{array}\right]$

where Template:Math, Template:Math, Template:Math and Template:Math are of dimensions $m\times m$. Partition similarly all Template:Math and Template:Math into

$\left[\begin{array}{c}{𝑿}_j^{(t)}\\ {{𝑿}_j}^{\prime }\\ {𝑿}_j^{(b)}\end{array}\right]$ and $\left[\begin{array}{c}{𝑮}_j^{(t)}\\ {{𝑮}_j}^{\prime }\\ {𝑮}_j^{(b)}\end{array}\right].$

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

$\left[\begin{array}{ccc}{𝑰}_m& 0& {𝑽}_1^{(t)}\\ 0& {𝑰}_m& {𝑽}_1^{(b)}& 0\\ 0& {𝑾}_2^{(t)}& {𝑰}_m& 0& {𝑽}_2^{(t)}\\ & {𝑾}_2^{(b)}& 0& {𝑰}_m& {𝑽}_2^{(b)}& 0\\ & & \ddots & \ddots & \ddots & \ddots & \ddots \\ & & & 0& {𝑾}_{p-1}^{(t)}& {𝑰}_m& 0& {𝑽}_{p-1}^{(t)}\\ & & & & {𝑾}_{p-1}^{(b)}& 0& {𝑰}_m& {𝑽}_{p-1}^{(b)}& 0\\ & & & & & 0& {𝑾}_p^{(t)}& {𝑰}_m& 0\\ & & & & & & {𝑾}_p^{(b)}& 0& {𝑰}_m\end{array}\right]\left[\begin{array}{c}{𝑿}_1^{(t)}\\ {𝑿}_1^{(b)}\\ {𝑿}_2^{(t)}\\ {𝑿}_2^{(b)}\\ ⋮\\ {𝑿}_{p-1}^{(t)}\\ {𝑿}_{p-1}^{(b)}\\ {𝑿}_p^{(t)}\\ {𝑿}_p^{(b)}\end{array}\right]=\left[\begin{array}{c}{𝑮}_1^{(t)}\\ {𝑮}_1^{(b)}\\ {𝑮}_2^{(t)}\\ {𝑮}_2^{(b)}\\ ⋮\\ {𝑮}_{p-1}^{(t)}\\ {𝑮}_{p-1}^{(b)}\\ {𝑮}_p^{(t)}\\ {𝑮}_p^{(b)}\end{array}\right]\text{,}$

which we call the reduced system and denote by Template:Math.

Once all Template:Math and Template:Math are found, all Template:Math can be recovered with perfect parallelism via

$\{\begin{array}{c}{{𝑿}_1}^{\prime }={{𝑮}_1}^{\prime }-{{𝑽}_1}^{\prime }{𝑿}_2^{(t)}\text{,}\\ {{𝑿}_j}^{\prime }={{𝑮}_j}^{\prime }-{{𝑽}_j}^{\prime }{𝑿}_{j+1}^{(t)}-{{𝑾}_j}^{\prime }{𝑿}_{j-1}^{(b)}\text{,}& j=2,\dots ,p-1\text{,}\\ {{𝑿}_p}^{\prime }={{𝑮}_p}^{\prime }-{𝑾}_p{𝑿}_{p-1}^{(b)}\text{.}\end{array}$

SPIKE as a polyalgorithmic banded linear system solver

Despite being logically divided into two stages, computationally, the SPIKE algorithm comprises three stages:

1. factorizing the diagonal blocks,
2. computing the spikes,
3. solving the reduced system.

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 Template:Math. 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 Template:Math denote a configurable "machine zero". In each step of LU factorization, we require that the pivot satisfy the condition

Template:Math.

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

$\mathrm{p}\mathrm{i}\mathrm{v}\mathrm{o}\mathrm{t}=\{\begin{array}{cc}\mathrm{p}\mathrm{i}\mathrm{v}\mathrm{o}\mathrm{t}+ϵ\Vert {𝑨}_j{\Vert }_1& \text{if }\mathrm{p}\mathrm{i}\mathrm{v}\mathrm{o}\mathrm{t}\ge 0\text{,}\\ \mathrm{p}\mathrm{i}\mathrm{v}\mathrm{o}\mathrm{t}-ϵ\Vert {𝑨}_j{\Vert }_1& \text{if }\mathrm{p}\mathrm{i}\mathrm{v}\mathrm{o}\mathrm{t}<0\end{array}$

where Template:Math 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.

Postprocessing stage

The two-partition case

In the two-partition case, i.e., when Template:Math, the reduced system Template:Math has the form

$\left[\begin{array}{ccc}{𝑰}_m& 0& {𝑽}_1^{(t)}\\ 0& {𝑰}_m& {𝑽}_1^{(b)}& 0\\ 0& {𝑾}_2^{(t)}& {𝑰}_m& 0\\ & {𝑾}_2^{(b)}& 0& {𝑰}_m\end{array}\right]\left[\begin{array}{c}{𝑿}_1^{(t)}\\ {𝑿}_1^{(b)}\\ {𝑿}_2^{(t)}\\ {𝑿}_2^{(b)}\end{array}\right]=\left[\begin{array}{c}{𝑮}_1^{(t)}\\ {𝑮}_1^{(b)}\\ {𝑮}_2^{(t)}\\ {𝑮}_2^{(b)}\end{array}\right]\text{.}$

An even smaller system can be extracted from the center:

$\left[\begin{array}{cc}{𝑰}_m& {𝑽}_1^{(b)}\\ {𝑾}_2^{(t)}& {𝑰}_m\end{array}\right]\left[\begin{array}{c}{𝑿}_1^{(b)}\\ {𝑿}_2^{(t)}\end{array}\right]=\left[\begin{array}{c}{𝑮}_1^{(b)}\\ {𝑮}_2^{(t)}\end{array}\right]\text{,}$

which can be solved using the block LU factorization

$\left[\begin{array}{cc}{𝑰}_m& {𝑽}_1^{(b)}\\ {𝑾}_2^{(t)}& {𝑰}_m\end{array}\right]=\left[\begin{array}{c}{𝑰}_m\\ {𝑾}_2^{(t)}& {𝑰}_m\end{array}\right]\left[\begin{array}{cc}{𝑰}_m& {𝑽}_1^{(b)}\\ & {𝑰}_m-{𝑾}_2^{(t)}{𝑽}_1^{(b)}\end{array}\right]\text{.}$

Once Template:Math and Template:Math are found, Template:Math and Template:Math can be computed via

Template:Math,

Template:Math.

The multiple-partition case

Assume that Template:Math is a power of two, i.e., Template:Math. Consider a block diagonal matrix

Template:Math

where

${\widetilde{𝑫}}_k^{[1]}=\left[\begin{array}{ccc}{𝑰}_m& 0& {𝑽}_{2k-1}^{(t)}\\ 0& {𝑰}_m& {𝑽}_{2k-1}^{(b)}& 0\\ 0& {𝑾}_{2k}^{(t)}& {𝑰}_m& 0\\ & {𝑾}_{2k}^{(b)}& 0& {𝑰}_m\end{array}\right]$

for Template:Math. Notice that Template:Math essentially consists of diagonal blocks of order Template:Math extracted from Template:Math. Now we factorize Template:Math as

Template:Math.

The new matrix Template:Math has the form

$\left[\begin{array}{ccc}{𝑰}_{3m}& 0& {𝑽}_1^{[2](t)}\\ 0& {𝑰}_m& {𝑽}_1^{[2](b)}& 0\\ 0& {𝑾}_2^{[2](t)}& {𝑰}_m& 0& {𝑽}_2^{[2](t)}\\ & {𝑾}_2^{[2](b)}& 0& {𝑰}_{3m}& {𝑽}_2^{[2](b)}& 0\\ & & \ddots & \ddots & \ddots & \ddots & \ddots \\ & & & 0& {𝑾}_{p/2-1}^{[2](t)}& {𝑰}_{3m}& 0& {𝑽}_{p/2-1}^{[2](t)}\\ & & & & {𝑾}_{p/2-1}^{[2](b)}& 0& {𝑰}_m& {𝑽}_{p/2-1}^{[2](b)}& 0\\ & & & & & 0& {𝑾}_{p/2}^{[2](t)}& {𝑰}_m& 0\\ & & & & & & {𝑾}_{p/2}^{[2](b)}& 0& {𝑰}_{3m}\end{array}\right]\text{.}$

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

Template:Math

and

Template:Math.

Such factorization steps can be performed recursively. After Template:Math steps, we obtain the factorization

Template:Math,

where Template:Math has only two spikes. The reduced system will then be solved via

Template:Math.

The block LU factorization technique in the two-partition case can be used to handle the solving steps involving Template:Math, ..., Template:Math and Template:Math for they essentially solve multiple independent systems of generalized two-partition forms.

Generalization to cases where Template:Math is not a power of two is almost trivial.

Truncated SPIKE

When Template:Math is diagonally-dominant, in the reduced system

$\left[\begin{array}{ccc}{𝑰}_m& 0& {𝑽}_1^{(t)}\\ 0& {𝑰}_m& {𝑽}_1^{(b)}& 0\\ 0& {𝑾}_2^{(t)}& {𝑰}_m& 0& {𝑽}_2^{(t)}\\ & {𝑾}_2^{(b)}& 0& {𝑰}_m& {𝑽}_2^{(b)}& 0\\ & & \ddots & \ddots & \ddots & \ddots & \ddots \\ & & & 0& {𝑾}_{p-1}^{(t)}& {𝑰}_m& 0& {𝑽}_{p-1}^{(t)}\\ & & & & {𝑾}_{p-1}^{(b)}& 0& {𝑰}_m& {𝑽}_{p-1}^{(b)}& 0\\ & & & & & 0& {𝑾}_p^{(t)}& {𝑰}_m& 0\\ & & & & & & {𝑾}_p^{(b)}& 0& {𝑰}_m\end{array}\right]\left[\begin{array}{c}{𝑿}_1^{(t)}\\ {𝑿}_1^{(b)}\\ {𝑿}_2^{(t)}\\ {𝑿}_2^{(b)}\\ ⋮\\ {𝑿}_{p-1}^{(t)}\\ {𝑿}_{p-1}^{(b)}\\ {𝑿}_p^{(t)}\\ {𝑿}_p^{(b)}\end{array}\right]=\left[\begin{array}{c}{𝑮}_1^{(t)}\\ {𝑮}_1^{(b)}\\ {𝑮}_2^{(t)}\\ {𝑮}_2^{(b)}\\ ⋮\\ {𝑮}_{p-1}^{(t)}\\ {𝑮}_{p-1}^{(b)}\\ {𝑮}_p^{(t)}\\ {𝑮}_p^{(b)}\end{array}\right]\text{,}$

the blocks Template:Math and Template:Math are often negligible. With them omitted, the reduced system becomes block diagonal

$\left[\begin{array}{c}{𝑰}_m\\ & {𝑰}_m& {𝑽}_1^{(b)}\\ & {𝑾}_2^{(t)}& {𝑰}_m\\ & & & {𝑰}_m& {𝑽}_2^{(b)}\\ & & & \ddots & \ddots & \ddots \\ & & & & {𝑾}_{p-1}^{(t)}& {𝑰}_m\\ & & & & & & {𝑰}_m& {𝑽}_{p-1}^{(b)}\\ & & & & & & {𝑾}_p^{(t)}& {𝑰}_m\\ & & & & & & & & {𝑰}_m\end{array}\right]\left[\begin{array}{c}{𝑿}_1^{(t)}\\ {𝑿}_1^{(b)}\\ {𝑿}_2^{(t)}\\ {𝑿}_2^{(b)}\\ ⋮\\ {𝑿}_{p-1}^{(t)}\\ {𝑿}_{p-1}^{(b)}\\ {𝑿}_p^{(t)}\\ {𝑿}_p^{(b)}\end{array}\right]=\left[\begin{array}{c}{𝑮}_1^{(t)}\\ {𝑮}_1^{(b)}\\ {𝑮}_2^{(t)}\\ {𝑮}_2^{(b)}\\ ⋮\\ {𝑮}_{p-1}^{(t)}\\ {𝑮}_{p-1}^{(b)}\\ {𝑮}_p^{(t)}\\ {𝑮}_p^{(b)}\end{array}\right]$

and can be easily solved in parallel Template:Ref.

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 Template:Ref 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 Template:Ref. A SPIKE-based algorithm for the GPU that is based on a special block diagonal pivoting strategy is described in Template:Ref.

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 Template:Math using a SPIKE-preconditioned iterative solver, one extracts center bands from Template:Math to form a banded preconditioner Template:Math and solves linear systems involving Template:Math 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 Template:Math close to the diagonal so that they are covered by the preconditioner. This can be accomplished by computing the weighted spectral reordering of Template:Math.

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 Template:Ref. Tridiagonal solvers have also been developed for the NVIDIA GPU Template:RefTemplate:Ref and the Xeon Phi co-processors. The method in Template:Ref 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

1. Template:NoteTemplate:Cite journal
2. Template:NoteTemplate:Cite journal
3. Template:NoteTemplate:Cite journal
4. Template:NoteTemplate:Cite book
5. Template:NoteTemplate:Cite web
6. Template:NoteTemplate:Cite journal
7. Template:NoteTemplate:Cite journal
8. Template:NoteTemplate:Cite journal

Further reading

• Template:Cite book