Basic Linear Algebra Subroutines Research Articles

Restructuring Software: A Case Study

We use knot count and path count metrics to identify which routines in the Level 1 basic linear algebra subroutines (BLAS) might benefit from code restructuring. We then consider how logical restructuring and the improvements in the facilities available from successive versions of Fortran have allowed us to improve the complexity of the code as measured by knot count, path count and cyclomatic complexity, and the user interface of one of the identified routines which computes the Euclidean norm of a vector. With these reductions in complexity we hope that we have contributed to improvements in the maintainability and clarity of the code. Software complexity metrics and the control graph are used to quantify and provide a visual guide to the quality of the software, and the performance of two Fortran code restructuring tools is reported. Finally, we give some indication of the cost of the extra numerical robustness offered by the BLAS routine over the use of new Fortran 90 intrinsic functions.

Computing the MDM T decomposition

The MDM T factorization of an n×n symmetric indefinite matrix A can be used to solve a linear system with A as the coefficient matrix. This factorization can be computed efficiently using an algorithm given in 1977 by Bunch and Kaufman. The LAPACK project has been implementing block versions of well-known algorithms for solving dense linear systems and eigenvalue problems. The block version of the MDM T decomposition algorithm in LAPACK requires the user to specify a block size b by supplying an n×b scratch array. It then makes ( n/b −2) 2 /2 invocations of a matrix-matrix product subroutine with one matrix no larger than b×b, n 2 /(2 b )− b invocations of a matrix-vector product routine with a matrix no larger than b×b , and between n − n/b and 2( n − n/b ) invocations of a matrix-vector product routine with matrices with less than b columns. Because the user can query LAPACK about an optimal block size, our concern is focused on users who cannot change the amount of available scratch space or who neglect to use this facility and are unaware of a performance degradation with small block sizes. This article suggests two alternative algorithms. The first is a block algorithm requiring b×b scratch space and is about 5% slower than LAPACK's current block algorithm with large b . The user does not have to specify the block size. The second algorithm is a rejuvenation of an old implementation of the MDM T decomposition algorithm that requires n matrix-vector products. The performance of the various algorithms on a specific machine is dependent on the manufacturer's implementation of the different basic linear algebra subroutines that they each invoke. Our data indicate that on the Cray Y-MP, Alliant, and Convex the time for the rejuvenated algorithm is either less than or within 10% of that of LAPACK's block algorithm with large b .

Multiplication of matrices of arbitrary shape on a data parallel computer

Some level-2 and level-3 Distributed Basic Linear Algebra Subroutines (DBLAS) that have been implemented on the Connection Machine system CM-200 are described. No assumption is made on the shape or size of the operands. For matrix-matrix multiplication, both the nonsystolic and the systolic algorithms are outlined. A systolic algorithm that computes the product matrix in-place is described in detail. We show that a level-3 DBLAS yields better performance than a level-2 DBLAS. On the Connection Machine system CM-200, blocking yields a performance improvement by a factor of up to three over level-2 DBLAS. For certain matrix shapes the systolic algorithms offer both improved performance and significantly reduced temporary storage requirements compared to the nonsystolic block algorithms. We show that, in order to minimize the communication time, an algorithm that leaves the largest operand matrix stationary should be chosen for matrix-matrix multiplication. Furthermore, it is shown both analytically and experimentally that the optimum shape of the processor array yields square stationary submatrices in each processor, i.e. the ratio between the length of the axes of the processing array must be the same as the ratio between the corresponding axes of the stationary matrix. The optimum processor array shape may yield a factor of five performance enhancement for the multiplication of square matrices. For rectangular matrices a factor of 30 improvement was observed for an optimum processor array shape compared to a poorly chosen processor array shape.

Block-Cyclic Dense Linear Algebra

Block-cyclic order elimination algorithms for LU and OR factorization and solve routines are described for distributed memory architectures with processing nodes configured as two-dimensional arrays of arbitrary shape. The cyclic-order elimination, together with a consecutive data allocation, yields good load balance for both the factorization and solution phases for the solution of dense systems of equations by LU and OR decomposition. Blocking may offer a substantial performance enhancement on architectures for which the level-2 or level-3 BLAS (basic linear algebra subroutines) are ideal for operations local to a node. High-rank updates local to a node may have a performance that is a factor of four or more higher than a rank-1 update.This paper shows that in many parallel implementations, the $O(N^2 )$ work in the factorization may be of the same significance as the $O(N^3 )$ work, even for large matrices. The $O(N^2 )$ work is poorly load balanced in two-dimensional nodal arrays, which are shown to b...

