Code:

#include <stdio.h>
#include <time.h>
/* Global variable declarations.*/
double a[5000][5000]; /* The two-dimensional matrix A. */
double b[5000]; /* The right-hand side vector b. */
double x[5000]; /* The solution vector x. */
int n; /* The dimension (n x n) of the matrix system. */
/* This is the main program. */
int main(void)
{
/* Local variable declarations. */
char fileaname[100]; /* The name of the input file containing the matrix A. */
char filebname[100]; /* The name of the input file containing the vector b. */
char filexname[100]; /* The name of the input file containing the vector x. */
int gausselimination(void); /* Declare our gausselimination function. */
char line[100]; /* The input line from stdin. */
clock_t time1; /* The current time before the call to gausselimination. */
clock_t time2; /* The current time after the call to gausselimination. */
FILE *filea; /* The input file containing the matrix A. */
FILE *fileb; /* The input file containing the vector b. */
FILE *filex; /* The input file containing the vector x. */
int i; /* Loop index. */
int j; /* Loop index. */
int ntemp; /* A temporary n that we'll use to check the input files. */
/* Get the file name for the matrix A. */
printf("\n");
printf("Enter the name of the file containing the matrix A: ");
fgets(line, sizeof(line), stdin);
sscanf(line, "%s", &fileaname);
/* Open the file for the matrix A. If the file cannot be open, then send
* an error message to the user. */
filea = fopen(fileaname, "r");
if(filea == NULL)
{
printf("ERROR: Cannot open %s.\n", fileaname);
return(8);
}
printf("Reading the matrix A ... ");
/* Get the dimension n of the matrix A. */
fgets(line, sizeof(line), filea);
sscanf(line, "%d", &ntemp);
n = ntemp;
/* Loop over the rows using the index i. Then loop over the columns using
* the index j. At each row, read the entry a[i][j]. */
for(i = 0; i <= n-1; ++i)
{
for(j = 0; j <= n-1; ++j)
{
fgets(line, sizeof(line), filea);
sscanf(line, "%lf", &a[i][j]);
}
}
/* Close the file. */
fclose(filea);
printf("Done.\n");
/* Get the file name for the vector b. */
printf("\n");
printf("Enter the name of the file containing the vector b: ");
fgets(line, sizeof(line), stdin);
sscanf(line, "%s", &filebname);
/* Open the file for the vector b. If the file cannot be open, then send
* an error message to the user. */
fileb = fopen(filebname, "r");
if(fileb == NULL)
{
printf("ERROR: Cannot open %s.\n", fileaname);
return(8);
}
printf("Reading the vector b ... ");
/* Read the dimension n of the vector b. If this n does not match the n
* that we read above, then alert the user. */
fgets(line, sizeof(line), fileb);
sscanf(line, "%d", &ntemp);
if(ntemp != n)
{
printf("ERROR: The dimension of b does not match the dimension of A.\n");
return(0);
}
n = ntemp;
/* Loop over the rows using the index i. At each row, read the entry b[i]. */
for(i = 0; i <= n-1; ++i)
{
fgets(line, sizeof(line), fileb);
sscanf(line, "%lf", &b[i]);
}
/* Close the file. */
fclose(fileb);
printf("Done.\n");
printf("\n");
printf("Performing the Gauss elimination ... ");
/* We are going to time how long it takes to perform the Gauss elimination.
* The variable time1 contains the current time before the call to the subroutine. */
time1 = clock();
/* This is the call to the Gauss elimination function. */
gausselimination();
/* The variable time2 contains the current time after the call to the subroutine. */
time2 = clock();
/* Compute how long it took to perform the Gauss elimination, and print the answer
* to the user. */
printf("Done.\n");
printf("The elapsed CPU time was %.8f seconds.\n", (float) (time2-time1)/CLOCKS_PER_SEC);
/* Get the file name for the vector x. */
printf("\n");
printf("Enter the name of the file to write the vector x: ");
fgets(line, sizeof(line), stdin);
sscanf(line, "%s", &filexname);
/* Open the file for the vector x. If the file exists already, then it will be
* written over. */
filex = fopen(filexname, "w");
/* Print the dimension n to the file containing the vector x. */
sprintf(line, "%d\n", n);
fputs(line, filex);
/* Loop over the rows using the index i. At each row, print the entry x[i]
* to the file. */
for(i = 0; i <= n-1; ++i)
{
sprintf(line, "%.14lf\n", x[i]);
fputs(line, filex);
}
/* Close the file. */
fclose(filex);
return(0);
}
/* This function does the Gauss elimination on A and b and returns
* the solution vector x. */
int gausselimination(void)
{
/* Local variable declarations. */
float elementbelowpivot; /* The elements below the pivot element but in the same column. */
float pivot; /* The value of the pivot element. */
float tempsum; /* A temporary sum for the backward substitution. */
int i; /* Loop index for rows. */
int ii; /* Loop index for rows. */
int j; /* Loop index for columns. */
/* First triangularize the coefficient matrix A.
* Loop over all of the rows. */
for(i = 0; i <= n-1; ++i)
{
/* Divide every element in the row by the pivot element.
* First save the pivot element. */
pivot = a[i][i];
/* Then loop over all of the columns after the pivot element.*/
for(j = i; j <= n-1; ++j)
{
/* Then divide each element by the pivot element. */
a[i][j] = a[i][j] / pivot;
}
/* Divide the element in the vector b by the pivot element. */
b[i] = b[i] / pivot;
/* Loop over the remaining rows. Multiply the current row
* by each element in the column below the pivot element,
* and then subtract the current row from each remaining row.
* We are zero-ing out the lower triangular portion of the
* matrix A. First loop over all of the rows below the
* current row. */
for(ii = i+1; ii <= n-1; ++ii)
{
/* Save the element below the pivot but in this row. */
elementbelowpivot = a[ii][i];
/* Then loop over the columns, starting at the column that
* contains the pivot element. */
for(j = i; j <= n-1; ++j)
{
/* Multiply the pivot row by the element below the pivot,
* and then subtract that quantity from the element in
* the row below the pivot. */
a[ii][j] = a[ii][j] - a[i][j] * elementbelowpivot;
}
/* Don't forget to do the same thing to the vector b. */
b[ii] = b[ii] - b[i] * elementbelowpivot;
}
}
/* The coefficient matrix A has been triangularized. Now it should be
* easy to use back substitution to solve for the vector x.
* Loop over the rows, starting at the bottom and going up. */
for(i = n-1; i >= 0; --i)
{
/* Each element in x is equal to the corresponding element in b
* minus the sum of the products of any a's and x's below the current
* x. (See page 2.10 in Dr. Westerink's lecture notes for a better
* description of backward substitution.) So we need to sum the
* products of the coefficients a and any already-solved-for x.
* Zero out our temporary sum. */
tempsum = 0.0;
/* Loop over any columns/rows to the right/below the current x. */
for(j = i+1; j <= n-1; ++j)
{
tempsum = tempsum + a[i][j] * x[j];
}
/* Solve for the current x. */
x[i] = b[i] - tempsum;
}
return(0);
}