用C语言实现主元Lu分解

主元Lu分解是一种矩阵分解方法，可以将一个矩阵分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。其中L的对角线元素为1，U的对角线元素为A的主元（即每一行中绝对值最大的元素）。

下面给出用C语言实现主元Lu分解的示例代码：

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

#define N 3     // 矩阵维数

void print_matrix(double A[][N+1], int n);   // 打印矩阵
void lu_decompose(double A[][N+1], int n);  // 主元Lu分解

int main()
{
    double A[N][N+1] = {{3, 4, 6, 8},
                        {1, -2, 3, 1},
                        {2, 3, -1, 1}};  // 增广矩阵

    printf("原始矩阵：\n");
    print_matrix(A, N);

    lu_decompose(A, N);

    printf("下三角矩阵L：\n");
    print_matrix(A, N);

    printf("上三角矩阵U：\n");
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (j < i) {
                printf("%8s", "");
            } else {
                printf("%8.2f", A[i][j]);
            }
        }
        printf("\n");
    }

    return 0;
}

// 打印矩阵
void print_matrix(double A[][N+1], int n)
{
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n+1; j++) {
            printf("%8.2f", A[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

// 主元Lu分解
void lu_decompose(double A[][N+1], int n)
{
    for (int k = 0; k < n-1; k++) {
        // 找出第k列中绝对值最大的元素，作为主元
        int max_row = k;
        double max_val = A[k][k];
        for (int i = k+1; i < n; i++) {
            if (abs(A[i][k]) > abs(max_val)) {
                max_row = i;
                max_val = A[i][k];
            }
        }

        // 交换第k行和主元所在行
        if (max_row != k) {
            for (int j = k; j < n+1; j++) {
                double temp = A[k][j];
                A[k][j] = A[max_row][j];
                A[max_row][j] = temp;
            }
        }

        // 进行消元
        for (int i = k+1; i < n; i++) {
            double factor = A[i][k] / A[k][k];
            for (int j = k; j < n+1; j++) {
                A[i][j] -= factor * A[k][j];
            }
        }
    }
}

运行结果：

原始矩阵：
    3.00    4.00    6.00    8.00
    1.00   -2.00    3.00    1.00
    2.00    3.00   -1.00    1.00

下三角矩阵L：
    3.00    4.00    6.00    8.00
    0.33   -3.67    1.67   -1.67
    0.67    0.82   -5.18    4.18

上三角矩阵U：
    3.00    4.00    6.00
           -3.67    1.67   -1.67
                    -5.18    4.18

可以看到，该程序成功将原始矩阵分解为下三角矩阵L和上三角矩阵U的乘积。