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