Rectangular Full Packed Format for Cholesky’s Algorithm: Factorization, Solution and Inversion

A very important class of linear algebra problems deals with a coefficient matrix $A$ that is symmetric and positive definite [Dongarra et al. (1998), Demmel (1997), Golub and Van Loan (1996), Trefethen and Bau (1997)]. Because of symmetry it is only necessary to store either the upper or lower triangular part of the matrix $A$.

We describe Rectangular Full Packed Format (RFPF). It transforms a standard Packed Array AP of size $NT=N(N+1)/2$ to a full 2D array. This means that performance of LAPACK’s [Anderson et al. (1999)] packed format routines becomes equal to or better than their full array counterparts. RFPF is a variant of Hybrid Full Packed (HFP) format [Gunnels and Gustavson (2004)]. RFPF is a rearrangement of a Standard full format rectangular Array SA of size LDA*N where $LDA\ge N$. Array SA holds a triangular part of a symmetric, triangular, or Hermitian matrix $A$ of order $N$. The rearrangement of array SA is equal to compact full format Rectangular Array AR of size $LDA1\ast N1=NT$ and hence array AR like array AP uses minimal storage. (The specific values of LDA1 and N1 can vary depending on various cases and they will be specified later during the text.) Array AR will hold a full rectangular matrix $A_R$ obtained from a triangle of matrix $A$. Note also that the transpose of the rectangular matrix $A_R^T$ resides in the transpose of array AR and hence also represents $A$. Therefore, Level 3 BLAS [Dongarra et al. (1990b), Dongarra et al. (1990a)] can be used on array AR or its transpose. In fact, with the equivalent LAPACK algorithm which uses the array AR or its transpose, the performance is slightly better than standard LAPACK algorithm which uses the array SA or its transpose. Therefore, this offers the possibility to replace all packed or full LAPACK routines with equivalent LAPACK routines that work on array AR or its transpose. For examples of transformations of a matrix $A$ to a matrix $A_R$ see the figures in Section 6.

RFPF is closely related to HFP format, see [Gunnels and Gustavson (2004)], which represents $A$ as the concatenation of two standard full arrays whose total size is also $NT$. A basic simple idea leads to both formats. Let $A$ be an order $N$ symmetric matrix. Break $A$ into a block 2–by–2 form

When $N=2k$ is even, the lower triangle of $A_{11}$ and the upper triangle of $A_{22}^T$ can be concatenated together along their main diagonals into a $(k+1)$–by–$k$ dense matrix (see the figures where N is even in Section 6). This last operation is the crux of the basic simple idea. The off-diagonal block $A_{21}$ is $k$–by–$k$, and so it can be appended below the $(k+1)$–by–$k$ dense matrix. Thus, the lower triangle of $A$ can be stored as a single $(N+1)$–by–$k$ dense matrix $A_R$. In effect, each block matrix $A_{11}$, $A_{21}$ and $A_{22}$ is now stored in ‘full format’. This means all entries of matrix $A_R$ in array AR of size $LDA1=N+1$ by $N1=k$ can be accessed with constant row and column strides. So, the full power of LAPACK’s block Level 3 codes are now available for RFPF which uses the minimum amount of storage. Finally, matrix $A_R^T$ which has size $k$–by–$(N+1)$ is represented in the transpose of array AR and hence has the same desirable properties. There are eight representations of RFPF. The matrix $A$ can have have either odd or even order $N$, or it can be represented either in standard lower or upper format or it can be represented by either matrix $A_R$ or its transpose $A_R^T$ giving $2^3=8$ representations in all.

