We present a new set of codes for solving almost block diagonal systems of linear equations and for performing multiplicative operations with matrices represented using the same data structures. These data structures arise when solving ordinary differential equation boundary value problems with non-separated boundary conditions by finite differences, and when using spline collocation methods. Our codes are written in a modular form using the BLAS and are intended to take advantage of vector architecture and, to a limited extent, parallelism.
1. Introduction
Recently, there have appeared implementations of a number of algorithms for solving almost block diagonal (ABD) systems of linear equations Ax=b, where A eR nΓn and x, beRn.
(1)
The ABD form is illustrated in Fig. 1. The precise definition is given in Ref.
[1] and in the prologue to the software described in this paper. A is said to be stored in ABD format if it is stored in an array of one dimension block by block, each block being stored column by column. Systems of this type arise in a wide variety of applications. In particular, they arise directly when solving ordinary differential equation boundary value problems with separated boundary conditions, using any "global method" based on local discretization. For example, they arise straightforwardly from finite-difference discretizations and also in spline collocation techniques [2,3]. An analysis of the cost of using ABD software in this context is given in Ref. [4]. ABDs also arise when discretizing time-dependent partial differential equations in one space dimension by the method of lines. Indeed, Keast and Muir [5] have modified the well-known PDECOL code [6] to use ABD software for this data structure. The software is entirely appropriate to this application. The original PDECOL used banded matrix structures hence automatically carrying significant fill-in, implying an inefficiency in both memory and arithmetic operations. Algorithms for solving equation (1) when A is ABD have been available for some time (see, for example Ref. [7]) but portable FORTRAN 77 implementations in the public domain are of more recent origin. The code SOLVEBLOK [8] used row eliminations, with row interchanges to ensure numerical stability. This algorithm clearly leads to some fill-in. In contrast, the ARCECO family of codes [1] uses row and column eliminations and row and column interchanges. This family of codes is designed to produce no fill-in. Similarly designed is the unpublished code LAMPAK, discussed in Ref. [1], for the same problem, which employs row eliminations and row and column interchanges. Indeed, for both the latter sets of codes, whether to use a row or column interchange at the next stage of elimination is determined as the unique choice which avoids fill-in tThe work of the second author was partly carried out at the Numerical Algorithms Group Ltd, and at the