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

Question

I'm trying to take the choleskey decomposition of a 1D double array using Lapack's dpotrf function. It seems to work for some cases, but others it has strange behavior.

Here is my code:

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;
}

and is built via:

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

The 2x2 case gives:

1.000000 1.000000
1.000000 1.732051

Technically, this is incorrect, as it isn't lower triangular, but the lower part is correct

It correctly (perhaps by accident), identifies that it cannot operate on the commented out A3 matrix.

It correctly indicates that it can operate on the valid A3, but the resulting value is:

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

The valid Cholesky decomposition is (according to Matlab):

1.4142         0         0
-0.7071    1.2247        0
      0   -0.8165    1.1547

I suspect the function is not broken, but I don't have the right arguments.
I suspect it might be the LDA variable, as I have no idea what this is, and other LAPACK functions typically just have this set to the same value of N.

Does anyone know what might be causing this issue?

Answer 1

this is incorrect, as it isn't lower triangular

It is correct since the original data in the upper triangular is not overwritten.

In the 3x3 case you are printing the wrong elements. Take a close look.