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

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

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

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

const int MAX = 2;

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[2][2] = {   {  2.0,   3.0 },
                          {  1.0,  4.0 }};

    luDecomposition(mat, 2);
    return 0;
}

輸出畫面

        Lower Triangular           Upper Triangular

1.0000 0.0000 2.0000 3.0000

0.5000 1.0000 0.0000 2.5000