A library of linear algebra utilities.

Support for linear algebra comes in the following components:

provides a suite of wrappers for the basic arithmetic operations on members of a field capable of subsuming the reals. This is provided for the purpose of enabling an application to be compiled with a local replacement of scalar.h which must provide the same functionality but may provide it for a different scalar field, such as complex, and still be able to compile up the other linear algebra code.
provides low-level array/matrix operations. By convention, most routines take a first argument (usually called `spare') which may either be null or point to a spare object of the type (and size) of the function return: this is to allow re-use of dead variables so as to reduce the number of creation and destruction actions required.
wraps up array in a way designed to make clear the distinction between (`column') vectors and (`row') co-vectors, especially at the `matrix' level where it is common practice to use one representation (the matrix) for four singularly distinct kinds of things, only two of which it resembles.
provides I/O routines for the array and vector modules.

My intent in designing the algorithms used here has been to make them reasonably straightforward and robust. These algorithms have no pretensions to being especially fast, except the Slatec pieces listed below, though some contain (probably naive) attempts at optimisation (to suppress these, -DDefinitive). Eventually, I would like to implement a compile-time switch to enable the code to use calls to the Basic Linear Algebra Subprograms (BLAS) or some other library of highly-optimised code with similar functionality, when available: the current implementation should be regarded as a fall-back option in the absence of such a library.

As an exception to the above, I have translated some of the public-domain code provided as the Slatec library. In time, I may realise my goal of doing things with BLAS by the simple expedient of translating all the Slatec implementations. My translation of the ForTran (itself, in some cases, obtained by translating Algol) may have added a few bugs and is likely to have reduced the efficiency. As derivative works of the Slatec library, these are public domain.

Singular Value Decomposition is used to obtain inverses to matrices. This means that matrix division and inversion will succeed even when the matrices are singular: SVD delivers the means to produce a `best attempt' at inverting matrices. (DSVDC.f)
Eigen-Basis Decomposition is a special case of SVD for use with symmetric matrices: special functions are provided for division by and inversion of symmetric matrices. (TRED2.f, IMTQL2.f)

Spare pointers

Most functions in this module take a `spare' first argument. The purpose of this is to save on calls to the allocation (*_new) and release (*_kill) functions: however, your code should always assume that any value other than null passed as a `spare' pointer will be released and the return from the function will be some totally different, freshly allocated pointer. At the same time, any value other than null which is passed as a `spare' pointer must actually be valid for the function return: this just means that it must be a heap pointer with the correct type and dimension(s).

It should be noted, in conjunction with this, that the other (pointer) arguments to such functions are always read-only. This does not prevent you from passing a pointer both as a data argument and as the `spare': however, I have not yet done the necessary tampering with my algorithms to ensure that this will work. In general, the purely first-rank operations and the simplest (add, subtract, scale) second rank operations are safe. Transposition definitely isn't safe, nor is any contraction involving second rank entities.

Maintained by Eddy.
$Id: linear.html,v 1.2 2001-11-04 18:28:34 eddy Exp $