C语言三角分解非奇异矩阵为L、U矩阵

#include<stdio.h>
# define N 4
int main (void)

{

double A[N][N]  = {0};
double L[N][N] = {0};
double U[N][N] = {0};
double M[N][N] = {0};
double P[N][N] = {0};
double sum1 = 0,sum2 = 0;
int i = 0,j = 0;
int r,k;

printf("请输入矩阵A：\n");
for(i=0;i<N;i++)
{
    printf("第%d行：",i+1);
    for(j=0;j<N;j++)
    {
        scanf("%lf",&A[i][j]);
    }

}
printf("你输入的矩阵是：\n");
for(i = 0;i<N;i++)
{
    for(j = 0;j<N;j++)
    {
        printf("%lf\t",A[i][j]);
    }
    printf("\n");
}

/*计算U和L*/
/*U的第一行*/
for(j = 0;j<N;j++)
{
    U[0][j] = A[0][j];
}

/*L的第一列*/
for(i = 0;i<N;i++)
{
    L[i][0] = A[i][0]/U[0][0];
}

/*U的第r+1行、L的第r+1列*/
for(r=1;r<N;r++)
{
    for(j=0;j<N;j++)
    {
        for(k = 0;k<r;k++)
        {
            sum1 += L[r][k]*U[k][j];
        }
        M[r][j] = sum1;
        U[r][j] = A[r][j] - M[r][j];    
        sum1 = 0;
    }
    if(j == 4)
    {
        j = r;
        for(r = 0;r<N;r++)
        {
            for(k = 0;k<j;k++)
            {
                sum2 += L[r][k]*U[k][j];
            }
            P[r][j] = sum2;
            L[r][j] = (A[r][j] - P[r][j])/U[j][j]; 
            sum2 = 0;
        }
    }
    r = j;
}

printf("矩阵U是：\n");
for(i = 0;i<N;i++)
{
    for(j = 0;j<N;j++)
    {
        printf("%.3lf\t",U[i][j]);
    }
    printf("\n");
}

printf("矩阵L是：\n");
for(i = 0;i<N;i++)
    {
    for(j = 0;j<N;j++)
    {
        printf("%.3lf\t",L[i][j]);
    }
    printf("\n");
}


return 0;
}
/*咱就说：就这？！*/
————————————————
版权声明：本文为CSDN博主「amen_yanwang」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/amen_yanwang/article/details/123731959