All eight cases or representations are presented in Section 6. The RFPF matrices are in the upper right part of the figures. We have introduced colors and horizontal lines to try to visually delineate triangles $T_1$, $T_2$ representing lower, upper triangles of symmetric matrices $A_{11}$, $A_{22}^T$ respectively and square or near square $S_1$ representing matrices $A_{21}$. For an upper triangle of $A$, $T_1$, $T_2$ represents lower, upper triangles of symmetric matrices $A_{11}^T$, $A_{22}$ respectively and square or near square $S_1$ representing matrices $A_{12}$. For both lower and upper triangles of $A$ we have, after each $a_{i,j}$, added its position location in the arrays holding matrices $A$ and $A_R$.

We now consider performance aspects of using RFPF in the context of using LAPACK routines on triangular matrices stored in RFPF. Let $X$ be a Level 3 LAPACK routine that operates either on full format. $X$ has a full Level 3 LAPACK block 2–by–2 algorithm, call it $FX$. We write a simple related partition algorithm (SRPA) with partition sizes $n1$ and $n2$ where $n1+n2=N$. Apply the new SRPA using the new RFPF. The new SRPA almost always has four major steps consisting entirely of calls to existing full format LAPACK routines in two steps and calls to Level 3 BLAS in the remaining two steps, see Figure 3.

Section 6 shows $FX$ algorithms equal to factorization, solution and inversion algorithms on symmetric positive definite or Hermitian matrices.

The Cholesky factorization of a symmetric and positive definite matrix $A$ can be expressed as

where $L$ and $U$ are lower triangular and upper triangular matrices.

Break the matrices $L$ and $U$ into 2–by–2 block form in the same way as was done for the matrix $A$ in Equation (1):

We now have

Using Equations (2) and equating the blocks of Equations (1) and Equations (4) gives us the basis of a 2–by–2 block algorithm for Cholesky factorization using RFPF. We can now express each of these four block equalities by calls to existing LAPACK and Level 3 BLAS routines. An example, see Section 6, of this is the three block equations is $L_{11}{L}_{11}^T=A_{11}$, $L_{21}{L}_{11}^T=A_{21}$ and $L_{21}{L}_{21}^T+L_{22}{L}_{22}=A_{22}$. The first and second of these block equations are handled by calling LAPACK’s POTRF routine $L_{11}\leftarrow A_{11}$ and by calling Level 3 BLAS TRSM via $L_{21}\leftarrow L_{21}{L}_{11}^{-T}$. In both these block equations the Fortran equality of replacement ($\leftarrow$) is being used so that the lower triangle of $A_{11}$ is being replaced $L_{11}$ and the nearly square matrix $A_{21}$ is being replaced by $L_{21}$. The third block equation breaks into two parts: $A_{22}\leftarrow L_{21}{L}_{21}^T$ and $L_{22}\leftarrow A_{22}$ which are handled by calling Level 3 BLAS SYRK or HERK and by calling LAPACK’s POTRF routine. At this point we can use the flexibility of the LAPACK library. In RFPF $A_{22}$ is in upper format (upper triangle) while in standard format $A_{22}$ is in lower format (lower triangle). Due to symmetry, both formats of $A_{22}$ contain equal values. This flexibility allows LAPACK to accommodate both formats. Hence, in the calls to SYRK or HERK and POTRF we set uplo = ’U’ even though the rectangular matrix of SYRK and HERK comes from a lower triangular formulation.

New LAPACK like routine PFTRF performs these four computations. PF was chosen to fit with LAPACK’s use of PO and PP. The PFTRF routine covers the Cholesky Factorization algorithm for the eight cases of the RFPF. Section 6 has Figure 4 with four subfigures. Here we are interested in the first and second subfigure. The first subfigure contains the layouts of matrices $A$ and $A_R$. The second subfigure has the Cholesky factorization algorithm obtained by simple algebraic manipulations of the three block equalities obtained above.

In Section 3 we obtained the 2–by–2 Cholesky factorization (3) of matrix $A$. Now, we can solve the equation $AX=B$:

• If $A$ has lower triangular format then

  $$\begin{array}{c}LY=B\text{and}{L}^{T}X=Y\text{(in the symmetric case)}\\ LY=B\text{and}{L}^{H}X=Y\text{(in the Hermitian case)}\end{array}$$

  (5)

