Logic error for Gauss elimination

Question

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...

Answer 1

This probably doesn't answer your question in the way you intended, but programming your own numerically-stable matrix algorithms is about as well-advised as do-it-yourself surgery.

There's a very nice library called TNT/JAMA from a reputable source (NIST) which does elementary matrix math in C++. To solve Ax=B, first factor A (the QR decomposition is a good general method, you can use LU but it's less numerically stable), then call solve(B). This works both for square matrices, where it's exact (subject to numerical computation issues), and overdetermined systems, where you get a least-squares answer.

Answer 2

//eliminate the coefficients of X(I)
for (i = k; i < n; i++)

Should this maybe be

for (i = k + 1; i < n; i++)

The way it is now, I believe this will cause you to divide the pivot row by itself, zeroing it out.

Answer 3

why the answer is The solution is x(1) = 3.000000 x(2) = -6.142857 x(3) = 0.167910 Determinant = -102.000000

can anyone help me why the answer of x1,x2,and x3 is not same with the correct answer.