The C language Realization __c language of Lu decomposition method

The product of the transformation of coefficient matrix A into equivalence two matrices L and U
, where L and u are respectively the lower triangular matrix and the upper triangular matrix. When the master of all orders of A is not 0 o'clock, matrix A can be decomposed into a=lu (all order master is not 0, the matrix is not necessarily not be able to do LU decomposition). Where l is the lower triangular matrix, and U is the upper triangular matrix.

In essence, Lu decomposition is an expression of Gauss elimination method. Essentially, a transformation matrix is transformed into an upper triangular matrix through elementary row transformations, and the transform matrices are the lower triangular matrices of a unit. This is the so-called Durrit algorithm (Doolittle algorithm): From the bottom to the matrix A primary row transformation, the diagonal left element into 0, and then prove that the effect of these row transformations is equivalent to the left by a series of units under the triangular matrix, The inverse of the product of the triangle matrix in this series of units is the L matrix, which is also a unit lower triangular matrix.

–

Construction Matrix A,b

                {1,2,3,4,}
    a=[         {1,4,2,-8,}
                {1,-1,4,1}   ]
                {1,3,5,2}

b = [14,-17,2,8]^t

Code implementation:

#include <stdio.h> #include <stdlib.h>//lu decomposition method for solving linear equations//copyright @ Mryang double Sumu (double l[4][4], DOUBL
    e U[4][4], int i, int j) {Double SU = 0.0;
    for (int k = 1; k <= i-1; k++) {SU + = l[i-1][k-1] * U[k-1][j-1];
return SU;
    }//calculates the sum 1 double suml (double l[4][4], double u[4][4], int i, int j) {Double SL = 0.0;
    for (int k = 0; k <= j-1; k++) {SL = l[i-1][k-1] * U[k-1][j-1];

return SL;
    }//calculates the sum 2 double SumY (double l[4][4], double y[4],int i) {double sy=0.0;
    for (int k = 1; k <= i-1; k++) {SY + l[i-1][k-1] * y[k-1];
return SY;
    }//calculates the sum 3 double sumx (double u[4][4], double x[4],int i, int m) {double SX = 0.0;
    for (int k = i+1 k <= m; k++) {SX = u[i-1][k-1] * X[k-1];
return SX; }//calculates sum 4 int main () {double A[4][4] = {{1,2,3,1,}, {1,4,1,-1,}, {1,-1,-2,3,} , {1,3,-1,2}};The coefficients are stored in a two-dimensional array double l[4][4] = {0};
    Double U[4][4] = {0};//initialization part double b[4] = {8,8,12,19};
    int n = 4;//n order//Output [Ab] printf ("[a]:\n"); for (int i = 1; I <= n; i++) {for (int j = 1; J <= N; j) {printf ("%f\t", a[i-1][
        J-1]);
    printf ("\ n");  ///Compute l,u for (int i = 1; I <= n; i++) {l[i-1][i-1] = 1;//diagonal element is 1 for (int j = i; J <= N

            J + +) {//due to array subscript starting from 0 so i-1,j-1 u[i-1][j-1] = a[i-1][j-1]-Sumu (L,U,I,J); if (j+1 <= n) l[j][i-1] = (A[j][i-1]-suml (l,u,j+1,i))/u[i-1][i-1];//i to j+1,j variable i}}//Output U pri
    NTF ("u:\n"); for (int i = 1; I <= n; i++) {for (int j = 1; J <= N; j) {printf ("%f\t", u[i-1][j
        -1]);
    printf ("\ n");
    //Output L printf ("l:\n");
    for (int i = 1; I <= n; i++) {for (int j = 1; J <= N; j +) {        printf ("%f\t", l[i-1][j-1]);
    printf ("\ n");
    ///by Ly=b y double y[4] = {0.0};

    Y[0] = b[0];//y (1) = B (1);
    for (int i = 2; I <= n; i++) {y[i-1] = b[i-1]-SumY (l,y,i);

    }///Ux=y x double x[4] = {0.0};
    for (int i = n; I >= 1; i--) {x[i-1] = (Y[i-1]-sumx (u,x,i,n))/u[i-1][i-1];
    //Output y printf ("y:\n");
    for (int i = 0; i < n; i++) {printf ("%f\n", Y[i]);
    printf ("\ n");
    Output x printf ("x:\n");
    for (int i = 0; i < n; i++) {printf ("%f\n", X[i]);
    printf ("\ n");
    System ("pause");
return 0;
 }

Run Result:

So
X1=1
X2= 2
X3=-1
x4 = 3