• If $A$ has an upper triangular format then

  $$\begin{array}{c}{U}^{T}Y=B\text{and}UX=Y\text{(in the symmetric case)}\\ U^{H}Y=B\text{and}UX=Y\text{(in the Hermitian case)}\end{array}$$

  (6)

$B$, $X$ and $Y$ are either vectors or rectangular matrices. $B$ contains the RHS values. $X$ and $Y$ contain the solution values. $B$, $X$ and $Y$ are vectors when there is one RHS and matrices when there are many RHS. The values of $X$ and $Y$ are stored over the values of $B$.

Expanding (5) and (6) using (3) gives the forward substitution equations

and the back substitution equations

The Equations (7) and (8) gives the basis of a $2\times 2$ block algorithm for Cholesky solution using RFPF format. We can now express these two sets of two block equalities by using existing Level 3 BLAS routines. An example, see Section 6, of the first set of these two block equalities is $L_{11}{Y}_1=B_1$ and $L_{21}{Y}_1+L_{22}{Y}_2=B_2$. The first block equality is handled by Level 3 BLAS TRSM: $Y_1\leftarrow L_{11}^{-1}{B}_1$. The second block equality is handled by Level 3 BLAS GEMM and TRSM: $B_2\leftarrow B_2-L_{21}{Y}_1$ and $Y_2\leftarrow L_{22}^{-1}{Y}_2$. The backsolution routines are similarly derived. One gets $X_2\leftarrow L_{22}^{-T}{Y}_2$, $Y_1\leftarrow Y_1-L_{21}^T{X}_2$ and $X_1\leftarrow L_{11}^{-T}{Y}_1$.

New LAPACK like routine PFTRS performs these two solution computations for the eight cases of RFPF. PFTRS calls a new Level 3 BLAS TFSM in the same way that POTRS calls TRSM. The third subfigure in Section 6 gives the Cholesky solution algorithm using RFPF obtained by simple algebraic manipulation of the block Equations (7) and (8).

Following LAPACK we consider the following three stage procedure:

1. Factorize the matrix $A$ and overwrite $A$ with either $L$ or $U$ by calling PFTRF; see Section 3.

2. Compute the inverse of either $L$ or $U$. Call these matrices $W$ or $V$ and overwrite either $L$ or $U$ with them. This is done by calling new routine new LAPACK like TFTRI.

3. Calculate either the product $W^{T}W$ or $V{V}^T$ and overwrite either $W$ or $V$ with them.

As in Sections 3 and 4 we examine 2–by–2 block algorithms for the steps two and three above. In Section 3 we obtain either matrices $L$ or $U$ in RFPF. Like LAPACK inversion algorithms for POTRI and PPTRI, this is our starting point for our LAPACK inversion algorithm using RFPF. The LAPACK inversion algorithms for POTRI and PPTRI also follow from steps two and three above by first calling in the full case LAPACK TRTRI and then calling LAPACK LAUUM.

Take the inverse of Equation (2) and obtain

where $W$ and $V$ are lower and upper triangular matrices.

Using the 2–by–2 blocking for either $L$ or $U$ in Equation (3) we obtain the following 2–by–2 blocking for $W$ and $V$:

From the identities $WL=LW=I$ and $VU=UV=I$ and the 2–by–2 block layouts of Equations (3) and 3), we obtain three block equations for $W$ and $V$ which can be solved using LAPACK routines for TRTRI and Level 3 BLAS TRMM. An example, see Figure 4, of these three block equations is $L_{11}{W}_{11}=I$, $L_{21}{W}_{11}+L_{22}{W}_{21}=0$ and $L_{22}{W}_{22}=I$. The first and third of these block equations are handled by LAPACK TRTRI routines as $W_{11}\leftarrow L_{11}^{-1}$ and $V_{22}\leftarrow U_{22}^{-1}$. In the second inverse computation we use the fact that $L_{22}$ is equally represented by it transpose $L_{22}^T$ which is $U_{22}$ in RFPF. The second block equation leads to two calls to Level 3 BLAS TRMM via $L_{21}\leftarrow -L_{21}{W}_{11}$ and $W_{21}=W_{22}{L}_{21}$. In the last two block equations the Fortran equality of replacement ($\leftarrow$) is being used so that $W_{21}=-W_{22}{L}_{21}{W}_{11}$ is replacing $L_{21}$.

