A C program for solving banded, symetric, positive definite systems

Abstract

Banded, symmetric, positive definite systems occur frequently. Although they can be solved by any program that incorporates Cholesky decomposition, or even LU decomposition, for greatest efficiency account should taken of the specialized nature of banded, symmetric, positive definite systems: elements outside the band need not be stored; because of symmetry nearly half of the non-null elements are superfluous and likewise need not be stored; positive definiteness obviates row and column interchanging. A C program, employing square root-free Cholesky decomposition, is given for solving such systems. For sparse matrix processing C offers advantages over FORTRAN: when the two-dimensional array is represented by a vector, C's conditional operator permits a compact way to index elements in the vector; the vector itself may be allocated dynamically and stores no null elements; divisions by powers of two may be replaced by bit shifts to the right; and preprocessor macros, substituting in-line code for subroutine calls, implement functions. An example, a 40 x 40 pentadiagonal system, illustrates use of the program. In this example the C code for banded, symmetric, positive definite systems is four hundred times faster than code for general Cholesky decomposition. Although accumulation of matrix inner products to high precision is usually not required for banded systems, a function subroutine is nevertheless presented to sum floating-point numbers with cascaded accumulators.