Code:
#include <stdio.h>
#include <stdlib.h>
/* MAX_SIZE is the maximum number of rows and columns.
Note that the matrices are stored one size larger because we want to
index them starting from 1, not 0, hence the extra constant MAX_SIZE_P1 */
#define MAX_SIZE 10
#define MAX_SIZE_P1 (MAX_SIZE + 1)
/* Function headers */
double msolve(int n,double a[MAX_SIZE_P1][MAX_SIZE_P1],double x[MAX_SIZE_P1],
double b[MAX_SIZE_P1]);
void read_input(FILE *f,double a[MAX_SIZE_P1][MAX_SIZE_P1],
double acopy[MAX_SIZE_P1][MAX_SIZE_P1],
double b[MAX_SIZE_P1],int *n);
void vprint(int n,double b[MAX_SIZE_P1]);
void atimesx(int n,double a[MAX_SIZE_P1][MAX_SIZE_P1],double x[MAX_SIZE_P1],
double bres[MAX_SIZE_P1]);
main()
{
FILE *f;
char fname[100];
int n,i,j;
double a[MAX_SIZE_P1][MAX_SIZE_P1],b[MAX_SIZE_P1],
x[MAX_SIZE_P1],acopy[MAX_SIZE_P1][MAX_SIZE_P1];
double bres[MAX_SIZE_P1];
printf("Enter input file name: ");
scanf("%s",fname);
f = fopen(fname,"r");
if (f == NULL) {
printf("Could not open data file: %s\n",fname);
exit(1);
}
read_input(f,a,acopy,b,&n); /* Read the A matrix, b vector, and a copy of A */
printf("Input b vector:\n");
vprint(n,b); /* Output the b vector */
printf("\nSolving system...");
msolve(n,a,x,b); /* Solve the system of equations */
printf("done. Solution is:\n");
vprint(n,x); /* Print the solution */
atimesx(n,acopy,x,bres); /* Substitute solution back into original system */
printf("\nSubstituting the result back into the original system gives:\n");
vprint(n,bres); /* Output the result from the substitution */
}
/* This function reads from the input file the a matrix and the b vector.
It also makes a copy of the a matrix for later use */
void read_input(FILE *f,double a[MAX_SIZE_P1][MAX_SIZE_P1],
double acopy[MAX_SIZE_P1][MAX_SIZE_P1],double b[MAX_SIZE_P1],int *n)
{
int i,j;
fscanf(f,"%d",n); /* Get the size of the matrix */
printf("Reading the A matrix and the b vector...");
for (i = 1; i <= *n; i++) {
for (j = 1; j <= *n; j++) {
fscanf(f,"%lf",&a[i][j]);
acopy[i][j] = a[i][j];
}
fscanf(f,"%lf",&b[i]);
}
printf("done.\n");
}
/* This function prints a vector b of size n */
void vprint(int n,double b[MAX_SIZE_P1])
{
int i;
for (i = 1; i <= n; i++)
printf("%lf\n",b[i]);
}
/* Performs basic Gaussian elimination */
double msolve(int n,double a[MAX_SIZE_P1][MAX_SIZE_P1],
double x[MAX_SIZE_P1],double b[MAX_SIZE_P1])
{
int k,i,j;
double factor,sum;
for (k = 1; k <= n-1; k++) {
for (i = k+1; i <= n; i++) {
factor = a[i][k]/a[k][k];
for (j = k+1; j <= n; j++) {
a[i][j] = a[i][j] - factor*a[k][j];
}
b[i] = b[i] - factor*b[k];
}
}
x[n] = b[n]/a[n][n];
for (i = n-1; i >= 1; i--) {
sum = b[i];
for (j = i + 1; j <= n; j++) {
sum = sum - a[i][j]*x[j];
}
x[i] = sum/a[i][i];
}
}
/* Multiplies A times x and stores the result in bres */
void atimesx(int n,double a[MAX_SIZE_P1][MAX_SIZE_P1],double x[MAX_SIZE_P1],
double bres[MAX_SIZE_P1])
{
int i,j;
for (i = 1; i <= n; i++) {
bres[i] = 0;
for (j = 1; j <= n; j++)
bres[i] = bres[i] + a[i][j]*x[j];
}
system("pause");
}