使用LAPACK sgetrf进行LU分解时为何出现重复枢轴值？

我正在使用LAPACK库调用sgetrf routine计算LU分解（预期公式：A = L*U*P），C代码得到的LU矩阵与MATLAB结果一致，但枢轴数组出现重复值：[3, 2, 2, 3]。该枢轴数组用于指示LU矩阵的行读取顺序以求解线性方程组Ax=b，但重复值导致C语言求解结果与MATLAB的linsolve结果不一致。

相关代码及输出 #

LU分解C代码 #

bool lup(float A[], float LU[], size_t P[], size_t row) {
    integer m = row, lda = row, n = row;
    size_t rowrow = row * row;
    integer* ipiv = (integer*)malloc(row * sizeof(integer));
    integer info;
    memcpy(LU, A, row * row * sizeof(float));
    tran(LU, row, row); /* Transpose - Make it column major */
    sgetrf_(&m, &n, LU, &lda, ipiv, &info);
    tran(LU, row, row); /* Transpose - Make it row major */
    size_t i;
    printf("P:\n");
    for (i = 0; i < row; i++) {
        P[i] = ipiv[i] - 1;
        printf("%i ", P[i]);
    }
    printf("\n");
    printf("LU:\n");
    print(LU, row, row);

    free(ipiv);
    return info == 0;
}

输出结果 #

P:
3 2 2 3

LU:
0.6280800       0.7259100       0.2024400       0.9674300
0.6982390       0.4815813       0.3990585       -0.4290273
0.5180869       0.2486716       0.3052040       -0.1796059
0.7556680       0.4116502       -0.0234923      0.0720639

线性方程组求解C代码 #

/*
 * This solves Ax=b with LUP-decomposition
 * A [m*n]
 * x [n]
 * b [m]
 * n == m
 * Returns true == Success
 * Returns false == Fail
 */
bool linsolve_lup(float A[], float x[], float b[], size_t row) {
    int32_t i, j;

    float* LU = (float*)malloc(row * row * sizeof(float));
    size_t* P = (size_t*)malloc(row * sizeof(size_t));
    bool ok = lup(A, LU, P, row);

    /* forward substitution with pivoting */
    for (i = 0; i < row; ++i) {
        x[i] = b[P[i]];
        for (j = 0; j < i; ++j) {
            x[i] = x[i] - LU[row * P[i] + j] * x[j];
        }
    }

    /* backward substitution with pivoting */
    for (i = row - 1; i >= 0; --i) {
        for (j = i + 1; j < row; ++j) {
            x[i] = x[i] - LU[row * P[i] + j] * x[j];
        }
        x[i] = x[i] / LU[row * P[i] + i];
    }

    free(LU);
    free(P);

    return ok;
}

MATLAB对比代码及输出 #

A = [0.47462,   0.74679,   0.31008,   0.63073,
        0.32540,   0.49584,   0.50932,   0.21492,
        0.43855,   0.98844,   0.54041,   0.24647,
        0.62808,   0.72591,   0.20244,   0.96743];

    b = [1.588964,
         0.901248,
         0.062029,
         0.142180];

    x = linsolve(A, b)

-44.1551
    6.1363
   15.1259
   21.0440

问题分析与解决方案 #

问题根源 #

1. 转置操作导致枢轴含义错位
   你对矩阵执行了两次转置，相当于在转置后的矩阵上调用sgetrf。sgetrf输出的ipiv数组记录的是转置矩阵的行交换操作，对应原矩阵的列交换操作，但你却将其当作原矩阵的行交换索引使用，这完全违背了sgetrf的输出定义，因此出现重复的枢轴值。

2. 分解公式的误解
   LAPACK的sgetrf实现的是行置换分解：PA=LU（置换矩阵P左乘A得到LU），而你代码中假设的是A=L*U*P（列置换矩阵P右乘），两者的置换矩阵逻辑完全相反，这是求解结果不一致的核心原因。

修正方案 #

方案1：修正LU分解函数，正确处理存储格式

直接将行优先的输入矩阵转换为列优先格式调用sgetrf，并正确转换ipiv的含义（对应原矩阵的行交换）：

bool lup(float A[], float LU[], size_t P[], size_t row) {
    integer m = row, lda = row, n = row;
    integer* ipiv = (integer*)malloc(row * sizeof(integer));
    integer info;

    // 行优先转列优先
    for (size_t i = 0; i < row; i++) {
        for (size_t j = 0; j < row; j++) {
            LU[j * row + i] = A[i * row + j];
        }
    }

    sgetrf_(&m, &n, LU, &lda, ipiv, &info);

    // 列优先转回行优先
    float* temp = (float*)malloc(row * row * sizeof(float));
    memcpy(temp, LU, row * row * sizeof(float));
    for (size_t i = 0; i < row; i++) {
        for (size_t j = 0; j < row; j++) {
            LU[i * row + j] = temp[j * row + i];
        }
    }
    free(temp);

    // 1-based转0-based行置换索引
    printf("P:\n");
    for (size_t i = 0; i < row; i++) {
        P[i] = ipiv[i] - 1;
        printf("%zu ", P[i]);
    }
    printf("\n");
    printf("LU:\n");
    print(LU, row, row);

    free(ipiv);
    return info == 0;
}

方案2：修正求解逻辑，匹配PA=LU分解

原求解代码错误地使用P[i]索引LU矩阵的行，正确的逻辑应基于PA=LU的分解步骤：先对b应用行置换，再执行前向/后向替换：

bool linsolve_lup(float A[], float x[], float b[], size_t row) {
    int32_t i, j;

    float* LU = (float*)malloc(row * row * sizeof(float));
    size_t* P = (size_t*)malloc(row * sizeof(size_t));
    bool ok = lup(A, LU, P, row);

    // 对b应用行置换P
    float* b_perm = (float*)malloc(row * sizeof(float));
    for (i = 0; i < row; ++i) {
        b_perm[i] = b[P[i]];
    }

    // 前向替换求解Ly = b_perm（L是单位下三角，存储在LU下三角）
    float* y = (float*)malloc(row * sizeof(float));
    for (i = 0; i < row; ++i) {
        y[i] = b_perm[i];
        for (j = 0; j < i; ++j) {
            y[i] -= LU[i * row + j] * y[j];
        }
    }

    // 后向替换求解Ux = y（U是上三角，存储在LU上三角）
    for (i = row - 1; i >= 0; --i) {
        x[i] = y[i];
        for (j = i + 1; j < row; ++j) {
            x[i] -= LU[i * row + j] * x[j];
        }
        x[i] /= LU[i * row + i];
    }

    free(LU);
    free(P);
    free(b_perm);
    free(y);

    return ok;
}

验证 #

修正后，枢轴数组将不再出现重复值，求解结果会与MATLAB的linsolve(A,b)完全一致。

内容的提问来源于stack exchange，提问作者euraad