http://www.ccl.net/cca/software/SOURCES/C/kinetics1/dgefa.shtml
CCL dgefa
```extern void daxpy(), dscal();
extern int idamax();

void dgefa( a, n, ipvt, info )

double **a;
int n, *ipvt, *info;

/*
Purpose : dgefa factors a double matrix by Gaussian elimination.

dgefa is usually called by dgeco, but it can be called directly
with a saving in time if rcond is not needed.
(Time for dgeco) = (1+9/n)*(time for dgefa).

This c version uses algorithm kji rather than the kij in dgefa.f.
Note that the fortran version input variable lda is not needed.

On Entry :

a   : double matrix of dimension ( n+1, n+1 ),
the 0-th row and column are not used.
a is created using NewDoubleMatrix, hence
lda is unnecessary.
n   : the row dimension of a.

On Return :

a     : a lower triangular matrix and the multipliers
which were used to obtain it.  The factorization
can be written a = L * U where U is a product of
permutation and unit upper triangular matrices
and L is lower triangular.
ipvt  : an n+1 integer vector of pivot indices.
*info : = 0 normal value,
= k if U[k][k] == 0.  This is not an error
condition for this subroutine, but it does
indicate that dgesl or dgedi will divide by
zero if called.  Use rcond in dgeco for
a reliable indication of singularity.

Notice that the calling program must use &info.

BLAS : daxpy, dscal, idamax
*/

{
int j, k, i;
double t;

/* Gaussian elimination with partial pivoting.   */

*info = 0;
for ( k = 1 ; k <= n - 1 ; k++ ) {
/*
Find j = pivot index.  Note that a[k]+k-1 is the address of
the 0-th element of the row vector whose 1st element is a[k][k].
*/
j = idamax( n-k+1, a[k]+k-1, 1 ) + k - 1;
ipvt[k] = j;
/*
Zero pivot implies this row already triangularized.
*/
if ( a[k][j] == 0. ) {
*info = k;
continue;
}
/*
Interchange if necessary.
*/
if ( j != k ) {
t = a[k][j];
a[k][j] = a[k][k];
a[k][k] = t;
}
/*
Compute multipliers.
*/
t = -1. / a[k][k];
dscal( n-k, t, a[k]+k, 1 );
/*
Column elimination with row indexing.
*/
for ( i = k + 1 ; i <= n ; i++ ) {
t = a[i][j];
if ( j != k ) {
a[i][j] = a[i][k];
a[i][k] = t;
}
daxpy( n-k, t, a[k]+k, 1, a[i]+k, 1 );
}
}                     /*  end k-loop  */

ipvt[n] = n;
if ( a[n][n] == 0. )
*info = n;

}

```
 Modified: Mon May 13 16:00:00 1991 GMT Page accessed 5905 times since Sat Apr 17 21:32:45 1999 GMT