通过 Lapack dpotrf 在 C++ 中进行 Choleskey 分解给出无效结果

Question

我正在尝试使用 Lapack 的 dpotrf 函数对一维双数组进行 choleskey 分解。它似乎适用于某些情况，但在其他情况下它有奇怪的行为。

这是我的代码：

int main(){
    int N=2;
    int INFO;
    double A2 [2*2] = {
            1,1,
            1,4
        };
        printf("%f %f\n", A2[0], A2[1]);
        printf("%f %f\n", A2[2], A2[3]);
        printf("\n");

    char UPLO='L';
    LAPACK_dpotrf(&UPLO,&N,A2,&N,&INFO);

    printf("%f %f\n", A2[0], A2[1]);
    printf("%f %f\n", A2[2], A2[3]);
    printf("\n");

    printf("2x2 working?: ");
    printf("%i\n",INFO);

// double A3 [3*3] = {
    //     1,1,1,
    //     1,1,1,
    //     1,1,4
    // };

  double A3 [3*3] = {
        2,-1, 0,
        -1, 2, -1, 
        0, -1, 2
        };

    printf("%f %f %f\n", A3[0], A3[1], A3[2]);
    printf("%f %f %f\n", A3[4], A3[5], A3[6]);
    printf("%f %f %f\n", A3[7], A3[8], A3[9]);
    printf("\n");

    int N3=3;
    LAPACK_dpotrf(&UPLO,&N3,A3,&N3,&INFO);

    printf("%f %f %f\n", A3[0], A3[1], A3[2]);
    printf("%f %f %f\n", A3[4], A3[5], A3[6]);
    printf("%f %f %f\n", A3[7], A3[8], A3[9]);
    printf("\n");

    printf("%i\n",INFO);
    

    return 0;
}

并通过以下方式构建：

g++ main.cpp -lblas -llapack -llapacke

2x2 的情况给出：

1.000000 1.000000
1.000000 1.732051

从技术上讲，这是不正确的，因为它不是下三角，但下部是正确的

它正确地（也许是偶然的）识别出它不能对注释掉的 A3 矩阵进行操作。

它正确地表明它可以在有效的A3上进行操作，但结果值为：

1.414214 -0.707107 0.000000
1.224745 -0.816497 0.000000
-1.000000 1.154701 0.000000

有效的 Cholesky 分解是（根据 Matlab）：

1.4142         0         0
-0.7071    1.2247        0
      0   -0.8165    1.1547

我怀疑函数没有损坏，但我没有正确的参数。
我怀疑它可能是 LDA 变量，因为我不知道这是什么，而其他 LAPACK 函数通常只是将其设置为与 N 相同的值。

有谁知道可能导致此问题的原因吗？

Answer 1

this is incorrect, as it isn't lower triangular

正确，因为上三角的原始数据没有被覆盖。

在 3x3 的情况下，您打印了错误的元素。仔细看看。

通过 Lapack dpotrf 在 C++ 中进行 Choleskey 分解给出无效结果

Choleskey Decomposition in C++ via Lapack dpotrf gives invalid result

c++

matrix

lapack

lapacke

matrix-factorization