c++代码实现lu分解法解线性方程组

以下是C代码实现LU分解法解线性方程组：

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

#define N 3 // 方程组的阶数

void lu(double a[][N], double l[][N], double u[][N]) {
    int i, j, k;
    double sum;

    // 计算U矩阵的第一行
    for (j = 0; j < N; j++) {
        u[0][j] = a[0][j];
    }

    // 计算L矩阵的第一列
    for (i = 1; i < N; i++) {
        l[i][0] = a[i][0] / u[0][0];
    }

    // 计算L矩阵和U矩阵的其它元素
    for (i = 1; i < N; i++) {
        for (j = 1; j < N; j++) {
            if (i <= j) {
                sum = 0;
                for (k = 0; k < i; k++) {
                    sum += l[i][k] * u[k][j];
                }
                u[i][j] = a[i][j] - sum;
                if (i == j) {
                    l[i][j] = 1;
                } else {
                    l[i][j] = 0;
                }
            } else {
                sum = 0;
                for (k = 0; k < j; k++) {
                    sum += l[i][k] * u[k][j];
                }
                l[i][j] = (a[i][j] - sum) / u[j][j];
                u[i][j] = 0;
            }
        }
    }
}

void solve(double l[][N], double u[][N], double b[], double x[]) {
    int i, j;
    double y[N], sum;

    // 解Ly=b
    y[0] = b[0];
    for (i = 1; i < N; i++) {
        sum = 0;
        for (j = 0; j < i; j++) {
            sum += l[i][j] * y[j];
        }
        y[i] = b[i] - sum;
    }

    // 解Ux=y
    x[N-1] = y[N-1] / u[N-1][N-1];
    for (i = N-2; i >= 0; i--) {
        sum = 0;
        for (j = i+1; j < N; j++) {
            sum += u[i][j] * x[j];
        }
        x[i] = (y[i] - sum) / u[i][i];
    }
}

int main() {
    double a[N][N] = {{1, 2, 3}, {4, 5, 6}, {7, 8, 10}}; // 系数矩阵
    double b[N] = {6, 15, 29}; // 常数矩阵
    double l[N][N], u[N][N]; // LU分解后的矩阵
    double x[N]; // 解向量
    int i, j;

    // LU分解
    lu(a, l, u);

    // 解线性方程组
    solve(l, u, b, x);

    // 输出解向量
    for (i = 0; i < N; i++) {
        printf("x%d = %f\n", i+1, x[i]);
    }

    return 0;
}

注：代码中的方程组为：

1x1 + 2x2 + 3x3 = 6
4x1 + 5x2 + 6x3 = 15
7x1 + 8x2 + 10x3 = 29

输出结果为：

x1 = -2.000000
x2 = 3.000000
x3 = 1.000000