high-performance computers. The BLAS enable LAPACK routines to achieve high performance with portable code.

The BLAS are not strictly speaking part of LAPACK, but Fortran 77 code for the BLAS is distributed with LAPACK, or can be obtained separately from netlib where “model implementations” are found.

The model implementation is not expected to perform as well as a specially tuned implementation on most high-performance computers – on some machines it may give much worse performance – but it allows users to run LAPACK codes on machines that do not o er any other implementation of the BLAS.

The complete LAPACK package or individual routines from LAPACK are freely available from the World Wide Web or by anonymous ftp. The LAPACK homepage can be accessed via the URL http://www.netlib.org/lapack/.

5.2BLAS

By 1976 it was clear that some standardization of basic computer operations on vectors was needed [9]. By then it was already known that coding procedures that worked well on one machine might work very poorly on others. In consequence of these observations, Lawson, Hanson, Kincaid and Krogh proposed a limited set of Basic Linear Algebra Subprograms (BLAS) to be (hopefully) optimized by hardware vendors, implemented in assembly language if necessary, that would form the basis of comprehensive linear algebra packages [9]. These so-called Level 1 BLAS consisted of vector operations and some attendant co-routines. The first major package which used these BLAS kernels was LINPACK [3]. Soon afterward, other major software libraries such as the IMSL library and NAG rewrote portions of their existing codes and structured new routines to use these BLAS. Early in their development, vector computers saw significant optimizations using the BLAS. Soon, however, such machines were clustered together in tight networks and somewhat larger kernels for numerical linear algebra were developed [6, 7] to include matrix-vector operations (Level 2 BLAS). Additionally, FORTRAN compilers were by then optimizing vector operations as e ciently as hand coded Level 1 BLAS. Subsequently, in the late 1980s, distributed memory machines were in production and shared memory machines began to have significant numbers of processors. A further set of matrix-matrix operations was proposed [4] and soon standardized [5] to form a Level 3. The first major package for linear algebra which used the Level 3 BLAS was LAPACK [1] and subsequently a scalable (to large numbers of processors) version was released as ScaLAPACK [2]. Vendors focused on Level 1, Level 2, and Level 3 BLAS which provided an easy route to optimizing LINPACK, then LAPACK. LAPACK not only integrated pre-existing solvers and eigenvalue routines found in EISPACK [10] (which did not use the BLAS) and LINPACK (which used Level 1 BLAS), but incorporated the latest dense and banded linear algebra algorithms available. It also used the Level 3 BLAS which were optimized by much vendor e ort. Later, we will illustrate several BLAS routines. Conventions for di erent BLAS are indicated by

• A root operation. For example, axpy for the operation

•A prefix (or combination prefix) to indicate the datatype of the operands, for example saxpy for single precision axpy operation, or isamax for the index of the maximum absolute element in an array of type single.

•a su x if there is some qualifier, for example cdotc or cdotu for conjugated or unconjugated complex dot product, respectively:

where both x, y are vectors of complex elements.

Tables 5.1 and 5.2 give the prefix/su x and root combinations for the BLAS, respectively.

Table 5.1: Basic Linear Algebra Subprogram prefix/su x conventions.

5.2.1Typical performance numbers for the BLAS

Let us look at typical representations of all three levels of the BLAS, daxpy, ddot, dgemv, and dgemm, that perform some basic operations. Additionally, we look at the rank-1 update routine dger. An overview on the number of memory accesses and floating point operations is given in Table 5.3. The Level 1 BLAS comprise basic vector operations. A call of one of the Level 1 BLAS thus gives rise to O(n) floating point operations and O(n) memory accesses. Here, n is the vector length. The Level 2 BLAS comprise operations that involve matrices and vectors. If the involved matrix is n-by-n then both the memory accesses and the floating point operations are of O(n2). In contrast, the Level 3 BLAS have a higher order of floating point operations than memory accesses. The most prominent operation of the Level 3 BLAS, matrix-matrix multiplication costs O(n3) floating point operations while there are only O(n2) reads and writes. The last column in Table 5.3 shows the crucial di erence between the Level 3 BLAS and the rest.

Table 5.4 gives some performance numbers for the five BLAS of Table 5.3. Notice that the timer has a resolution of only 1 µsec! Therefore, the numbers in Table 5.4 have been obtained by timing a loop inside of which the respective function is called many times. The Mflop/s rates of the Level 1 BLAS ddot and daxpy quite precisely reflect the ratios of the memory accesses of the two routines, 2n vs. 3n. The high rates are for vectors that can be held in the on-chip cache of 512 MB. The low 240 and 440 Mflop/s with the very long vectors are related to the memory bandwidth of about 1900 MB/s.

The Level 2 BLAS dgemv has about the same performance as daxpy if the matrix can be held in cache (n = 100). Otherwise it is considerably reduced. dger has a high volume of read and write operations, while the number of floating point operations is limited.