SGA - Single precision general matrix functions adapted from LINPACK FORTRAN:

sgefa	Gaussian elimination to obtain the LU factorization of a matrix.
sgeco	Gaussian elimination to obtain the LU factorization and 
	condition number of a matrix.
sgesl	Solve linear system Ax = b or A'x = b after LU factorization.

Function Prototypes:
void sgefa (float **a, int n, int *ipvt, int *info);
void sgeco (float **a, int n, int *ipvt, float *rcond, float *z);
void sgesl (float **a, int n, int *ipvt, float *b, int job);

sgefa:
Input:
a		matrix[n][n] to be factored (see notes below)
n		dimension of a

Output:
a		matrix[n][n] factored (see notes below)
ipvt		indices of pivot permutations (see notes below)
info		index of last zero pivot (or -1 if no zero pivots)

sgeco:
Input:
a		matrix[n][n] to be factored (see notes below)
n		dimension of a

Output:
a		matrix[n][n] factored (see notes below)
ipvt		indices of pivot permutations (see notes below)
rcond		reciprocal condition number (see notes below)

Workspace:
z		array[n]

sgesl
Input:
a		matrix[n][n] that has been LU factored (see notes below)
n		dimension of a
ipvt		indices of pivot permutations (see notes below)
b		right-hand-side vector[n]
job		=0 to solve Ax = b
		=1 to solve A'x = b

Output:
b		solution vector[n]

Notes:
These functions were adapted from LINPACK FORTRAN.  Because two-dimensional 
arrays cannot be declared with variable dimensions in C, the matrix a
is actually a pointer to an array of pointers to floats, as declared
above and used below.

Elements of a are stored as follows:
a[0][0]    a[1][0]    a[2][0]   ... a[n-1][0]
a[0][1]    a[1][1]    a[2][1]   ... a[n-1][1]
a[0][2]    a[1][2]    a[2][2]   ... a[n-1][2]
.                                       .
.             .                         .
.                        .              .
.                                       .
a[0][n-1]  a[1][n-1]  a[2][n-1] ... a[n-1][n-1]

Both the factored matrix a and the pivot indices ipvt are required
to solve linear systems of equations via sgesl.

sgeco:
Given the reciprocal of the condition number, rcond, and the float
epsilon, FLT_EPSILON, the number of significant decimal digits, nsdd,
in the solution of a linear system of equations may be estimated by:
	nsdd = (int)log10(rcond/FLT_EPSILON)

This function was adapted from LINPACK FORTRAN.  Because two-dimensional 
arrays cannot be declared with variable dimensions in C, the matrix a
is actually a pointer to an array of pointers to floats, as declared
above and used below.

Author:  Dave Hale, Colorado School of Mines, 10/01/89