[net.sources] Matrix solution with Gaussian elimination and LU decomposition

ken@turtleva.UUCP (Ken Turkowski) (08/11/83)

echo x - matrix_solve
mkdir matrix_solve
echo x - matrix_solve/READ_ME
cat >matrix_solve/READ_ME <<'!Funky!Stuff!'
Several weeks ago, umcp-cs!fred sent out a few matrix manipulation
routines.  These are at least theoretically correct, and work for
"nice" matrices.  However, they do involve more computation than
necessary, and may not be so stable when the matrices are not so
well-conditioned.

In this packet, I include two routines which are more robust for
solving simultaneous equations and matrix inversion.  Straightforward
solution by Gaussian elimination can be found in gausselim.c, and
solution with LU decomposition (more suitable for matrix inversion and
multiple equations with the same coefficient matrix) can be found in
lu.c.

			    Ken Turkowski
			CADLINC, Palo Alto
		{decvax,ucbvax}!decwrl!turtlevax!ken
		{allegra, ihnp4}!amd70!turtlevax!ken
!Funky!Stuff!
echo x - matrix_solve/gausselim.c
cat >matrix_solve/gausselim.c <<'!Funky!Stuff!'
/* Gausselim() solves the nth order linear equation Ax=b, using Gaussian
 * elimination with partial pivoting and row equilibration.
 * The arguments are:
 *
 *	A - the (n x n) coefficient matrix
 *	x - the (1 x n) solution vector
 *	b - the (1 x n) constant vector
 *	n - the order of the equation
 *
 *  1 is returned if a solution was obtained,
 *  and 0 if either the matrix A is singular.
 *
 *	Ken Turkowski, CADLINC, Inc. (turtlevax!ken)
 *	written 2/16/79, revised and enhanced 8/9/83.
 */

# define MAXN 15	/* Size of working storage: n <= MAXN */
# define Ael(i,j)	A[(i*n)+j]

gausselim(A, x, b, n)
double *A;	/* n x n */
double *x, *b;	/* 1 x n */
register int n;
{
	register int i, j, k;
	int tempi;
	int ps[MAXN];			/* Pivoting sequence */
	double AB[MAXN][MAXN+1];	/* Working storage: augmented matrix */
	double mult, biggest, tempf;
	double sum;
	double fabs();

	/* Make augmented matrix AB */
	for (i = 0; i < n; i++) {	/* For each row ... */
	    /* Find the largest element in each row */
	    biggest = 0.0;
	    for (j = 0; j <= n; j++)
		if (biggest < (tempf = fabs(Ael(i,j))))
		    biggest = tempf;

	    if (biggest == 0.0)
		return(0);	/* Zero row: singular matrix */

	    /* Normalize each row by its largest element */
	    mult = 1.0 / biggest;
	    for (j = 0; j < n; j++)
		AB[i][j] = Ael(i,j) * mult;
	    AB[i][n] = b[i] * mult;

	    /* Initialize the pivoting sequence */
	    ps[i] = i;
	}

	/* Gaussian elimination with partial pivoting */
	for (i = 0; i < n-1; i++) {	/* For each column ... */
	    /* Choose a pivot from the largest number in the remaining rows */
	    biggest = 0.0;
	    for (j = i; j < n; j++) {
		if (biggest < (tempf = AB[ps[j]][i])) {
		    biggest = tempf;
		    k = j;
		}
	    }

	    if (biggest == 0.0)	/* Oops, no nonzero pivot */
		return(0);	/* Gaussian elimination failed */

	    if (k != i) {	/* Update the pivot sequence */
		tempi = ps[i];
		ps[i] = ps[k];
		ps[k] = tempi;
	    }

	    /* Eliminate column ps[i] from the remaining equations */
	    for (j = i+1; j < n; j++) {
		mult = AB[ps[j]][i] / AB[ps[i]][i];
		AB[ps[j]][i] = 0.0;
		for (k = i+1; k <= n; k++)
		    AB[ps[j]][k] -= mult * AB[ps[i]][k];
	    }
	}

	if (AB[ps[n-1]][n-1] == 0.0)
	    return(0);		/* Singular matrix */

	/* Back substitution */
	for (i = n-1; i >= 0; i--) {
	    sum = AB[ps[i]][n];
	    for (j = i + 1; j < n; j++)
		sum -= AB[ps[i]][j] * x[j];
	    x[i] = sum / AB[ps[i]][i];
	}

	return(1);
}
!Funky!Stuff!
echo x - matrix_solve/lu.c
cat >matrix_solve/lu.c <<'!Funky!Stuff!'
/* This module solves linear equations in several variables (Ax = b) using
 * LU decomposition with partial pivoting and row equilibration.  Although
 * slightly more work than Gaussian elimination, it is faster for solving
 * several equations using the same coefficient matrix.  It is
 * particularly useful for matrix inversion, by sequentially solving the
 * equations with the columns of the unit matrix.
 *
 * lu_decompose() decomposes the coefficient matrix into the LU matrix,
 * and lu_solve() solves the series of matrix equations using the
 * previous LU decomposition.
 *
 *	Ken Turkowski, CADLINC, Inc. (turtlevax!ken)
 *	written 3/2/79, revised and enhanced 8/9/83.
 */

