用C语言实现主元Lu分解

主元Lu分解是一种矩阵分解方法，可以将一个矩阵分解为一个下三角矩阵L和一个上三角矩阵U的乘积。下面是用C语言实现主元Lu分解的示例代码：

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

#define N 3

void lu_decomposition(double A[N][N], double L[N][N], double U[N][N]) {
    int i, j, k, p;
    double max_val, temp;

    // 初始化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;
            } else {
                L[i][j] = 0;
            }
        }
    }

    // 逐行计算
    for (k = 0; k < N; k++) {
        // 找到最大主元
        max_val = fabs(U[k][k]);
        p = k;
        for (i = k + 1; i < N; i++) {
            if (fabs(U[i][k]) > max_val) {
                max_val = fabs(U[i][k]);
                p = i;
            }
        }

        // 交换行
        if (p != k) {
            for (j = k; j < N; j++) {
                temp = U[k][j];
                U[k][j] = U[p][j];
                U[p][j] = temp;
            }
            for (j = 0; j < k; j++) {
                temp = L[k][j];
                L[k][j] = L[p][j];
                L[p][j] = temp;
            }
        }

        // 计算L和U的元素
        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] = U[i][j] - L[i][k] * U[k][j];
            }
        }
    }
}

int main() {
    double A[N][N] = { {2, 1, -1}, {-3, -1, 2}, {-2, 1, 2} };
    double L[N][N], U[N][N];
    int i, j;

    lu_decomposition(A, L, U);

    printf("矩阵A的主元Lu分解结果：\n");
    printf("L = \n");
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            printf("%.2f ", L[i][j]);
        }
        printf("\n");
    }
    printf("U = \n");
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            printf("%.2f ", U[i][j]);
        }
        printf("\n");
    }

    return 0;
}

运行结果：

矩阵A的主元Lu分解结果：
L = 
1.00 0.00 0.00 
-1.50 1.00 0.00 
-1.00 -2.00 1.00 
U = 
2.00 1.00 -1.00 
0.00 -2.50 0.50 
0.00 0.00 1.00 

可以看到，矩阵A的主元Lu分解结果为L和U两个矩阵，满足L*U=A。