C语言实现LU三角分解法（计算方法）

#include <stdio.h>
#include <stdlib.h>
#include<stdbool.h>
#define MAX 10

double A[MAX][MAX];
double b[MAX];
double X[MAX];
double Y[MAX];
double U[MAX][MAX];
double L[MAX][MAX];
int NUM;

void Input_Matrix()//输入矩阵
{
    int i,j;
    printf("系数矩阵A的阶数:\n");
    scanf("%d",&NUM);
    for(i=1; i<=NUM; i++)
    {
        printf("系数矩阵A的第%d行元素:\n",i);
        for(j=1; j<=NUM; j++)
            scanf("%lf",&A[i-1][j-1]);
    }
    printf("右端项b:\n");
    for(i=1; i<=NUM; i++)
    {
        scanf("%lf",&b[i-1]);
    }
    printf("\n");
    printf("输入的系数矩阵A:\n");
    for(i=0; i<NUM; i++)
    {
        for(j=0; j<NUM; j++)
            printf("%.4lf\t",A[i][j]);
        printf("\n");
    }
    printf("\n");
    printf("输入的右端项b:\n");
    for(i=0; i<NUM; i++)
        printf("%.4lf\n",b[i]);
}

bool Direct()
{
    bool flag = true;
    int r,i,k;
    double sum_u, sum_l;
    for(r=0; r<NUM; r++)
        U[0][r]=A[0][r];
    for(r=1; r<NUM; r++)
        L[r][0] = A[r][0]/A[0][0];
    for(r=1; r<NUM; r++)
    {
        for(i=r; i<NUM; i++)
        {
            sum_u = 0;
            for(k=0; k<r; k++)
                sum_u += L[r][k]*U[k][i];
            U[r][i] = A[r][i]-sum_u;
        }
        for(i=r+1; i<NUM&&r!=NUM-1; i++)
        {
            sum_l = 0;
            for(k=0; k<r; k++)
                sum_l += L[i][k]*U[k][r];
            L[i][r] = (A[i][r]-sum_l)/U[r][r];
        }
    }
    for(i=0; i<NUM; i++)
        L[i][i]=1.0;
    for(r=0; r<NUM; r++)
    {
        if(U[r][r]==0) //等价于判断矩阵的顺序主子式是否为0
        {
            flag = false;
            return flag;
        }
    }
    //回代求解
    Y[0] = b[0];
    double sum;
    for(i=1; i<NUM; i++)
    {
        sum = 0;
        for(k=0; k<i; k++)
            sum += L[i][k]*Y[k];
        Y[i] = b[i]-sum;
    }
    X[NUM-1] = Y[NUM-1]/U[NUM-1][NUM-1];
    for(i=NUM-2; i>=0; i--)
    {
        sum = 0;
        for(k=i+1; k<NUM; k++)
            sum += U[i][k]*X[k];
        X[i] = (Y[i]-sum)/U[i][i];
    }
    return flag;
}

int main()
{
    int i,r;
    Input_Matrix();
    if (Direct())
    {
        printf("\n");
        printf("Result:\n");
        for (i=0; i<NUM; i++)
            printf("%.4lf\n",X[i]);
        printf("\n");
        printf("L:\n");
        for(i=0; i<NUM; i++)
        {
            for(r=0; r<NUM; r++)
                printf("%.4lf\t",L[i][r]);
            printf("\n");
        }
        printf("\n");
        printf("U:\n");
        for(i=0; i<NUM; i++)
        {
            for(r=0; r<NUM; r++)
                printf("%.4lf\t",U[i][r]);
            printf("\n");
        }
    }
    else
        printf("Failed.n");
    return 0;
}