C語言矩陣求逆：LU 分解法完整實現 (Numerical Recipes 風格)

通用矩陣求逆：LU 分解 + 前向/回向替換

算法流程

1. LU 分解：將矩陣 A 分解為 P*A = L*U (L 下三角, U 上三角, P 置換矩陣)，使用部分樞軸選擇保證數值穩定。
2. 求解多個右端向量：對單位矩陣 I 的每一列 e_j，解 A * x_j = e_j，即 L*U*x_j = P*e_j。
   先解 L*y = P*e_j (前向替換)，再解 U*x_j = y (回向替換)。
3. 組裝逆矩陣：各解向量 x_j 組成 A^{-1} 的列。

核心代碼結構 (適配任意階數，含 2x2)

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

#define TINY 1.0e-20

// 動態分配矩陣 (行優先連續內存)
double **matrix(int n) {
    double **m = malloc(n * sizeof(double*));
    double *data = malloc(n * n * sizeof(double));
    for (int i = 0; i < n; ++i) m[i] = data + i * n;
    return m;
}
void free_matrix(double **m) { free(m[0]); free(m); }

// LU 分解 (帶部分樞軸) - 修改傳入矩陣為 LU 合併存儲
void ludcmp(double **a, int n, int *indx, double *d) {
    double *vv = malloc(n * sizeof(double));
    *d = 1.0;
    for (int i = 0; i < n; ++i) {
        double big = 0.0;
        for (int j = 0; j < n; ++j)
            if (fabs(a[i][j]) > big) big = fabs(a[i][j]);
        if (big == 0.0) { fprintf(stderr, "奇異矩陣\n"); exit(1); }
        vv[i] = 1.0 / big;
    }
    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < j; ++i) {
            double sum = a[i][j];
            for (int k = 0; k < i; ++k) sum -= a[i][k] * a[k][j];
            a[i][j] = sum;
        }
        double big = 0.0; int imax = j;
        for (int i = j; i < n; ++i) {
            double sum = a[i][j];
            for (int k = 0; k < j; ++k) sum -= a[i][k] * a[k][j];
            a[i][j] = sum;
            if (vv[i] * fabs(sum) >= big) { big = vv[i] * fabs(sum); imax = i; }
        }
        if (j != imax) {
            for (int k = 0; k < n; ++k) { double tmp = a[imax][k]; a[imax][k] = a[j][k]; a[j][k] = tmp; }
            *d = -(*d); vv[imax] = vv[j];
        }
        indx[j] = imax;
        if (a[j][j] == 0.0) a[j][j] = TINY;
        if (j != n-1) {
            double dum = 1.0 / a[j][j];
            for (int i = j+1; i < n; ++i) a[i][j] *= dum;
        }
    }
    free(vv);
}

// 回代求解
void lubksb(double **a, int n, int *indx, double *b) {
    int ii = -1;
    for (int i = 0; i < n; ++i) {
        int ip = indx[i];
        double sum = b[ip];
        b[ip] = b[i];
        if (ii >= 0) for (int j = ii; j < i; ++j) sum -= a[i][j] * b[j];
        else if (sum) ii = i;
        b[i] = sum;
    }
    for (int i = n-1; i >= 0; --i) {
        double sum = b[i];
        for (int j = i+1; j < n; ++j) sum -= a[i][j] * b[j];
        b[i] = sum / a[i][i];
    }
}

// 矩陣求逆入口
void inverse(double **a, int n) {
    double **y = matrix(n);
    int *indx = malloc(n * sizeof(int));
    double *col = malloc(n * sizeof(double));
    double d;

    ludcmp(a, n, indx, &d);  // a 被修改為 LU

    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < n; ++i) col[i] = 0.0;
        col[j] = 1.0;
        lubksb(a, n, indx, col);
        for (int i = 0; i < n; ++i) y[i][j] = col[i];
    }
    // 複製回 a
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            a[i][j] = y[i][j];

    free_matrix(y); free(indx); free(col);
}

// 2x2 專用簡化版 (教學/嵌入式無 malloc 場景)
void inverse_2x2(const double a[2][2], double inv[2][2]) {
    double det = a[0][0]*a[1][1] - a[0][1]*a[1][0];
    if (fabs(det) < 1e-12) { printf("奇異\n"); return; }
    inv[0][0] =  a[1][1] / det;
    inv[0][1] = -a[0][1] / det;
    inv[1][0] = -a[1][0] / det;
    inv[1][1] =  a[0][0] / det;
}

使用示例

int main() {
    // 通用版
    int n = 3;
    double **A = matrix(n);
    A[0][0]=4; A[0][1]=7; A[0][2]=2;
    A[1][0]=3; A[1][1]=5; A[1][2]=1;
    A[2][0]=2; A[2][1]=6; A[2][2]=3;
    inverse(A, n);
    printf("逆矩陣:\n");
    for (int i=0;i