New generalized data structures for matrices lead to a variety of high-performance algorithms

We describe new data structures for full and packed storage of dense symmetric/triangular arrays that generalize both. Using the new data structures, one is led to several new algorithms that save "half" the storage and outperform the current blocked-based level-3 algorithms in LAPACK. We concentrate on the simplest forms of the new algorithms and show for Cholesky factorization they are a direct generalization of LINPACK. This means that level-3 BLAS's are not required to obtain level-3 performance. The replacement for level-3 BLAS are so-called kernel routines, and on IBM platforms they are producible from simple textbook type codes, by the XLF Fortran compiler. In the sequel I will label these "vanilla" codes. The results for Cholesky, on Power3 with a peak performance of 800 Mflop/s at n greater than or equal to 200 is over 720 MFlop/s and reaches 735 MFlop/s. Using conventional full-format LAPACK DPOTRF with ESSL BLAS's, one first gets 600 MFlop/s at n 600 and only reaches a peak of 620 MFlop/s. We have also produced simple square blocked full-matrix data formats where the blocks themselves are stored in column-major (Fortran) order or row-major (C) format. The simple algorithms of LU factorization with partial pivoting for this new data format is a direct generalization of LINPACK algorithm DGEFA. Again, no conventional level-3 BLAS's are required; the replacements are again so-called kernel routines, Programming far squared blocked full-matrix format can be accomplished in standard Fortran through the use of three- and four-dimensional arrays. Thus, no new compiler support is necessary. Finally we mention that other more complicated algorithms are possible, for example, recursive ones. The recursive algorithms are also easily programmed via the use of tables that address where the blocks are stored in the two-dimensional recursive block array.