http://www.ccl.net/cca/software/SOURCES/C/kinetics1/dgesl.shtml
CCL dgesl
```extern double ddot();
extern void daxpy();

void dgesl( a, n, ipvt, b, job )

double **a, *b;
int n, *ipvt, job;

/*
Purpose : dgesl solves the linear system
a * x = b or Transpose(a) * x = b
using the factors computed by dgeco or degfa.

On Entry :

a    : double matrix of dimension ( n+1, n+1 ),
the output from dgeco or dgefa.
The 0-th row and column are not used.
n    : the row dimension of a.
ipvt : the pivot vector from degco or dgefa.
b    : the right hand side vector.
job  : = 0       to solve a * x = b,
= nonzero to solve Transpose(a) * x = b.

On Return :

b : the solution vector x.

Error Condition :

A division by zero will occur if the input factor contains
a zero on the diagonal.  Technically this indicates
singularity but it is often caused by improper argments or
improper setting of the pointers of a.  It will not occur
if the subroutines are called correctly and if dgeco has
set rcond > 0 or dgefa has set info = 0.

BLAS : daxpy, ddot
*/

{
int nm1, k, j;
double t;

nm1 = n - 1;

/*
Job = 0, solve a * x = b.
*/
if ( job == 0 ) {
/*
First solve L * y = b.
*/
for ( k = 1 ; k <= n ; k++ ) {
t = ddot( k-1, a[k], 1, b, 1 );
b[k] = ( b[k] - t ) / a[k][k];
}
/*
Now solve U * x = y.
*/
for ( k = n - 1 ; k >= 1 ; k-- ) {
b[k] = b[k] + ddot( n-k, a[k]+k, 1, b+k, 1 );
j = ipvt[k];
if ( j != k ) {
t = b[j];
b[j] = b[k];
b[k] = t;
}
}
return;
}

/*
Job = nonzero, solve Transpose(a) * x = b.

First solve Transpose(U) * y = b.
*/
for ( k = 1 ; k <= n - 1 ; k++ ) {
j = ipvt[k];
t = b[j];
if ( j != k ) {
b[j] = b[k];
b[k] = t;
}
daxpy( n-k, t, a[k]+k, 1, b+k, 1 );
}
/*
Now solve Transpose(L) * x = y.
*/
for ( k = n ; k >= 1 ; k-- ) {
b[k] = b[k] / a[k][k];
t = -b[k];
daxpy( k-1, t, a[k], 1, b, 1 );
}

}

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