矩阵Crout分解的C语言实现

//A = LU，其中L为下三角阵，U为单位上三角阵

#include <stdio.h>
#include <math.h>

void print_matrix(int row, int col, double matrix[row][col])
{
    for(int i = 0;i < row; i++)
    {
        for(int j = 0;j < col; j++)
        {
            printf("%lf\t",matrix[i][j]);
        }    
        printf("\n");
    }
    printf("----------------------------------------------");
    printf("\n");
}

int main()
{
    int n = 3;
    double A[3][3] = {{2, 1, 1}, {1, 3, 2}, {1, 2, 2}};//A[n][n]
    double U[n][n];
    for (int j = 0; j < n; j++)//initialize U = I
    {
        for (int k = 0; k < n; k++)
        {
            if (k == j)
            {
                U[j][k] = 1;
            }
            else
            {
                U[j][k] = 0;
            }
        }        
    }
    double L[n][n];
    for (int i = 0; i < n; i++)//initialize L = 0
    {
        for (int j = 0; j < n; j++)
        {
            L[i][j] = 0;
        }        
    }

    //get L, U
    for (int k = 0; k < n; k++)
    {
        if (k == 0)
        {            
            for (int i = k; i < n; i++)
            {
                L[i][k] = A[i][k];
            }
            for (int j = k+1; j < n; j++)
            {         
                U[k][j] = A[k][j]/L[k][k];
            }
        }
        else
        {            
            for (int i = k; i < n; i++)
            {
                double sum2 = 0;
                for (int r = 0; r <= k-1; r++)
                {
                    sum2 = sum2 + L[i][r]*U[r][k];
                }            
                L[i][k] = A[i][k] - sum2;
            }
            for (int j = k+1; j < n; j++)
            {
                double sum1 = 0;
                for (int r = 0; r <= k-1; r++)
                {
                    sum1 = sum1 + L[k][r]*U[r][j];
                }            
                U[k][j] = (A[k][j] - sum1)/L[k][k];
            }
        }       
    }
    //print A, L, U
    printf("A = \n");
    print_matrix(n, n, A);
    printf("L = \n");
    print_matrix(n, n, L);
    printf("U = \n");
    print_matrix(n, n, U);

    return 0;
}