Logic error problem with the Gaussian Elimination code...This code was from my Numerical Methods text in 1990's. The code is typed in from the book- not producing correct output...
Sample Run:
SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS
USING GAUSSIAN ELIMINATION
This program uses Gaussian Elimination to solve the
system Ax = B, where A is the matrix of known
coefficients, B is the vector of known constants
and x is the column matrix of the unknowns.
Number of equations: 3
Enter elements of matrix [A]
A(1,1) = 0
A(1,2) = -6
A(1,3) = 9
A(2,1) = 7
A(2,2) = 0
A(2,3) = -5
A(3,1) = 5
A(3,2) = -8
A(3,3) = 6
Enter elements of [b] vector
B(1) = -3
B(2) = 3
B(3) = -4
SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS
The solution is
x(1) = 0.000000
x(2) = -1.#IND00
x(3) = -1.#IND00
Determinant = -1.#IND00
Press any key to continue . . .
The code as copied from the text...
//Modified Code from C Numerical Methods Text- June 2009
#include <stdio.h>
#include <math.h>
#define MAXSIZE 20
//function prototype
int gauss (double a[][MAXSIZE], double b[], int n, double *det);
int main(void)
{
double a[MAXSIZE][MAXSIZE], b[MAXSIZE], det;
int i, j, n, retval;
printf("\n \t SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS");
printf("\n \t USING GAUSSIAN ELIMINATION \n");
printf("\n This program uses Gaussian Elimination to solve the");
printf("\n system Ax = B, where A is the matrix of known");
printf("\n coefficients, B is the vector of known constants");
printf("\n and x is the column matrix of the unknowns.");
//get number of equations
n = 0;
while(n <= 0 || n > MAXSIZE)
{
printf("\n Number of equations: ");
scanf ("%d", &n);
}
//read matrix A
printf("\n Enter elements of matrix [A]\n");
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
{
printf(" A(%d,%d) = ", i + 1, j + 1);
scanf("%lf", &a[i][j]);
}
//read {B} vector
printf("\n Enter elements of [b] vector\n");
for (i = 0; i < n; i++)
{
printf(" B(%d) = ", i + 1);
scanf("%lf", &b[i]);
}
//call Gauss elimination function
retval = gauss(a, b, n, &det);
//print results
if (retval == 0)
{
printf("\n\t SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS\n");
printf("\n\t The solution is");
for (i = 0; i < n; i++)
printf("\n \t x(%d) = %lf", i + 1, b[i]);
printf("\n \t Determinant = %lf \n", det);
}
else
printf("\n \t SINGULAR MATRIX \n");
return 0;
}
/* Solves the system of equations [A]{x} = {B} using */
/* the Gaussian elimination method with partial pivoting. */
/* Parameters: */
/* n - number of equations */
/* a[n][n] - coefficient matrix */
/* b[n] - right-hand side vector */
/* *det - determinant of [A] */
int gauss (double a[][MAXSIZE], double b[], int n, double *det)
{
double tol, temp, mult;
int npivot, i, j, l, k, flag;
//initialization
*det = 1.0;
tol = 1e-30; //initial tolerance value
npivot = 0;
//mult = 0;
//forward elimination
for (k = 0; k < n; k++)
{
//search for max coefficient in pivot row- a[k][k] pivot element
for (i = k + 1; i < n; i++)
{
if (fabs(a[i][k]) > fabs(a[k][k]))
{
//interchange row with maxium element with pivot row
npivot++;
for (l = 0; l < n; l++)
{
temp = a[i][l];
a[i][l] = a[k][l];
a[k][l] = temp;
}
temp = b[i];
b[i] = b[k];
b[k] = temp;
}
}
//test for singularity
if (fabs(a[k][k]) < tol)
{
//matrix is singular- terminate
flag = 1;
return flag;
}
//compute determinant- the product of the pivot elements
*det = *det * a[k][k];
//eliminate the coefficients of X(I)
for (i = k; i < n; i++)
{
mult = a[i][k] / a[k][k];
b[i] = b[i] - b[k] * mult; //compute constants
for (j = k; j < n; j++) //compute coefficients
a[i][j] = a[i][j] - a[k][j] * mult;
}
}
//adjust the sign of the determinant
if(npivot % 2 == 1)
*det = *det * (-1.0);
//backsubstitution
b[n] = b[n] / a[n][n];
for(i = n - 1; i > 1; i--)
{
for(j = n; j > i + 1; j--)
b[i] = b[i] - a[i][j] * b[j];
b[i] = b[i] / a[i - 1][i];
}
flag = 0;
return flag;
}
The solution should be: 1.058824, 1.823529, 0.882353 with det as -102.000000
Any insight is appreciated...
