This is gauss.c in view mode; [Download] [Up]
#include "lpc.h"
gauss(aold, bold, b)
double aold[POLEMAX][POLEMAX], *bold, *b;
{
double amax, dum, pivot;
double c[POLEMAX], a[POLEMAX][POLEMAX];
int i, j, k, l, istar, ii, lp;
/* aold and bold untouched by this subroutine */
for (i=0; i < NPOLE ;++i) {
c[i] = bold[i];
for (j=0; j < NPOLE ;++j)
a[i][j] = aold[i][j];
}
/* eliminate i-th unknown */
for (i=0; i < NPOLE - 1 ;++i) { /* find largest pivot */
amax = 0.0;
for (ii=i; ii < NPOLE ;++ii) {
if (fabs(a[ii][i]) >= amax) {
istar = ii;
amax = fabs(a[ii][i]);
}
}
if (amax < 1e-20) {
fprintf(stderr, "gauss: ill-conditioned\n");
return(0); /* ill-conditioned */
}
for (j=0; j < NPOLE ;++j) { /* switch rows */
dum = a[istar][j];
a[istar][j] = a[i][j];
a[i][j] = dum;
}
dum = c[istar];
c[istar] = c[i];
c[i] = dum;
/* pivot */
for (j=i+1; j < NPOLE ;++j) {
pivot = a[j][i] / a[i][i];
c[j] = c[j] - pivot * c[i];
for (k=0; k < NPOLE ;++k)
a[j][k] = a[j][k] - pivot * a[i][k];
}
} /* return if last pivot is too small */
if (fabs(a[NPOLE-1][NPOLE-1]) < 1e-20) {
fprintf(stderr, "gauss: ill-conditioned.\n");
return(0); /* ill-conditioned */
}
*(b+NPOLE-1) = c[NPOLE-1] / a[NPOLE-1][NPOLE-1];
/* back substitute */
for (k=0; k<NPOLE-1; ++k) {
l = NPOLE-1 -(k+1);
*(b+l) = c[l];
lp = l + 1;
for (j = lp; j<NPOLE; ++j)
*(b+l) += -a[l][j] * *(b+j);
*(b+l) /= a[l][l];
}
return(1);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.