/*****************************************************
Expanded Name: gauss-jordan fortran style
Global Function List: gauss_jordan
Portability: Standard C
*****************************************************/
/* Standard Library */
#include <float.h>
#include <math.h>
#include <stdlib.h>
/* Types and Prototypes */
#include <matrix_ t.h>
/* Own */
#include <gs_jrd.h>
/*****************************************************
Name: gauss_jordan
Parameters: A matrix of coeficients
n number of rows, cols in A
B matrix of constant vectors
m number of vectors, cols, in B
Return: SUCCESS, 0, or FAIL, -1, directly
Inverse of A in A indirectly
Solution vectors X in B indirectly
Description: Gauss-Jordan elimination with partial
pivoting is used to invert a matrix and
solve a system of linear equations.
The equations are expressed in the matrix
form as [A] [X] = [B]. The A matrix is an
n by n matrix. The B matrix is an n by m
matrix. B is replaced with the solution
vectors X. A is replaced with the
inverse of A.
*****************************************************/
int gauss_jordan( matrix_t A, size_t n,
matrix_t B, size_t m )
{
size_t i, j, k, PivotRow;
double Pivot, *Ai, *Aj;
/* Loop through all the rows. */
for ( i = 0; i < n; i++ )
{
/* Set the default pivot to the
** current row. */
PivotRow = i;
Pivot = A[i] [i];
/* Loop through the rows below the current
** row Looking for a large diagonal to
** pivot */
for ( j = i + 1; j < n; j++ )
{
/* Find the pivot row */
if ( fabs( A[j][i] ) >= fabs( Pivot ) )
{
PivotRow = j;
Pivot= A[j] [i];
} /* if */
} /* for j */
/* Check to make sure we are not pivoting
** around the current row */
if ( PivotRow != i )
{
/* Swap the current row with the
** pivot row */
/* Swap the pointers to the rows. */
Ai = A[i];
A[i] = A[PivotRow];
A[PivotRow] = Ai;
/* Swap the pointer to the rows. */
Ai = B[i];
B[i] = B[PivotRow];
B[PivotRow] = Ai;
} /* if PivotRow */
/* Check for a divide by zero error */
if ( fabs( Pivot ) <= DBL_MIN )
{
return ( -1 );
} /* if */
/* Invert the pivot so we can use * operator
** instead of / operator */
Pivot = 1.0 / Pivot;
/* Set the diagonal to 1.0 */
Ai = A[i]; Ai[i] = 1.0;
/* Divide by the pivot */
for ( j = 0; j < n; j++ )
{
Ai[j] *= Pivot;
} /* for j */
Ai = B[i];
/* Divide by the pivot */
for ( j = 0; j < m; j++ )
{
Ai[j] *= Pivot;
} /* for j */
/* Loop through all the rows doing row
** reduction */
for ( j = 0; j < n; j++ )
{
if ( j != i )
{
/* Reuse the Pivot variable as a
** dummy */
Ai = A[i]; Aj = A[j]; Pivot = Aj[i];
Aj[i] = 0.0;
/* Row reduce the matrix */
for ( k = 0; k < n; k++ )
{
Aj[k] -= Ai[k] * Pivot;
} /* for k */
Ai = B[i]; Aj = B[j];
/* Row reduce the solution vectors */
for ( k = 0; k < m; k++ )
{
Aj[k] -= Ai[k] * Pivot;
} /* for k */
} /* if j */
} /* for j */
} /* for i */
return ( 0 );
} /* function gauss_jordan */
/* End of File */