C語言 例題6-4 LU分解法

C語言 例題6-4 LU分解法

A  =    {   { -6.0, 14,  2 },
                {  4,     4,   9 },
                { -3,     2,   13} };


// C Program to decompose a matrix into
// lower and upper traingular matrix

#include <stdio.h>
#include <stdlib.h>

const int MAX = 3;
 
void luDecomposition(float mat[][MAX], int n)
{
    float lower[n][n], upper[n][n];
    for (int i=0;i<=n;i++)
    {
        for (int j=0;j<=n;j++)
        {   
            lower[i][j]=0.0;
            upper[i][j]=0.0;
        }   
    }
    //memset(lower, 0, sizeof(lower));
    //memset(upper, 0, sizeof(upper));
 
    // Decomposing matrix into Upper and Lower
    // triangular matrix
    for (int i = 0; i < n; i++) {
 
        // Upper Triangular
        for (int k = i; k < n; k++) {
 
            // Summation of L(i, j) * U(j, k)
            float sum = 0;
            for (int j = 0; j < i; j++)
                sum += (lower[i][j] * upper[j][k]);
 
            // Evaluating U(i, k)
            upper[i][k] = mat[i][k] - sum;
        }
 
        // Lower Triangular
        for (int k = i; k < n; k++) {
            if (i == k)
                lower[i][i] = 1; // Diagonal as 1
            else {
 
                // Summation of L(k, j) * U(j, i)
                float sum = 0;
                for (int j = 0; j < i; j++)
                    sum += (lower[k][j] * upper[j][i]);
 
                // Evaluating L(k, i)
                lower[k][i] = (mat[k][i] - sum) / upper[i][i];
            }
        }
    }
 
    // setw is for displaying nicely
    printf("         Lower Triangular                       Upper Triangular\n");
 
    // Displaying the result :
    for (int i = 0; i < n; i++) {
        // Lower
        for (int j = 0; j < n; j++)
            printf("\t%.4lf\t",  lower[i][j]); 
        printf("\t");
 
        // Upper
        for (int j = 0; j < n; j++)
            printf("\t%.4lf\t",  upper[i][j]);
        printf("\n"); 
    }
}
 
// Driver code
int main()
{
    float mat[3][3] = {   { -6.0, 14, 2 },
                          {  4,  4, 9 },
                          { -3,  2, 13} };
 
    luDecomposition(mat, 3);
    return 0;
}

輸出畫面

         Lower Triangular                       Upper Triangular
1.0000 0.0000 0.0000 -6.0000 14.0000 2.0000
-0.6667 1.0000 0.0000 0.0000 13.3333 10.3333
0.5000 -0.3750 1.0000 0.0000 0.0000 15.8750