为什么我会得到两个相同类型的主元值进行 LU 分解？ LAPACK sgetrf

Question

我正在使用 LAPACK，我得到了两种相同的主值。我正在使用 LAPACK 例程

sgetrf

来计算 LU 分解

A = L*U*P

下面的

C

代码给出的结果与我的 MATLAB 代码相同。除了一件事之外，主元数组有两个相同的值。

>> A
A =

   0.4746   0.7468   0.3101   0.6307
   0.3254   0.4958   0.5093   0.2149
   0.4385   0.9884   0.5404   0.2465
   0.6281   0.7259   0.2024   0.9674

>> LU = lu(A)
LU =

   0.628080   0.725910   0.202440   0.967430
   0.698239   0.481581   0.399058  -0.429027
   0.518087   0.248672   0.305204  -0.179606
   0.755668   0.411650  -0.023492   0.072064

>>


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

问题：

为什么我会得到两个相同的值作为主元？这没有道理。

P:
3 2 2 3

枢轴应该索引我将如何读取

LU

矩阵的行。 有了这个，我可以求解线性方程组

Ax = b

使用此

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) {
    /* Decleration */
    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 */
    free(LU);
    free(P);

    return ok;
}

但问题是，仅仅因为枢轴

P

，我不会得到与此 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)

Output: 

  -44.1551
    6.1363
   15.1259
   21.0440

Answer 1

“一个常见的错误是将 piv 视为排列向量 - 事实并非如此。相反，它是分解过程中发生的交换的顺序记录。可接受的 piv 向量可以具有相同的所有值，即方程的数量。因此，第一行与最后一行交换。接下来，第二行与最后一行交换。然后是第三行，依此类推”。

本文复制自英特尔聊天论坛上有关 MKL 和 getrf 的现已失效的网络链接。 （即我没有发布损坏的链接）

下面的伪代码应该有助于阐明如何将行更改的“序列”转换为（唯一）索引向量（您正在寻找的东西）。

    CVector<int, N> ExtractAsIndices ( )
    {
        // Here, we are generating from the "sequential record", the indices that mirror the 1's in the permutation matrix

        CVector<int, N> piv;
        LinSpace ( piv, 0, N-1 ) ;
        for ( int n = 0; n < N; ++n )
        {
            const int i = GetPivot ( n ) ;   // GetPivot accesses your sequential record (ie 3 2 2 3 in your example)
            _ASSERTE ( i >= n ) ;   // We never Pivot with a Row prior to the current one
            if ( i != n )
                swap ( piv ( n ), piv ( i ) ) ;
        }
        return piv;
    }

此示例返回一个置换矩阵

    CMatrix<T, N, N> ExtractAsPermutation ( )
    {
        CMatrix<T, N, N> piv ( 0 ) ;
        const CVector<int, N> indices = ExtractAsIndices ( ) ;
        for ( int n = 0; n < N; ++n )
            piv ( n, indices ( n ) ) = T ( 1 ) ;
        return piv;
    }

注意：这些示例使用基于 0 的索引（来自 C++），但 Matlab 使用基于 1 的索引，因此您需要考虑到这一点。