# define MAXN 16	/* Size of working storage: n <= MAXN */

/* lu_decompose() decomposes the coefficient matrix A into upper and lower
 * triangular matrices, the composite being the LU matrix.
 *
 * The arguments are:
 *
 *	a - the (n x n) coefficient matrix
 *	n - the order of the matrix
 *
 *  1 is returned if the decomposition was successful,
 *  and 0 is returned if the coefficient matrix is singular.
 */

# define ael(i,j)	a[i*n+j]

static int ps[MAXN];
static double lu[MAXN][MAXN];

lu_decompose(a, n)
register double *a;
register int n;
{
	double scales[MAXN];
	register int i, j, k;
	int pivotindex;
	double pivot, biggest, mult, tempf;
	double fabs();

	for (i = 0; i < n; i++) {	/* For each row */
	    /* Find the largest element in each row for row equilibration */
	    biggest = 0.0;
	    for (j = 0; j < n; j++)
		if (biggest < (tempf = fabs(lu[i][j] = ael(i,j))))
		    biggest  = tempf;
	    if (biggest != 0.0)
		scales[i] = 1.0 / biggest;
	    else {
		scales[i] = 0.0;
		return(0);	/* Zero row: singular matrix */
	    }

	    ps[i] = i;		/* Initialize pivot sequence */
	}

	for (k = 0; k < n-1; k++) {	/* For each column */
	    /* Find the largest element in each column to pivot around */
	    biggest = 0.0;
	    for (i = k; i < n; i++) {
		if (biggest < (tempf = fabs(lu[ps[i]][k]) * scales[ps[i]])) {
		    biggest = tempf;
		    pivotindex = i;
		}
	    }
	    if (biggest == 0.0)
		return(0);	/* Zero column: singular matrix */
	    if (pivotindex != k) {	/* Update pivot sequence */
		j = ps[k];
		ps[k] = ps[pivotindex];
		ps[pivotindex] = j;
	    }

	    /* Pivot, eliminating an extra variable  each time */
	    pivot = lu[ps[k]][k];
	    for (i = k+1; i < n; i++) {
		lu[ps[i]][k] = mult = lu[ps[i]][k] / pivot;
		if (mult != 0.0) {
		    for (j = k+1; j < n; j++)
			lu[ps[i]][j] -= mult * lu[ps[k]][j];
		}
	    }
	}

	if (lu[ps[n-1]][n-1] == 0.0)
	    return(0);		/* Singular matrix */
}

/* lu_solve() solves the linear equation (Ax = b) after the matrix A has
 * been decomposed with lu_decompose() into the lower and upper triangular
 * matrices L and U.
 *
 * The arguments are:
 *
 *	x - the solution vector
 *	b - the constant vector
 *	n - the order of the equation
*/

lu_solve(x, b, n)
register double *x, *b;
{
	register int i, j;
	double dot;

	/* Vector reduction using U triangular matrix */
	for (i = 0; i < n; i++) {
	    dot = 0.0;
	    for (j = 0; j < i; j++)
		dot += lu[ps[i]][j] * x[j];
	    x[i] = b[ps[i]] - dot;
	}

	/* Back substitution, in L triangular matrix */
	for (i = n-1; i >= 0; i--) {
	    dot = 0.0;
	    for (j = i+1; j < n; j++)
		dot += lu[ps[i]][j] * x[j];
	    x[i] = (x[i] - dot) / lu[ps[i]][i];
	}
}
!Funky!Stuff!