so i have a c program lab with all my comments in there and now i need throw in some code for LU decomposition but i just can't figure it out, can anyone help, im gonna attach my lab along with a test file
Printable View
so i have a c program lab with all my comments in there and now i need throw in some code for LU decomposition but i just can't figure it out, can anyone help, im gonna attach my lab along with a test file
I can pretty much speak for every regular poster here when I say that it's preferred that you post your code (and only the important code) in the post rather than attach a source file which we have to download and open.
oppps here it is
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");
}
Numerical recipes in C may be of some help to you. Especially, section 2.3