A recursive formulation of Cholesky factorization of a matrix in packed storage.
A new compact way to store a symmetric or triangular matrix called RPF for Recursive Packed Format is fully described. Novel ways to transform RPF to and from standard packed format are included. A new algorithm, called RPC for Recursive Packed Cholesky, that operates on the RPG format is presented. ALgorithm RPC is basd on level-3 BLAS and requires variants of algorithms TRSM and SYRK that work on RPF. We call these RP_TRSM and RP_SYRK and find that they do most of their work by calling GEMM. It follows that most of the execution time of RPC lies in GEMM. The advantage of this storage scheme compared to traditional packed and full storage is demonstrated. First, the RPC storage format uses the minimal    amount of storage for the symmetric or triangular matrix. Second, RPC gives a level-3 implementation of Cholesky factorization whereas standard packed implementations are only level 2. Hence, the performance of our RPC implementation is decidedly superior. Third, unlike fixed block size algorithms, RPC, requires no block size tuning parameter. We present performance measurements on several current architectures that demonstrate improvements over the traditional packed routines. Also MSP parallel computations on the IBM SMP computer are made. The graphs that are attached in Section 7 show that the RPC algorithms are superior by a factor between 1.6 and 7.4 for order around 1000, and between 1.9 and 10.3 for order around 3000 over the traditional packed algorithms. For some    architectures, the RPC performance results are almost the same or even better than the traditional full-storage algorithms results.
Introduction
A very important class of linear algebra problems are those in which the coefcient
matrix A is symmetric and positive denite [5, 11, 23]. Because of the
symmetry it is only necessary to store either the upper or lower triangular part
of the matrix.
Lower triangular caseB
7 14 21 28
Upper triangular caseB
43
9

Figure

1: The mapping of 7 7 matrix for the LAPACK Cholesky Algorithm
using the full storage (LDA= 7 if in Fortran77)
Lower triangular caseB
3 9 14
Upper triangular caseB

Figure