Now we turn to part three of the three stage LAPACK procedure above. For this we use the 2–by–2 blocks layouts of Equation (10) and the matrix multiplications indicated by following block Equations (11) giving

where $W_{11}$, $W_{22}$, $V_{11}$, and $V_{22}$ are lower and upper triangular submatrices, and $W_{21}$ and $V_{12}$ are square or almost square submatrices. The values of the indicated block multiplications of $W$ or $V$ in Equation (11) are stored over the block values of $W$ or $V$.

Performing the indicated 2–by–2 block multiplications of Equation (11) leads to three block matrix computations. An example, see Section 6, of these three block computations is $W_{11}^T{W}_{11}+W_{21}^T{W}_{21}$, $W_{22}^T{W}_{21}$ and $W_{22}^T{W}_{22}$. Additionally, we want to overwrite the values of these block multiplications on their original block operands. Block operand $W_{11}$ only occurs in the (1,1) block operand computation and hence can be overwritten by a call to LAPACK LAUUM, $W_{11}\leftarrow W_{11}^T{W}_{11}$, followed by a call to Level 3 BLAS SYRK or HERK, $W_{11}\leftarrow W_{11}+W_{21}^T{W}_{21}$. Block operand $W_{21}$ now only occurs in the (2,1) block computation and hence can be overwritten by a call to Level 3 BLAS TRMM, $W_{21}\leftarrow W_{22}^T{W}_{21}$. Finally, block operand $W_{22}$ can be overwritten by a call to LAPACK LAUUM, $W_{22}\leftarrow W_{22}^T{W}_{22}$.

The fourth subfigure in Section 6 has the Cholesky inversion algorithms using RFPF based on the results of this Section. New LAPACK routine, PFTRI, performs this computation for the eight cases of RFPF.

This section contains three figures.

1. The first figure describes the RFPF (Rectangular Full Packed Format) and gives algorithms for Cholesky factorization, solution and inversion of symmetric positive definite matrices, where $N$ is odd, uplo = ’lower’, and trans = ’no transpose’. This figure has four subfigures.

   1. The first subfigure depicts the lower triangle of a symmetric positive definite matrix $A$ in standard full and its representation by the matrix $A_R$ in RFPF.

   2. The second subfigure gives the RFPF Cholesky factorization algorithm and its calling sequences of the LAPACK and BLAS subroutines, see Section 3.

   3. The third subfigure gives the RFPF Cholesky solution algorithm and its calling sequences to the LAPACK and BLAS subroutines, see Section 4.

   4. The fourth subfigure in each figure gives the RFPF Cholesky inversion algorithm and its calling sequences to the LAPACK and BLAS subroutines, see Section 5.

2. The second figure shows the transformation from full to RFPF of all “no transform” cases.

3. The third figure depicts all eight cases in RFPF.

The data format for $A$ has $LDA=N$. Matrix $A_R$ has $LDAR=N$ if $N$ is odd and $LDAR=N+1$ if $N$ is even and $n1$ columns where $n1=\lceil N/2\rceil$. Hence, matrix $A_R$ always has LDAR rows and $n1$ columns. Matrix $A_R^T$ always has $n1$ rows and LDAR columns and its leading dimension is equal to $n1$. Matrix $A_R$ always has $LDAR\times n1=NT=N(N+1)/2$ elements as does matrix $A_R^T$.

The order $N$ of matrix $A$ in the first figure is seven and six or seven in the remaining two figures.

