home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
vis-ftp.cs.umass.edu
/
vis-ftp.cs.umass.edu.tar
/
vis-ftp.cs.umass.edu
/
pub
/
Software
/
ASCENDER
/
ascender.tar.Z
/
ascender.tar
/
Triangulate
/
gaussj_check.c
< prev
next >
Wrap
C/C++ Source or Header
|
1995-04-13
|
4KB
|
137 lines
/* @(#)gaussj.c 1.5 6/17/94 */
/* Gauss Jordan Elimination */
#include "cvar.h"
#include <math.h>
#include "alloc.h"
#include "error_message.s"
#define SWAP(a,b) {double temp = (a); (a)=(b); (b) = temp;}
FUNCTION_DEF ( int gaussj_check, (
double **a,
int n,
double *b))
{
/* Linear equation solution to Gauss-Jordan elimination.
a[1..n][1..n] is an input matrix of n by n elements.
b[1..n] is an input matrix of size n containing the
right-hand side vector.
On output, a is replaced by its matrix inverse, and b is replaced
by the corresponding solution vector.
*/
int *indxc, *indxr, *ipiv; /* Integer arrays used for book-keeping on the
* pivoting */
int i, icol, irow, j, k, l, ll;
double big, dum, pivinv;
/* Declare some temporary arrays */
indxc = VECTOR(1, n, int);
indxr = VECTOR(1, n, int);
ipiv = VECTOR(1, n, int);
for (j=1; j<=n; j++) ipiv[j] = 0;
for (i=1; i<=n; i++)
{
/* This is the main loop over the columns to be reduced */
big = 0.0;
for (j=1; j<=n; j++)
{
/* This is the outer loop of the search for the pivot element */
if (ipiv[j] != 1)
{
for (k=1; k<=n; k++)
{
if (ipiv[k] == 0)
{
if (fabs(a[j][k]) >= big)
{
big = fabs(a[j][k]);
irow = j;
icol = k;
}
}
else if (ipiv[k] > 1)
{
error_message("GAUSSJ : Singular Matrix - 1");
FREE (indxc, 1);
FREE (indxr, 1);
FREE (ipiv, 1);
return(FALSE);
}
}
}
}
++(ipiv[icol]);
/* We now have the pivot element, so we interchange rows, if
* needed, to put the pivot element on the diagonal. The columns
* are not physically interchanged, only relabelled : indxr[i],
* the colun of the i-th pivot element, is the i-th column, while
* indxr[i] is the row in which that pivot element was originally
* located. If indxr[i] != indxc[i] there is an implied column
* interchange. With this form of book-keeping, the solutions b
* will end up in the correct order and the inverse matrix will
* be scrambled by columns.
*/
if (irow != icol)
{
for (l=1; l<=n; l++) SWAP(a[irow][l], a[icol][l]);
SWAP(b[irow], b[icol]);
}
/* We are now ready to divide the pivot row by
* the pivot element, located ate irow, icol */
indxr[i] = irow;
indxc[i] = icol;
if (a[icol][icol] == 0.0)
{
error_message("GAUSSJ : Singular Matrix - 2");
FREE (indxc, 1);
FREE (indxr, 1);
FREE (ipiv, 1);
return(FALSE);
}
pivinv = 1.0 / a[icol][icol];
a[icol][icol] = 1.0;
for (l=1; l<=n; l++) a[icol][l] *= pivinv;
b[icol] *= pivinv;
for (ll=1; ll<=n; ll++)
{
/* Next we reduce the rows */
if (ll != icol)
{
/* .. except for the pivot one of course */
dum = a[ll][icol];
a[ll][icol] = 0.0;
if (dum != 0.0)
{
for (l=1; l<=n; l++) a[ll][l] -= a[icol][l]*dum;
b[ll] -= b[icol]*dum;
}
}
}
}
/* This is the end of the main loop over columns of the reduction.
* It only remains to unscramble the solution in view of the column
* interchanges. We do this by interchanging pairs of columns in the
* inverse order that the permutation was built up */
for (l=n; l>=1; l--)
{
if (indxr[l] != indxc[l])
for (k=1; k<=n; k++)
SWAP(a[k][indxr[l]], a[k][indxc[l]]);
}
/* Free temporaries */
FREE (indxc, 1);
FREE (indxr, 1);
FREE (ipiv, 1);
return(TRUE);
}