Solving a linear system with dpotrs (Cholesky factorization)

Question

I am trying to solve a linear system of equations with clapack.

My code is as follows:

//ATTENTION: matrix in column-major
double A[3*3]={  2.0, -1.0,  0.0,
                0.0,  2.0, -1.0,
                0.0,  0.0,  2.0},
b[3]={1.0,2.0,3.0};

integer n=3,info,nrhs=1; char uplo='L';

dpotrf_("L", &n, A, &n, &info);
dpotrs_("L", &n, &nrhs, A, &n, b, &n, &info);

printf("Solution: %10.4f %10.4f %10.4f",b[0], b[1], b[2]);

As asked in this question, it is necessary to first factorize the matrix. dpotrf is supposed to factorize, dpotrs solves the system using the factorized matrix.

However, my result

2.5   4.0   3.5

is clearly wrong, I checked it here with WolframAlpha

Where is my mistake?

Answer 1

Curiously, 2.5 4.0 3.5 is the solution of rhs 1 2 3...if the matrix is

 2   -1  0
 -1  2   -1
 0   -1  2

dpotrf and dpotrs are made for symetric positive definite matrices...

I would suggest to use dgetrf and dgetrs instead. In your case :

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

int main()
{
    //ATTENTION: matrix in column-major
    double A[3*3]={  2.0, -1.0,  0.0,
            0,  2.0, -1.0,
            0.0,  0,  2.0},
            b[3]={1.0,2,3.0};

    int n=3,info,nrhs=1; char uplo='L';
    int ipvs[3];
    dgetrf_(&n, &n, A, &n,ipvs, &info);
    printf("info %d\n",info);
    dgetrs_("T", &n, &nrhs, A, &n,ipvs, b, &n, &info);
    printf("info %d\n",info);
    printf("Solution: %10.4f %10.4f %10.4f\n",b[0], b[1], b[2]);

    return 0;
}

The solution i got is 1.3750 1.75 1.5. If it is not the column major order, switch to "N" instead of "T". Then, the solution becomes 0.5 1.25 2.125

Choose your favorite !

Bye,

Francis