LU分解求矩阵逆 C语言实现

LU分解是一种求解矩阵逆的方法之一，下面是一个使用C语言实现的LU分解求矩阵逆的示例代码：

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

#define N 3

void LUdecomposition(double A[N][N], double L[N][N], double U[N][N]) {
    int i, j, k;
    double sum;

    // 初始化U矩阵为A矩阵
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            U[i][j] = A[i][j];
        }
    }

    // 初始化L矩阵为单位矩阵
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            if (i == j) {
                L[i][j] = 1.0;
            } else {
                L[i][j] = 0.0;
            }
        }
    }

    // LU分解
    for (k = 0; k < N - 1; k++) {
        for (i = k + 1; i < N; i++) {
            L[i][k] = U[i][k] / U[k][k];
            for (j = k; j < N; j++) {
                U[i][j] -= L[i][k] * U[k][j];
            }
        }
    }
}

void inverse(double A[N][N], double A_inv[N][N]) {
    int i, j, k;
    double L[N][N], U[N][N], y[N][N], sum;

    // LU分解
    LUdecomposition(A, L, U);

    // 求解Ly=I
    for (k = 0; k < N; k++) {
        for (i = 0; i < N; i++) {
            sum = 0.0;
            for (j = 0; j < i; j++) {
                sum += L[i][j] * y[j][k];
            }
            if (i == k) {
                y[i][k] = 1.0 - sum;
            } else {
                y[i][k] = -sum;
            }
        }
    }

    // 求解Ux=y
    for (k = 0; k < N; k++) {
        for (i = N - 1; i >= 0; i--) {
            sum = 0.0;
            for (j = i + 1; j < N; j++) {
                sum += U[i][j] * A_inv[j][k];
            }
            A_inv[i][k] = (y[i][k] - sum) / U[i][i];
        }
    }
}

int main() {
    double A[N][N] = {{2.0, -1.0, 0.0},
                      {-1.0, 2.0, -1.0},
                      {0.0, -1.0, 2.0}};
    double A_inv[N][N];
    int i, j;

    inverse(A, A_inv);

    printf("A_inv:\n");
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            printf("%.2f\t", A_inv[i][j]);
        }
        printf("\n");
    }

    return 0;
}

该代码中，`LUdecomposition()`函数实现了LU分解，`inverse()`函数实现了矩阵逆的求解，`main()`函数中给出了一个示例矩阵A，调用`inverse()`函数求解其逆矩阵并输出。