Parallel Solution of Hierarchical Symmetric Positive Definite Linear Systems

In this section we initially provide a definition for the Cholesky factorization. Next, we present a blocked algorithm that will be later evolved in order to obtain a procedure for the decomposition of hierarchical matrices.

There exist three algorithmic variants (procedures) to compute the Cholesky factorization [5]. We next review the blocked right-looking (BRL) procedure, which underlies the scheme to compute the ℋ-Cholesky factorization detailed later in this section.

Let us assume that the n × n SPD matrix A is partitioned into n_t × n_t blocks of dimension t_s × t_s each where, for simplicity, we consider that n = n_tt_s. This consideration simplifies the description of the algorithms next. However, our codes operate on blocks on any dimension and this consideration has no effect on the performance results later in the paper. Let us denote the blocks defined by this partitioning as A_ij, with i, j ∊ {1, 2,...,n_t}. The BRL algorithm for the Cholesky factorization then commences by computing the factorization of the top-left block of matrix A:

A1,1=L1,1L1,1T.$$\begin{array}{} \displaystyle {A_{1,1}} = {L_{1,1}}L_{1,1}^T. \end{array}$$

Next, it solves n_t − 1 lower triangular systems (with t_s right-hand side vectors each) to obtain the first column of blocks of L:

Li,1 := Ai,1L1,1−T,i=2,3,...,nt;$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{L_{i,1}}{\rm{ : = }}{A_{i,1}}L_{1,1}^{ - T},}&{i = 2,3,...,{n_t};} \end{array} \end{array}$$

and then it performs the update of the trailing sub-matrix:

Ai,j := Ai,j−Li,1⋅L1,j−T,i=2,3,...,nt;j=2,3,...,i.$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{A_{i,j}}{\rm{ : = }}{A_{i,j}} - {L_{i,1}} \cdot L_{1,j}^{ - T},}&{i = 2,3,...,{n_t};} \end{array}}&{j = 2,3,...,i.} \end{array} \end{array}$$

Note that, because of the symmetry of the matrix, it is not necessary to compute the first row of blocks of L nor update those blocks in the upper triangular part of A. Therefore, these operations complete the first iteration of the BRL algorithm, and the process is then (recursively) repeated by computing the Cholesky factorization of the trailing sub-matrix A_2:nt,2:nt ; see Algorithm 1. A practical implementation of this algorithm overwrites the entries in the lower triangular part of A with those of L.

An inspection of the operations that compose the BRL procedure exposes a collection of data dependencies between the blocks of A. For example, L_2,1 cannot be computed until the factorization of A_1,1 is ready. Similarly, A_2,2 cannot be updated until L_2,1 has been computed. Figure 3 illustrates the data dependencies for the operations that compose the first iteration of the Cholesky factorization applied to a 4 × 4 (blocked) matrix. In addition to these dependencies, there exist others connecting the operations of the first iteration of loop k to those appearing in subsequent iterations of the procedure. For example, at the beginning of the second iteration A_2,2 cannot be factorized until it has been updated with respect to L_2,1 during the first iteration.

Fig. 3

The BRL algorithm for the Cholesky factorization specifies a unique ordering of the operations, different from that dictated by the two other algorithmic variants for this factorization [5]. However, this schedule can be modified, provided the data dependencies are fulfilled, while still computing the same result. For example, the solution of the n_t − 1 lower triangular systems for L_{i, 1}, i = 2, 3,...,n_t, during the first iteration of the BRL algorithm, can be obtained in any order, as these operations are independent from each other.

Algorithm 1

BRL algorithm for the Cholesky factorization.

We next present a procedure for the ℋ-Cholesky factorization that is obtained from a specialization of the BRL algorithm that takes into account the hierarchical structure of A. For this purpose, assume that A ∊ ℝ^n×n is a hierarchical SPD matrix, partitioned as in Figure 2, and consider that all blocks in that partitioning are dense. Then, the following sequence of operations O1–O10 computes the ℋ-Cholesky factorization of A:

The operations in this sequence correspond to three basic linear algebra building blocks (or computational kernels):

Cholesky factorization: O1, O4 O7, O10;

These kernels are also present in the BRL Cholesky factorization procedure introduced in Algorithm 1. The only difference is that some of the kernels in the ℋ-Cholesky factorization involve operands of different sizes.

Figure 4 provides a graphical representation of the ℋ-Cholesky factorization process and the corresponding data dependencies between blocks. The five plots in the figure are to be read from left to right, starting with those in the top row first. They respectively correspond to the groups of operations:

Fig. 4

that appeared separated with horizontal lines in the sequence of operations for the ℋ-Cholesky factorization listed above. In the figure, the top-left matrix specifies the sequence of dependencies O1→O2→O3; and the top-right matrix corresponds to the implicit dependency O5→O6. For simplicity, some of the data dependencies are not shown, but they can be easily derived from the context.

