Thread: Lu Decomposition Help Me Please

  1. #1
    Registered User
    Join Date
    Sep 2006
    Posts
    2

    Lu Decomposition Help Me Please

    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

  2. #2
    Devil's Advocate SlyMaelstrom's Avatar
    Join Date
    May 2004
    Location
    Out of scope
    Posts
    4,079
    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.
    Sent from my iPadŽ

  3. #3
    Registered User
    Join Date
    Sep 2006
    Posts
    2
    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");
                
    }

  4. #4
    Registered User
    Join Date
    Mar 2005
    Location
    Mountaintop, Pa
    Posts
    1,058
    Numerical recipes in C may be of some help to you. Especially, section 2.3

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. LU decomposition
    By khdani in forum C Programming
    Replies: 3
    Last Post: 10-15-2010, 08:08 PM
  2. LU decomposition
    By scottmanc in forum C++ Programming
    Replies: 1
    Last Post: 10-15-2003, 08:25 AM
  3. LUP Decomposition (simultaneous equations)
    By Markallen85 in forum C Programming
    Replies: 6
    Last Post: 08-24-2003, 02:08 AM