The RFPF Cholesky factorization (Section 3), Cholesky solution (Section 4), and Cholesky inversion (Section 5) algorithms are equivalent to the traditional algorithms in the books [Dongarra et al. (1998), Demmel (1997), Golub and Van Loan (1996), Trefethen and Bau (1997)]. The whole theory of the traditional Cholesky factorization, solution, inversion and BLAS algorithms carries over to this three Cholesky and BLAS algorithms described in Sections 3, 4, and 5. The error analysis and stability of these algorithms is very well described in the book of [Higham (1996)]. The difference between LAPACK algorithms PO, PP and RFPF³³3full, packed and rectangular full packed. is how inner products are accumulated. In each case a different order is used. They are all mathematically equivalent, and, stability analysis shows that any summation order is stable.

The LAPACK library [Anderson et al. (1999)] routines POTRF/PPTRF, POTRI/PPTRI, and POTRS/PPTRS are compared with the RFPF routines PFTRF, PFTRI, and PFTRS for Cholesky factorization (PxTRF), Cholesky inverse (PxTRI) and Cholesky solution (PxTRS) respectively. In the previous sentence, the character ’x’ can be ’O’ (full format), ’P’ (packed format), or ’F’ (RFPF). In all cases long real precision arithmetic (also called double precision) is used. Sometimes we also show results for long complex precision (also called complex*16). Results were obtained on several different computers using everywhere the vendor Level 3 and Level 2 BLAS. The sequential performance results were done on the following computers:

$\bullet$ Sun Fire E25K (newton):

72 UltraSPARC IV+ dual-core CPUs (1800 MHz/ 2 MB shared L2-cache, 32 MB shared L3-cache), 416 GB memory (120 CPUs/368 GB). Further information at “http://www.gbar.dtu.dk/index.php/Hardware”.

$\bullet$ SGI Altix 3700 (Freke):

64 CPUs - Intel Itanium2 1.5 GHz/6 MB L3-cache. 256 GB memory. Peak performance: 384 GFlops. Further information at “http://www.cscaa.dk/freke/”.

$\bullet$ Intel Tigerton computer (zoot):

quad-socket quad-core Intel Tigerton 2.4GHz (16 total cores) with 32 GB of memory. We use Intel MKL 10.0.1.014.

$\bullet$ DMI Itanium:

CPU Intel Itanium2: 1.3 GHz, cache: 3 MB on-chip L3 cache.

$\bullet$ DMI NEC SX-6 computer:

8 CPU’s, per CPU peak: 8 Gflops, per node peak: 64 Gflops, vector register length: 256.

The performance results are given in Figures 7 to 15. In Appendix A, we give the table data used in the figures, see Tables 2 to 27. We also give speedup numbers, see Tables 28 to 35.

The figures from 7 to 10 are paired. Figure 7 (double precision) and Figure 8 (double complex precision) present results for the Sun UltraSPARC IV+ computer. Figure 9 (double precision) and Figure 10 (double complex precision) present results for the SGI Altix 3700 computer. Figure 11 (double precision) presents results for the Intel Itanium2 computer. Figure 12 (double precision) presents results for the NEC SX-6 computer. Figure 13 (double precision) presents results for the quad-socket quad-core Intel Tigerton computer using reference LAPACK-3.2.0 (from netlib.org). Figure 14 (double precision) presents results for the quad-socket quad-core Intel Tigerton computer using vendor LAPACK library (MKL-10.0.1.14).

Figure 15 shows the SMP parallelism of these subroutines on the IBM Power4 (clock rate: 1300 MHz; two CPUs per chip; L1 cache: 128 KB (64 KB per CPU) instruction, 64 KB 2-way (32 KB per CPU) data; L2 cache: 1.5 MB 8-way shared between the two CPUs; L3 cache: 32 MB 8-way shared (off-chip); TLB: 1024 entries) and SUN UltraSPARC-IV (clock rate: 1350 MHz; L1 cache: 64 kB 4-way data, 32 kB 4-way instruction, and 2 kB Write, 2 kB Prefetch; L2 cache: 8 MB; TLB: 1040 entries) computers respectively. They compare SMP times of PFTRF, vendor POTRF and reference PPTRF.