We note that if any of the matrix blocks of A is compressed in low-rank factorized form, as is natural in for an ℋ-matrix, the operations involving this block will have to be performed in the appropriate ℋ-arithmetic. In any case, the dependencies remain the same. Conversely, in case a matrix block only contains zeros, some of the operations may become unnecessary. For example, if A_2,1 is “null” (that is A_2,1 ≡ 0), then L_2,1 ≡ 0, and O2, O3 do not need to be computed. Furthermore, in that case O1 and O4 can be computed in parallel as there exist no dependency connecting them.

The BRL algorithm decomposes the Cholesky factorization into a collection of four types of basic computational kernels operating on the matrix blocks: Cholesky factorization (of a diagonal block A_kk), triangular solve (line 4), symmetric rank-t_s update (line 8, when i = j), and a general matrix-matrix multiplication (line 8, when i ≠ j). Efficient implementations of these four building blocks are provided in vendor as well as open-source libraries for dense linear algebra such as Intel MKL, NVIDIA cuBLAS, IBM ESSL, GotoBLAS, OpenBLAS and BLIS. Most of these implementations are also multi-threaded yielding a direct parallel execution on a multicore architecture by simply linking any of them to the application that invokes these kernels. In particular, this is the case of an implementation of the BRL procedure (or its hierarchical generalization) that relies on a multi-threaded instance of the BLAS to perform these computations.

This first approach, based on a multi-threaded BLAS and transparent to the programmer (as it does not require any change in the original code), will only deliver a substantial fraction of the machine peak performance in case the matrix blocks involved in the operations performed inside the BLAS are sufficiently large with respect to the number of cores of the platform. Unfortunately, for hierarchical matrices, many of the blocks exhibit a reduced dimension, constraining the parallel efficiency of this simple approach.

The BRL algorithm consists of three nested loops, indexed by variables k, j, i. Among them, it was already argued that the triangular solves inside the loop in lines 3–5 are independent of each other; and the same property applies to the operations comprised by nested loops in lines 6–10. This implies that a direct approach to extract loop-parallelism from this code can employ OpenMP to parallelize these two cases using a parallel for directive. When applied to the first case, for example, this means that the iterations of loop j are distributed among the threads in charge of executing the code in parallel. The type of distribution can be adjusted via specific clauses. However, given that all iterations of the loop have the same theoretical cost (and can be expected to contribute a similar practical cost), a simple static schedule will be, in general, appropriate.

The parallelization of the second case, corresponding to the loops in lines 6–10, is slightly more complex because of the nested organization of the loops. Here, these loops cannot be tackled by using a collapse clause that internally transforms them into a single loop, because of the dependency between the counter of the innermost loop on the outermost one (index j iterates from k + 1 till i.) A solution is to parallelize the outermost loop only, but tune this with an appropriate scheduling clause that enforces a fair distribution of the workload among the threads.

In summary, this second approach provides an easy-to-program strategy to parallelize the BRL procedure (as well as the corresponding generalization to hierarchical matrices). However, the performance will depend on the number of iterations of the loops and may be constrained because, at execution, each loop parallelized in this manner implicitly requires a synchronization point for the threads upon completion.

The third parallelization approach is the most sophisticated strategy but also the option that can be expected to render higher performance. The idea is to rely on the “sequential” BRL procedure in order to identify the real data dependencies between the operations of the algorithm. Each one of these operations becomes then a task/node and the dependencies act as vertices connecting pairs of nodes of a task dependency graph (TDG). This graph defines a partial order that specifies multiple correct orderings for the operations. Concretely, it depends only on the actual dependencies of the mathematical operation but not on the algorithm itself. For example, when applied to all three known variants for the computation of the Cholesky factorization, the result is the same TDG. In other words, the execution of the operations comprised by the Cholesky factorization, executed in any order that fulfills the dependencies registered in the TDG, will produce the same correct solution as a sequential run (except for small variations due to the use of finite-precision arithmetic and rounding error).

This third approach aims to maximize the degree of concurrency exposed to the hardware, which offers multiple (correct) alternatives in the sequencing of the operations involved by the factorization. In principle, exploiting this type of parallelization approach may seem complex, because it requires to encode the algorithm, discover the data dependencies between operations (preferably, at execution time), and keep track of the order in which they are issued for execution. Fortunately, this burden can be greatly alleviated by employing a proper parallel programming model, assisted by a runtime with support for task-parallelism. This is the case of OmpSs, which was originally conceived to exploit this type of parallelism, and also for releases 3.0 and higher of OpenMP. The mechanism to automatically detect dependencies, the implications from the point of view of decomposing tasks into finer-grain sub-operations, and the consequences of the data storage layout are all discussed in [2], which was focussed on the implementations details.