DGE - Double precision Gaussian Elimination matrix subroutines adapted
from LINPACK FORTRAN:
dgefa Gaussian elimination to obtain the LU factorization of a matrix.
dgeco Gaussian elimination to obtain the LU factorization and condition
number of a matrix.
dgesl Solve linear system Ax = b or A'x = b after LU factorization.
Function Prototypes:
void dgefa (double **a, int n, int *ipvt, int *info);
void dgeco (double **a, int n, int *ipvt, double *rcond, double *z);
void dgesl (double **a, int n, int *ipvt, double *b, int job);
dgfa:
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)
dgeco:
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 of condition number (see notes below)
Workspace:
z array[n]
dgesl:
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 dgesl.
dgeco:
Given the reciprocal of the condition number, rcond, and the double
epsilon, DBL_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/DBL_EPSILON)
Author: Dave Hale, Colorado School of Mines, 1989