Local Basic Linear Algebra Subroutines (Lblas) for Distributed Memory Architectures and Languages With Array Syntax

We describe a subset of the level-1, level-2, and level-3 BLAS implemented for each node of the Connection Machine system CM-200. The routines, collectively called LBLAS, have interfaces consistent with lan guages with an array syntax such as Fortran 90. One novel feature, important for distributed memory archi tectures, is the capability of performing computations on multiple instances of objects in a single call. The number of instances and their allocation across mem ory units, and the strides for the different axes within the local memories, are derived from an array descrip tor that contains type, shape, and data distribution in formation. Another novel feature of the LBLAS is a se lection of loop order for rank-1 updates and matrix- matrix multiplication based on array shapes, strides, and DRAM page faults. The peak efficiencies for the routines are in excess of 75%. Matrix-vector multiplica tion achieves a peak efficiency of 92%. The optimiza tion of loop ordering has a success rate exceeding 99.8% for matrices for which the sum of the lengths of the axes is at most 60. The success rate is even higher for all possible matrix shapes. The performance loss when a nonoptimal choice is made is less than ∼15% of peak and typically less than 1% of peak. We also show that the performance gain for high-rank updates may be as much as a factor of 6 over rank-1 updates.

Parallel Decomposition of Multicommodity Network Flows Using a Linear-Quadratic Penalty Algorithm

The central theme of this paper is the design and evaluation of a parallel decomposition algorithm for multicommodity network flow problems based on the notion of piecewise linear-quadratic penalty (LQP) functions. The algorithm induces separability of both the constraint set and the objective function by commodity during the subproblem phase. Hence it is suitable for implementations on coarse grain parallel architectures. A master problem, of significantly smaller dimension than the original problem, uses dense linear algebra computations. Hence it exploits a vector architecture, but is also suited for parallel implementation and can be executed using Basic Linear Algebra (BLAS) subroutines. Computational results on a CRAY Y-MPE264 supercomputer with a set of large multicommodity network flow problems drawn from a military application are presented and analyzed. INFORMS Journal on Computing, ISSN 1091-9856, was published as ORSA Journal on Computing from 1989 to 1995 under ISSN 0899-1499.

Parallel matrix multiplication on networked microcomputers

We present methods for utilizing parallel processing capability of the idle microcomputers on a local area network to perform computationally intensive operations frequently encountered in linear algebra. By making use of the data communication properties of the network, parallel algorithms for multiplication of large matrices have been designed and implemented. The techniques presented in this paper can be used to obtain efficient implementations of basic linear algebra subroutines.

A parallel and vectorial implementation of basic linear algebra subroutines in iterative solving of large sparse linear systems of equations

Electromagnetic field analysis by finite element methods, which involve the solution of large sparse systems of linear equations, is discussed. Though no discernible structure for the distribution of nonzero elements can be found (e.g. multidiagonal structures), subsets of independent equations can be determined. Equations that are in the same subset are then solved in parallel. A good choice for the storage scheme of sparse matrices is also very important to speed up the resolution by vectorization. The modifications made to data structures are presented, and the possibility of using other schemes is discussed.< <ETX xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">&gt;</ETX>

LU decomposition of matrices with augmented dense constraints

AbstractSparse matrices composed of a central band and augmented dense rows and columns are becoming prevalent in the numerical solution of a large class of boundary and initial‐value problems. A Fortran Subroutine ARROW is presented for the LU decomposition and solution of linear equation systems with such a structure. The computational speed of the program is compared in MFLOPS (millions of floating point operations per second) to the LINPACK benchmark for the solution of a dense linear system and is found to be of comparable speed on both supercomputers and minicomputers. Use of the Basic Linear Algebra Subroutines (BLAS) available on most machines significantly enhances the speed of ARROW.