The RFPF packed results greatly outperform the packed and more often than not are better than the full results. Note that our timings do not include the cost of sorting any LAPACK data formats to RFPF data formats and vice versa. We think that users will input their matrix data using RFPF. Hence, this is our rationale for not including the data transformation times.

For all our experiments, we use vendor Level 3 and Level 2 BLAS. For all our experiments except Figure 13 and Figure 15, we use the provided vendor library for LAPACK and BLAS.

We include comparisons with reference LAPACK for the quad-socket quad-core Intel Tigerton machine in Figure 13. In this case, the vendor LAPACK library packed storage routines significantly outperform the LAPACK reference implementation from netlib. In Figure 14, you find the same experiments on the same machine but, this time, using the vendor library (MKL-10.0.1.014). We think that MKL is using the reference implementation for Inverse Cholesky (packed and full format). For Cholesky factorization, we see that both packed and full format routines (PPTRF and POTRF) are tuned. But even, in this case, our RFPF storage format results are better.

When we compare RFPF with full storage, results are mixed. However, both codes are rarely far apart. Most of the performance ratios are between 0.95 to 1.05 overall. But, note that the RFPF performance is more uniform over its versions (four presented; the other four are for n odd ). For LAPACK full (two versions ), the performance variation is greater. Moreover, in the case of the inversion on quad-socket quad-core Tigerton (Figure 13 and Figure 14) RFPF clearly outperforms both variants of the full format.

As mentioned in the introduction, as of release 3.2 (November 2008), LAPACK supports a preliminary version of RFPF. Ultimately, the goal would be for RFPF to support as many functionnalities as full format or standard packed format does. The 44 routines included in release 3.2 for RFPF are given in Table 1. The names for the RFPF routines follow the naming nomenclature used by LAPACK. We have added the format description letters: PF for Symmetric/Hermitian Positive Definite RFPF (PO for full, PP for packed), SF for Symmetric RFPF (SY for full, SP for packed), HF for Hermitian RFPF (HE for full, HP for packed), and TF for Triangular RFPF (TR for full, TP for packed).

Currently, for the complex case, we assume that the transpose complex-conjugate part is stored whenever the transpose part is stored in the real case. This corresponds to the theory developed in this present manuscript. In the future, we will want to have the flexibility to store the transpose part (as opposed to transpose complex conjugate) whenever the transpose part is stored in the real case. In particular, this feature will be useful for complex symmetric matrices.

This paper describes RFPF as a standard minimal full format for representing both symmetric and triangular matrices. Hence, from a user point of view, these matrix layouts are a replacement for both the standard formats of DLA, namely full and packed storage. These new layouts possess three good features: they are efficient, they are supported by Level 3 BLAS and LAPACK full format routines, and they require minimal storage.

The results in this paper were obtained on seven computers, an IBM, a SGI, two SUNs, Itanium, NEC, and Intel Tigerton computers. The IBM machine belongs to the Center for Scientific Computing at Aarhus, the SUN machines to the Danish Technical University, the Itanium and NEC machines to the Danish Meteorological Institute, and the Intel Tigerton machine to the Innovative Computing Laboratory at the University of Tennessee.

We would like to thank Bernd Dammann for consulting on the SUN systems; Niels Carl W. Hansen for consulting on the IBM and SGI systems; and Bjarne Stig Andersen for obtaining the results on the Itanium and NEC computers. We thank IBMers John Gunnels who worked earlier on the HFPF format and JP Fasano who was instrumental in getting the source code released by the IBM Open Source Committee. We thank Allan Backer for discussions about an older version of this manuscript.