C語言 例題6-5 LU分解法 LU Decomposition

C語言 例題6-5 LU分解法 LU Decomposition

/* ex6-5.c is the LU decomposition of a matrix A
 * A=LU.
 */
#include <stdio.h>
#include <math.h>
void lluu(int nn,double a[5][5],double l[10][10],double u[10][10]);

void main()
{
    int n,nn,i,j;
    double u[10][10],l[10][10];
    nn=4;
    double a[5][5] = {{0, 0. , 0.  , 0. ,  0.},
                      {0, 8. , 1.  , 3. ,  2.},
                      {0, 2. , 9.  ,-1. , -2.},
                      {0, 1. , 3.  , 2. , -1.},
                      {0, 1. , 0.  , 6. ,  4.}};
    printf("The A Matrix.\n");
    for(i=1;i<=nn;i++)
    {
        for(j=1;j<=nn;j++)
        {
        //scanf("%lf",&a[i][j]);
        printf("%8.4lf ",a[i][j]);
        }
        printf("\n");
    }
    lluu(nn,a,l,u); /* call function lluu for L and U Matrix */
    printf("The Lower Triangular Matrix\n");
   
    for(i=1;i<=nn;i++)
    {
        for(j=1;j<=nn;j++)
            printf("%8.4lf ",l[i][j]);
        printf("\n");
    }
    printf("The Upper Triangular Matrix\n");
    for(i=1;i<=nn;i++)
    {
        for(j=1;j<=nn;j++)
        printf("%8.4lf ",u[i][j]);
        printf("\n");
    }
}


void lluu(int nn,double a[5][5],double l[10][10],double u[10][10])
{
    int n,i,j,k;
    double temp1,temp2;
    for(i=1;i<=nn;i++)
        for(j=1;j<=nn;j++)
        {
        u[i][j]=0.0;
        l[i][j]=0.0;
        }
    for(n=1;n<=nn;n++)
    {
        for(j=n;j<=nn;j++)
        {
        temp1=0.0;
        for(k=1;k<=n-1;k++)
            temp1=temp1+l[n][k]*u[k][j];
         u[n][j]=a[n][j]-temp1;
        }
        for(i=n+1;i<=nn;i++)
        {
        temp2=0.0;
        for(k=1;k<=n-1;k++)
            temp2=temp2+(l[i][k]*u[k][n]);
        l[i][n]=(a[i][n]-temp2)/u[n][n];
        }
    }
   
    for(i=1;i<=nn;i++)
        l[i][i]=1.0;
    return;
}


輸出畫面
The A Matrix.
  8.0000   1.0000   3.0000   2.0000
  2.0000   9.0000  -1.0000  -2.0000
  1.0000   3.0000   2.0000  -1.0000
  1.0000   0.0000   6.0000   4.0000
The Lower Triangular Matrix
  1.0000   0.0000   0.0000   0.0000
  0.2500   1.0000   0.0000   0.0000
  0.1250   0.3286   1.0000   0.0000
  0.1250  -0.0143   2.5455   1.0000
The Upper Triangular Matrix
  8.0000   1.0000   3.0000   2.0000
  0.0000   8.7500  -1.7500  -2.5000
  0.0000   0.0000   2.2000  -0.4286
  0.0000   0.0000   0.0000   4.8052