2: The mapping of 7  7 matrix for the LAPACK Cholesky Algorithm
using the packed storage
1.1 LAPACK POTRF and PPTRF subroutines
The LAPACK library[3] oers two dierent kind of subroutines to solve the same
problems, for instance POTRF 1 and PPTRF both factorize symmetric, positive
Four names SPOTRF, DPOTRF, CPOTRF and ZPOTRF are used in LAPACK for real
symmetric and complex Hermitian matrices[3], where the rst character indicates the precision
and arithmetic versions: S { single precision, D { double precision, C { complex and Z { double
complex. LAPACK95 uses one name LA POTRF for all versions[7]. POTRF and/or PPTRF
express, in this paper, any precision, any arithmetic and any language version of the PO
and/or PP matrix factorization algorithms.
denite matrices by means of the Cholesky algorithm. The only dierence is
the way the triangular matrix is stored (see gures 1 and 2).
In the POTRF case the matrix is stored in one of the lower left or upper right
triangles of a full square matrix[16, page 64] 2 , the other triangle is wasted (see
gure 1). Because of the uniform storage scheme, blocking and level 3 BLAS[8]
subroutines can be employed, resulting in a high speed solution.
In the PPTRF case the matrix is kept in packed storage ([1], [16, page 74,
75]), which means that the columns of the lower or upper triangle are stored
consecutively in a one dimensional array (see gure 2). Now the triangular
matrix only occupies the strictly necessary storage space but the nonuniform
storage scheme means that use of full storage BLAS is impossible and only the
level 2 BLAS[20, 9] packed subroutines can be employed, resulting in a low
speed solution.
To summarize, there is a choice between high speed with waste of memory
versus low speed with no waste of memory.
1.2 A new Way of Storing Real Symmetric and Complex
Hermitian and, in either case, Positive Denite Matrice

Together with some new recursively formulated linear algebra subroutines, we
propose a new way of storing a lower or upper triangular matrix that solves this
dilemma[14, 24]. In other words we obtain the speed of POTRF with the amount
of memory used by PPTRF. The new storage scheme is named RPF, recursive
packed format(see gure 4), and it is explained below.
The benet of recursive formulations of the Cholesky factorization and the
LU decomposition is described in the works of Gustavson [14] and Toledo [22].
The symmetric, positive denite matrix in the Cholesky case is kept in full
matrix storage, and the emphasis in these works are the better data locality
and thus better utilization of the computers memory hierarchy, that recursive
formulations oer. However, the recursive packed formulation also has this
property.
We will provide a very short introduction on a computer memory hierarchy
and the Basic Linear Algebra Subprograms (BLAS) before describing the
Recursive Packed Cholesky (RPC) and the Recursive Packed Format (RPF).
1.3 The Rationale behind introducing our New Recursive
Algorithm, RPC and the New Recursive Data Format,
RPF
Computers have several levels of memory. The
ow of data from the memory
to the computational units is the most important factor governing performance
of engineering and scientic computations. The object is to keep the functional
units running at their peak capacity. Through the use of a memory hierarchy
In Fortran column major, in C row major.
system (see gure 3), high performance can be achieved by using locality of
reference within a program. In the present context this is called blocking.
Registers
Cache Level 1
Cache Level 2
Cache Level 3
Shared Memory
Distributed Memory
Secondary Storage 1
Secondary Storage 2
Faster,
smaller,
more
expensive
Slower,
larger,
cheaper

Figure

3: A computer memory hierarchy
At the top of the hierarchy is a Central Processing Unit (CPU). It communicates
directly with the registers. The number of the registers is usually very
small. A Level 1 cache is directly connected to the registers. The computer
will run with almost peak performance if we are able to deliver the data to the
(level 1) cache in such way that the CPU is permanently busy. There are
several books describing problems associated with the computer memory hier-
archy. The literature in [10, 5, 11] is adequate for Numerical Linear Algebra
specialists.
The memories near the CPU (registers and caches) have a faster access to
CPU than the memories further away. The fast memories are very expensive
and this is one of the reason that they are small. The register set is tiny. Cache
memories are much larger than the set of registers. However, L1 cache is still
not large enough for solving scientic problems. Even a subproblem like matrix
factorization does not t into cache if the order of the matrix is large.
A special set of Basic Linear Algebra Subprograms (BLAS) have been developed
to address the computer memory hierarchy problem in the area of Numerical
Linear Algebra. The BLAS are documented in [20, 9, 8, 6]. BLAS are
very well summarized and explained for Numerical Linear Algebra specialists
in [10, 5].
There are three levels of BLAS: Level 1 BLAS shows vector vector opera-
tions, Level 2 BLAS shows vector matrix (and/or matrix vector) operations,
and Level 3 BLAS shows matrix matrix operations.
For Cholesky factorization one can make the following three observations
with respect to the BLAS.
1. Level 3 implementations using full storage format run fast.
2. Level 3 implementations using packed storage format rarely exist. A level 3
implementation was previously used in [16], however, at great programming
cost. Conventional Level 2 implementations using packed storage
format run, for large problem sizes, considerably slower than the full storage
implementations.
3. Transforming conventional packed storage to RPF and using our RPC
algorithm produces a Level 3 implementation using the same amount of
storage as packed storage.
1.4 Overview of the Paper
Section 2 describes the new packed storage data format and the data transformations
to and from conventional packed storage. Section 2.1 describes conventional
lower and upper triangular packed storage. Section 2.2 discusses how
to transform in place either a lower or upper trapezoid packed data format to
recursive packed data format and vice versa. Section 2.3 describes the possibility
to transpose the matrix while it is reordered from packed to recursive
packed format and vice versa. Finally, in Section 2.4 the recursive aspects of
the data transformation is described. These four subsections describe the in
place transformation pictorially via several gures.
In Sections 3.1 and 3.2, recursive TRSM and SYRK, both which work on
RPF, are described. Both routines do almost all their required
oating point
operations by calling level 3 BLAS GEMM. Finally, in Section 3.3, the RPC
algorithm is described in terms of using the recursive algorithms of Sections 3.1
and 3.2. As in Section 2, all three algorithms are described pictorially via several
gures. Note that the RPC algorithm only uses one Level 3 BLAS subroutine,
namely GEMM. Usually the GEMM routine is very specialized, highly tuned
and done by the computer manufacturer. If not, the ATLAS[25] GEMM can be
used.
Section 4 explains that the RPC algorithm is numerically stable.
Section 5 describes performance graphs of the packed storage LAPACK[3]
algorithms and of our recursive packed algorithms on several computers; the
most typical computers like COMPACQ, HP, IBM SP, IBM SMP, INTEL Pen-
tium, SGI and SUN were considered (gures results
show that the recursive packed Cholesky factorization (RP PPTRF) and the
solution (RP PPTRS) are 4 { 9 times faster than the traditional packed sub-
routines. There are three more graphs. One demonstrates successful use of
OpenMP[17, 18] parallelizing directives (gure 18). The second graph shows
that the recursive data format is also eective for the complex arithmetic (g-
ure 19). The third one shows the performance of all three algorithms for the
factorization (POTRF, PPTRF and RP PPTRF) and the solution
(POTRS, PPTRS and RP PPTRS) (gure 20).
Section 6 discusses the most important developments in this paper.
2 The recursive packed storage
A new way to store triangular matrices in packed storage called recursive packed
is presented. This is a storage scheme by its own right, and a way to explain it,
is to describe the conversion from packed to recursive packed storage and vice
versa (see gures 2 and 4).
Lower triangular caseB
22
Upper triangular caseB
9 12 15
19 2022 24
26 271

Figure

4: The mapping of 7  7 matrix for the Cholesky Algorithm using the
recursive packed storage. The recursive block division is illustrated.
2.1 Lower and upper triangular packed storage
Symmetric, complex hermitian or triangular matrices may be stored in packed
storage form (see LAPACK Users' Guide [3], IBM ESSL Library manual[16,
pages 66{67] and gure 2). The columns of the triangle are stored sequentially
in a one dimensional array starting with the rst column. The mapping between
positions in full storage and in packed storage for a triangular matrix of size m
is,
A i;j
AP i+(j 1)(2m j)=2
'L'The advantage of this storage is the saving of almost half 4 the memory
compared to full storage.
3 For upper triangular and for lower triangular of A is stored.
4 At least m  (m 1)=2. This formulae is a function of LDA (leading dimension of
and m in Fortran77. The saving in Fortran77 is m
2.2 Reordering of a lower and upper trapezoid
Packed storage
packed storage memory map
buer p(p 1)=2 words
Recursive packed storage memory map

Figure

5: Reordering of the lower packed matrix. First, the last p 1 columns
of the leading triangle are copied to the buer. Then, in place, the columns of
the accentuated rectangle are assembled in the bottom space of the trapezoid.
Last, the buer is copied back to the top of the trapezoid.
It is assumed that the matrices are stored in column major order, but the
concepts in the paper are fully applicable also if the matrices are stored in row
major order. As an intermediate step to transform a packed triangular matrix
to a recursive packed matrix, the matrix is divided into two parts along
a column thus dividing the matrix in a trapezoidal and a triangular part as
shown in g. 5 and 6. The triangular part remains in packed form, the
trapezoidal part is reordered so it consists of a triangle in packed form, and
a rectangle in full storage form. The reordering demands a buer of the size
of the triangle minus the longest column. The reordering in the lower case,
g. 5, takes the following steps. First the columns of the triangular part of
the trapezoid are moved to the buer (note that the rst column is in correct
place), then the columns of the rectangular part of the trapezoid are
moved into consecutive locations and nally the buer is copied back to the
correct location in the reordered array. If p in gure 5 is chosen to bm=2c
the rectangular submatrix will be square or deviate from a square only by a
single column. The buer size is p(p 1)=2 and the addresses of the lead-
8Packed storage
packed storage memory map
buer
(m p)(m p 1)=2 words
Recursive packed storage memory map

Figure

Reordering of the upper packed matrix. First, the rst
columns of the trailing triangle are copied to the buer. Then, in place, the
columns of the accentuated rectangle are assembled in the top space of the
trapezoid. Last, the buer is copied back to the bottom of the trapezoid.
ing triangle, the rectangular submatrix and the trailing triangle are given by,
After the reordering the leading and trailing triangles are both in the same
lower or upper packed storage scheme as the original triangular matrix. The
reordering can be implemented as subroutines,
subroutine TPZ TO TR(m; p; AP)
and
subroutine TR TO TPZ (m; p; AP)
where TPZ TO TR means the reordering of the trapezoidal part from packed
format to the triangular-rectangular format just described. TR TO TPZ is the
opposite reordering.
2.3 Transposition of the rectangular part
The rectangular part of the reordered matrix are now kept in full matrix storage.
If desired, this oers an excellent opportunity to transpose the matrix while it is
transformed to recursive packed format. If the rectangular submatrix is square
the transposition can be done completely in-place. If it deviates from a square by
a column, a buer of the size of the columns is necessary to do the transposition,
for this purpose we can reuse the buer used for the reordering.
2.4 Recursive application of the reordering
The method of reordering is applied recursively to the leading and trailing triangles
which are still in packed storage, until nally the originally triangular
packed matrix is divided in rectangular submatrices of decreasing size, all in
full storage. The implementation of the complete transformation from packed
to recursive packed format, P TO RP is (compare the gures 2 and 4),
recursive subroutine P TO RP(m;AP)
if (m > 1) then
call TPZ TO TR(m; p; AP)
call P TO RP(p; AP)
call P TO RP(m
and the inverse transformation from recursive packed to packed, RP TO P is,
recursive subroutine RP TO P(m;AP)
if (m > 1) then
call RP TO P(p; AP)
call TR TO TPZ (m; p; AP)
call RP TO P(m
The examples shown here concerns the lower triangular matrix, but the upper
triangular transformation, and the transformation with transposition follows
the same pattern. The gure 7 illustrates the recursive division of small lower
and upper triangular matrices.
Figure

7: The lower and upper triangular matrices, in recursive packed storage
data format, for 20. The rectangular submatrices, shown in the gures,
are kept in full storage in column major order, in the array containing the whole
matrices.
Recursive formulation of the Cholesky algorithm
and its necessary BLAS
Two BLAS[6] operations, the triangular solver with multiple right hand sides,
TRSM 5 and the rank k update of a symmetric matrix, SYRK are needed for the
recursive Cholesky factorization and solution, RP PPTRF 6 and RP PPTRS[2].
In this section RP TRSM, RP SYRK, RP PPTRF and RP PPTRS are formulated
recursively and their use of recursive packed operands are explained.
TRSM, SYRK, PPTRF and PPTRS operate in various cases depending of the
operands and the order of the operands. In the following we only consider single
specic cases, but the deduction of the other cases follows the same guidelines.
All the computational work in the recursive BLAS routines RP TRSM and
RP SYRK (and also RP TRMM) is done by the non recursive matrix-matrix
multiply routine GEMM[19, 25]. This is a very attractive property, since GEMM
usually is or can be highly optimized on most current computer architectures.
The GEMM operation is very well documented and explained in [12, 13, 6].
The speed of our computation depends very much from the speed of a good
GEMM. Good GEMM implementations are usually developed by computer
manufacturers. The model implementation of GEMM can be obtained from
5 On naming of TRSM, SYRK, HERK and GEMM see footnote of POTRF on page 1.
6 The prex RP indicates that the subroutine belongs to the Recursive Packed library, for
example RP PPTRF is the Recursive Packed Cholesky factorization routine.
netlib [6]; it works correctly but slowly. The Innovative Computing Laboratory
at the University of Tennessee in Knoxville developed an automatic system
called which usually can produce a very fast GEMM subroutine. Another
automatic code generator scheme for GEMM was developed at Berkeley[4].
In ESSL, see [1], GEMM and all other BLAS are produced via blocking and
high performance kernel routines. For example, ESSL produces a single kernel
routine, DATB, which has the same function as the ATLAS on chip GEMM
kernel. The principles underlying the production of both kernels are similar.
The major dierence is that ESSL's GEMM code is written by hand whereas
ATLAS' GEMM code is parametrized and run over all parameter settings until
a best parameter setting is found for the particular machine.
3.1 Recursive TRSM based on non-recursive GEMM
Fig. 8 shows the splitting of the TRSM operands. The operation now consists
of the three suboperations,
Based on this splitting, the algorithm can be programmed as follows,
recursive subroutine RP TRSM (m; n;
if (n == 1) then
do
do
else
call RP TRSM (m;
call GEMM ( 0 N
call RP TRSM (m; n
3.2 Recursive SYRK based on non-recursive GEMM
Fig. 9 shows the splitting of the SYRK operands. The operation now consists
of the three suboperations,

A 21
A 11
A 22

Figure

8: The recursive splitting of the matrices in the RP TRSM operation for
the case where SIDE=Right, UPLO=Lower and TRANSA=Transpose.
Based on this splitting, the algorithm can be programmed as follows,
recursive subroutine RP SYRK (m; n;
if (m == 1) then
do
do
else
call RP SYRK (p; n;
call GEMM ( 0 N
call RP SYRK (m
3.3 Recursive PPTRF and PPTRS based on recursive
TRSM and recursive SYRK
Fig. 10 shows the splitting of the PPTRF operand. The operation now consists
of four suboperations,
M-P
C 22
M-P
bC 22
M-P

A 11
A 21
M-P

Figure

9: The recursive splitting of the matrices in the RP SYRK operation
for the case where UPLO=Lower and TRANS=No transpose.
22 RP PPTRF
Based on this splitting the algorithm can be programmed as follows,
recursive subroutine RP PPTRF (m; AP)
if (m == 1) then
else
call RP PPTRF (p; AP)
call RP TRSM (m
call RP SYRK (m
call RP PPTRF (m
The solution subroutine RP PPTRS performs consecutive triangular solutions
to the transposed and the non-transposed Cholesky factor. This routine
is not explicitly recursive, as it just calls the recursive RP TRSM twice.
4 Stability of the Recursive Algorithm
The paper [24] shows that the recursive Cholesky factorization algorithm is
equivalent to the traditional algorithms in the books[5, 11, 23]. The whole theAP
A 21
M-P
A 22
M-P
M-P

M-P

Figure

10: The recursive splitting of the matrix in the RP PPTRF operation
for the case where UPLO=Lower.
ory of the traditional Cholesky factorization and BLAS (TRSM and SYRK) algorithms
carries over to the recursive Cholesky factorization and BLAS (TRSM
and SYRK) algorithms described in Section 3. The error analysis and stability
of these algorithms is very well described in the book of Nicholas J. Higham[15].
The dierence between LAPACK algorithms PO, PP and RP 7 is how inner
products are accumulated. In each case a dierent order is used. They are all
mathematically equivalent, and, stability analysis shows that any summation
order is stable.
5 Performance results
SUN UltraSparc II @ 400 MHz
SGI R10000 @ 195 MHz
COMPAQ Alpha EV6 @ 500 MHz
HP PA-8500 @ 440 MHz
INTEL Pentium III @ 500 MHz

Table

1: Computer names
The new recursive packed BLAS (RP TRSM and RP SYRK), and the new
recursive packed Cholesky factorization and solution (RP PPTRF and RP PPTRS)
routines were compared to the traditional LAPACK subroutines, both concerning
the results and the performance. The comparisons were made on seven
7 full, packed and recursive packed.
dierent architectures, listed in the Table 1. The result graphs are attached in
the appendix of this paper. The double precision arithmetic in Fortran90[21]
was used in all cases.
IBM-PPC ESSL 3.1.0.0 -lesslsmp
IBM-PW2 ESSL 2.2.2.0 -lesslp2
SUN Sun Performance Library 2.0 -lsunperf
SGI Standard Execution Environment 7.3 -lblas
COMPAQ DXML V3.5 -ldxmp ev6
HP HP-UX PA2.0 BLAS Library 10.30 -lblas
INTEL ATLAS 3.0

Table

2: Computer library versions
The following procedure was used in carrying out our performance tests.
On each machine the recursive and the traditional routines were compiled
with the same compiler and compiler
ags and they call the same vendor
optimized, or otherwise optimized, BLAS library. The BLAS library
versions can be seen in Table 2.
The compared recursive and traditional routines received the same input
and produced the same output for each time measurement. The time
spent in reordering the matrix to and from 8 recursive packed format is
included in the run time for both RP PPTRF and RP PPTRS. For the
traditional routines there was no data transformation cost.
The CPU time is measured by the timing function ETIME except on the
PowerPC machine, which is a 4 way SMP. On this machine the run time
was measured by the wall clock time by means of a special IBM utility
function called RTC. Except for the operating system no other programs
were running during these test runs.
For each machine the timings were made for a sequence of matrix sizes
ranging from in steps of In case of the
HP and Intel machines the matrix size starts at We start at
because the resolution of the ETIME utility was too coarse. The
number of right hand sides were taken to be nrhs = n=10. Due to memory
limitations on the actual HP machine, this test series could only range to
The operation counts for Cholesky factorization and solution are
8 However it is only necessary to perform the to transformation in RP PPTRF and no
transformation in RP PPTRS, to get the correct results.
where n is the number of equations and nrhs the number of right hand
sides. These counts were used to convert run times to Flop rates.
Ten gures (gure show performance graph comparisons,
between the new RPC algorithms and the traditional LAPACK algorithms. The
RPC algorithms use the RPF data format in all comparisons. As mentioned the
cost of transforming from packed format to RPF and from RPF to packed format
is included in the both the recursive packed factor and solve routines. The
subroutines DPPTRF, ZPPTRF, DPPTRS and ZPPTRS use packed
data format, and DPOTRF and DPOTRS use full data format. Figure 20
compares all three algorithms RPC, LAPACK full storage and LAPACK packed
storage.
Every gure has two subgures and one caption. The upper subgure shows
comparison curves for Cholesky factorization. The lower subgures show comparison
curves of forward and backward substitutions. The captions describe
details of the performance gures. The rst seven gures
describe the same comparison of performance on several dierent computers.
5.1 The IBM SMP PowerPC

Figure

11 shows the performance on the the IBM 4-way PowerPC 604e 332 MHz
computer.
The LAPACK routine DPPTRF (the upper subgure) performs at about
100 MFlops. Performance of the 'U' graph is a little better than the 'L' graph.
Performance remains constant as the order of the matrix increases.
The performance of the RPC factorization routine increases as n increases.
The 'U' graph increases from 50 MFlops to almost 600 MFlops and the 'L' graph
from 200 MFlops to 650 MFlops. The 'U' graph performance is better than the
'L' graph performance. The relative ('U', 'L') RPC algorithm performance is
(4.9, 7.2) times better than the DPPTRF algorithm for large matrix sizes.
The performance of the RPC solution routine (the lower subgure) for the
'L' and 'U' graphs are almost equal. The DPPTRS routine performs about
100 MFlops for all matrix sizes. The RPC algorithm curve increases from 250
MFlops to almost 800 MFlops. The relative ('U', 'L') performance of the RPC
algorithm is (5.7, 5.5) times faster than the DPPTRS algorithm for large matrix
sizes.
The matrix size varies from 300 to 3000 on these subgures.
5.2 The IBM Power2

Figure

12 shows the performance on the IBM Power2 160 MHz computer.
The LAPACK routine DPPTRF (the upper subgure) 'U' graph performs
at about 200 MFlops, the 'L' graph performs at about 150 MFlops. There is no
increase in both graphs as the size of the matrix grows.
The performance graphs of the RPC factorization routine both increase, the
'U' graph from 300 to a little more than 400 MFlops, and the 'L' graph from
200 MFlops to 450 MFlops. The 'L' graph is better than the 'U' graph when
the matrix sizes are between 750 and 3000. The 'U' graph is better than the
'L' graph when the matrix sizes are between 300 and 750. Both graphs grow
very rapidly for matrix sizes between 300 and 500. The relative ('U', 'L') RPC
algorithm performance is (1.9, 3.1) times faster than the DPPTRF algorithm
for large matrix sizes.
The performance of the RPC solution routine (the lower subgure) for the
'L' and 'U' graphs are almost equal. The performance of the DPPTRS algorithm
stays constant at about 250 MFlops decreasing slightly as n ranges from 300 to
3000. The performance of the RPC algorithm increases from 350 to more than
500 MFlops. The relative ('U', 'L') RPC algorithm performance is (2.3, 2.3)
times faster than the DPPTRS algorithm for large matrix sizes.
The matrix size varies from 300 to 3000 on these subgures.
5.3 The Compaq Alpha EV6

Figure

13 shows the performance on the the COMPAQ Alpha EV6 500 MHz
computer.
The LAPACK routine DPPTRF (the upper subgure) 'U' graph performs
better than the 'L' graph. The dierence is about 50 MFlops. The performance
starts at about 300 MFlops, increases to 400 MFlops and than drops down to
about 200 MFlops.
The performance of the RPC factorization routine increases as n increases.
Both graphs (the 'U' and 'L' graphs) are almost equal. The 'U' graph is a
little higher for matrix sizes between 300 and 450. The relative ('U', 'L') RPC
algorithm performance is (3.4, 5.0) times faster than the DPPTRF algorithm
for large matrix sizes.
For the routine DPPTRS the shape of the solution performance curves (the
lower subgure) for the 'L' and 'U' graphs are almost equal. The performance
of the DPPTRS routine decreases from 450 to 250 MFlops as n increases from
300 to 3000. The RPC performance curves increases from about 450 MFlops
to more than 750 MFlops. The performance ('U', 'L') of the RPC algorithm is
(3.3, 3.3) times faster than DPPTRS algorithm for large matrix sizes.
The matrix size varies from 300 to 3000 on these subgures.
5.4 The SGI R10000

Figure

14 shows the performance on the the SGI R10000 195 MHz computer,
on one processor only.
The LAPACK routine DPPTRF (the upper subgure) 'U' graph performs
better than the 'L' graph for matrix sizes from 300 to about 2000, after which
both the 'U' and the 'L' graphs are the same. The DPPTRF performance slowly
decreases.
The performance of the RPC factorization routine ('U' and 'L' graphs) increases
from about 220 to 300 MFlops as n increases from 300 to about 1000,
and stays constant as n increase to 3000. The relative ('U', 'L') RPC algorithm
performance is (4.9, 4.9) times faster than the DPPTRF algorithm for large
matrix sizes.
For the routine DPPTRS the shape of the solution performance curves (the
lower subgure) for the 'L' and 'U' graphs are almost equal. The performance
of the DPPTRS routine decreases from 130 MFlops to 60 MFlops as n increases
from 300 to 3000. The Performance of the RPC solution routine increases in the
beginning, and then runs constantly at about 300 MFlops. The performance
('U', 'L') of the RPC algorithm is (5.1, 5.2) times faster than the DPPTRS
algorithm for large matrix sizes.
The matrix size varies from 300 to 3000 on these subgures.
5.5 The SUN UltraSparc II

Figure

15 shows the performance on the the SUN UltraSparc II 400 MHz computer

The LAPACK routine DPPTRF (the upper subgure) 'U' and 'L' graphs
show almost equal performance when n > 1500. These functions start between
200 and 225 MFlops and then decrease down to about 50 MFlops.
For the RPC factorization routine, the performance of the 'U' and 'L' graphs,
are also almost equal over the whole interval. Their function values start from
quickly rise to 350 MFlops and then slowly increase to about 450
MFlops. The RPC factorization ('U', 'L') algorithm is (9.7, 10.2) times faster
than the DPPTRF algorithm for large matrix sizes.
The performance of the RPC solution routine (the lower subgure) for the 'L'
and 'U' graphs are almost equal. The DPPTRS performance graphs decreases
from 225 to 50 MFlops. The performance for the RPC solution graphs increases
from 330 to almost 450 MFlops. The RPC solution ('U', 'L') algorithm is (10.0,
times faster than the DPPTRS algorithm for large matrix sizes.
The matrix size varies from 300 to 3000 on these subgures.
5.6 The HP PA-8500

Figure

shows the performance on the the HP PA-8500 440 MHz computer.
The LAPACK routine DPPTRF (the upper subgure) 'U' and 'L' graphs
are decreasing functions. The 'U' graph function values go from about 370 to
100 MFlops. The 'L' graph function goes from 280 to about 180 MFlops.
The performance of the RPC factorization graphs are increasing functions as
the matrix size increases from 1000 to 3000. The performance varies for matrix
sizes between 500 and 1500. The 'U' graph function values range from about
700 MFlops to almost 800 MFlops, the 'L' graph function values range from 600
MFlops to a little more than 700 MFlops. The RPC algorithm ('U', 'L') is (4.7,
times faster than the DPPTRF algorithm for large matrix sizes.
The performance of the RPC solution routine (the lower subgure) for the 'L'
and 'U' graphs are almost equal. The DPPTRS routine performance decreases
from 300 MFlops to 200 MFlops. The RPC algorithm curve increases from 550
MFlops to almost 810 MFlops. The RPC algorithm ('U', 'L') is (5.2, 5.0) times
faster than the DPPTRS algorithm for large matrices in the solution case.
The matrix size varies from 500 to 2500 on these subgures.
5.7 The INTEL Pentium III

Figure

17 shows the performance on the INTEL Pentium III 500 MHz computer.
The LAPACK routine DPPTRF (the upper subgure) 'U' and 'L' graphs
are decreasing functions. The 'U' graph function ranges from about 100 to 80
MFlops. The 'L' graph function ranges from less than 50 to about 25 MFlops.
For the RPC factorization routine the 'U' and the 'L' graphs are almost
equal. The graphs are increasing functions from about 200 to 310 MFlops.
The RPC factorization algorithm ('U', 'L') is (4.2, 9.2) times faster than the
DPPTRF algorithm for large matrices.
The performance of the RPC solution routine (the lower subgure) for the 'L'
and 'U' graphs are almost equal. The DPPTRS performance graphs decreases
from about 80 to about 50 MFlops. The RPC algorithm curves increases from
240 to about 330 MFlops. The RPC algorithm ('U', 'L') is (5.9, 6.0) times faster
than the DPPTRS algorithm for large matrices.
The matrix size varies from 500 to 3000 on these subgures.
5.8 The IBM SMP PowerPC with OpenMP directives

Figure

shows the performance on the the IBM 4-way PowerPC 604e 332 MHz
computer.
These graphs demonstrate successful use of OpenMP[17, 18] parallelizing
directives. The curves LAPACK(L), LAPACK(U), Recursive(L) and Recur-
are identical to the corresponding curves of gure 11. We compare curves
Recursive(L), Recursive(U), Rec.Par(L) and Rec.Par(U). The Rec.Par(L) and
Rec.Par(U) curves result from double parallelization. The RPC algorithms call
a parallelized DGEMM and they are parallelized themselves by the OpenMP
directives.
The Rec.Par(L) curve is not much faster than Recursive(L), sometimes it
is slower. The Rec.Par(U) is the fastest, specially for large size matrices. The
doubly parallelized RPC algorithm (Rec.Par(U)) is about 100 MFlops faster
than the ordinary RPC algorithm (Recursive(U)). The relative ('U', 'L') RPC
factorization algorithm performance is (5.6, 7.6) times faster than the DPPTRF
algorithm for large matrices.
The RPC double parallelization algorithm for the solution (lower subgure)
exceeds 800 MFlops. The relative ('U', 'L') RPC solution algorithm performance
is (6.5, 6.6) times faster than the DPPTRF algorithm for large matrices.
The matrix size varies from 300 to 3000 on these subgures.
5.9 The INTEL Pentium III running Complex Arithmetic

Figure

19 shows the performance on the INTEL Pentium III 500 MHz computer.
This gure demonstrate the successful use of RPC algorithm for Hermitian
positive denite matrices. The performance is measured in Complex MFlops.
To compare with the usual real arithmetic MFlops the Complex MFlops should
be multiplied by 4.
The LAPACK routine ZPPTRF (the upper subgure) 'U' graph performs a
little better than the 'L' graph. These routine performs at about 80 MFlops.
The RPC Hermitian factorization routine 'U' graph performs better than
the 'L' graph. The RPC performance graphs are increasing functions. They
go from 240 up to 320 MFlops. The RPC Hermitian factorization algorithm
('U', 'L') is (3.8, 4.3) times faster than the ZPPTRF algorithm for the large size
matrices.
The performance of the RPC solution routine (the lower subgure) for the
'L' and 'U' graphs are almost equal. The ZPPTRS performance decreases from
about 108 to 80 MFlops. The RPC solution algorithm increases from about 240
up to more than 320 MFlops. The RPC algorithm ('U', 'L') is (3.9, 3.7) times
faster than the ZPPTRS algorithm for large Hermitian matrices.
The matrix size varies from 500 to 3000 on these subgures.
5.10 The INTEL Pentium III with all three Cholesky Algorithm

Figure

20 shows the performance on the INTEL Pentium III 500 MHz computer.
The graphs on this gure depict all three Cholesky algorithms, the LAPACK
full storage (DPOTRF and DPOTRS) algorithms, the LAPACK packed storage
(DPPTRF and DPPTRS) algorithms and the RPC (factorization and solution)
algorithms.
The LAPACK packed storage algorithms (DPPTRF and DPPTRS) are previously
explained on gure 17.
The DPOTRF routine (the upper subgure), for both the 'U' and 'L' cases,
performs better than the RPC factorization routine for smaller matrices. For
larger matrices the RPC factorization algorithm performs equally well or slightly
better than the DPOTRF algorithm.
The performance of the DPOTRS algorithms ('U' and 'L' graphs) are better
than the RPC performance for this computer.
However, the POTRF and POTRS storage requirement is almost twice the
storage requirement of the RPC algorithms.
The matrix size varies from 500 to 3000 on these subgures.
6 Conclusion
We summarize and emphasize the most important developments described in
our paper.
A recursive packed Cholesky factorization algorithm based on BLAS Level 3
operations has been developed.
The RPC factorization algorithm works with almost the same speed as
the traditional full storage algorithm but occupies the same data storage
as the traditional packed storage algorithm. Also see bullet 4.
The user interface of the new packed recursive subroutines (RP PPTRF
and RP PPTRS) is exactly the same as the traditional LAPACK sub-routines
(PPTRF and PPTRS). The user will see identical data formats.
However, the new routines run much faster.
Two separate routines are described here: RP PPTRF and RP PPTRS.
The data format is always converted from LAPACK packed data format
to the recursive packed data format before the routine starts its operation
and converted back to LAPACK data format afterwards. The RP PPSV
subroutine exists in our package which is equivalent to the LAPACK PPSV
routine. In the RP PPSV subroutine the data is not converted between
the factorization and the solution.
New recursive packed Level 3 BLAS, RP TRSM and RP SYRK, written
in Fortran90[21] were developed. They only call the GEMM routine.
This GEMM subroutine can be developed either by the computer manufacturer
or generated by ATLAS system[25]. The ATLAS generated GEMM
subroutine is usually compatible with the manufacturer developed routine.

Acknowledgements

This research was partially supported by the LAWRA project, the UNIC collaboration
with the IBM T.J. Watson Research Center at Yorktown Heights.
The last two authors were also supported by the Danish Natural Science Re-search
Council through a grant for the EPOS project (E-cient Parallel Algorithms
for Optimization and Simulation).



--R

Exploiting functional parallelism on power2 to design high-performance numerical algorithms



Applied Numerical Linear Algebra.
BLAS (Basic Linear Algebra Subprograms).


An extended set of FORTRAN basic linear algebra subroutines.

Matrix Computations.


Recursion leads to automatic variable blocking for dense linear-algebra algorithms
Accuracy and Stability of Numerical Algorithms.




Basic linear algebra subprograms for Fortran usage.
FORTRAN 90/95 Explained.
Locality of Reference in LU Decomposition with Partial Pivoting.
Numerical Linear Algebra.

Automatically Tuned Linear Algebra Software (ATLAS).
--TR
An extended set of FORTRAN basic linear algebra subprograms
A set of level 3 basic linear algebra subprograms
Exploiting functional parallelism of POWER2 to design high-performance numerical algorithms
Matrix computations (3rd ed.)
Optimizing matrix multiply using PHiPAC
Applied numerical linear algebra
Locality of Reference in LU Decomposition with Partial Pivoting
Recursion leads to automatic variable blocking for dense linear-algebra algorithms
GEMM-based level 3 BLAS
Fortran 90/95 explained (2nd ed.)
Basic Linear Algebra Subprograms for Fortran Usage
Accuracy and Stability of Numerical Algorithms
Numerical Linear Algebra for High Performance Computers
Recursive Formulation of Cholesky Algorithm in Fortran 90
Superscalar GEMM-based Level 3 BLAS - The On-going Evolution of a Portable and High-Performance Library
Recursive Blocked Data Formats and BLAS''s for Dense Linear Algebra Algorithms

--CTR
Chan , Enrique S. Quintana-Orti , Gregorio Quintana-Orti , Robert van de Geijn, Supermatrix out-of-order scheduling of matrix operations for SMP and multi-core architectures, Proceedings of the nineteenth annual ACM symposium on Parallel algorithms and architectures, June 09-11, 2007, San Diego, California, USA
Steven T. Gabriel , David S. Wise, The Opie compiler from row-major source to Morton-ordered matrices, Proceedings of the 3rd workshop on Memory performance issues: in conjunction with the 31st international symposium on computer architecture, p.136-144, June 20-20, 2004, Munich, Germany
Dror Irony , Gil Shklarski , Sivan Toledo, Parallel and fully recursive multifrontal sparse Cholesky, Future Generation Computer Systems, v.20 n.3, p.425-440, April 2004
Bjarne S. Andersen , John A. Gunnels , Fred G. Gustavson , John K. Reid , Jerzy Waniewski, A fully portable high performance minimal storage hybrid format Cholesky algorithm, ACM Transactions on Mathematical Software (TOMS), v.31 n.2, p.201-227, June 2005
F. G. Gustavson, High-performance linear algebra algorithms using new generalized data structures for matrices, IBM Journal of Research and Development, v.47 n.1, p.31-55, January
complex Hermitian matrices;cholesky factorization and solution;real symmetric matrices;recursive algorithms;novel packed matrix data structures;BLAS;positive definite matrices