Cholesky decomposition

[woglinde]

Hello,
I had a question about the implementation of Cholesky decomposition in C using LAPACK. The LAPACK function for Cholesky decomposition that I'm using is spotrf.

The value of INFO is set to zero when the input matrix is symmetric positive definite and there are no errors. Suppose that the matrix for which I want to find the Cholesky decomposition is:

cov = [1 3 4;
3 13 16;
4 16 21]

This is symmetric and positive definite though when I pass it to spotrf, the value of INFO gets set to 1, signifying that the matrix cov is *not* positive definite. The situation is the same even the matrix that I pass (i.e., cov) is the identity matrix.

The relevant code that I'm using is:

int INFO, size;
char UPLO;
UPLO = 'L';
size = 3;
spotrf_(&UPLO, &size, cov, &size, &INFO);


I'm not sure why INFO is being set to 1 on return. Any help would be greatly appreciated.

Thanks,
Rahul.

[buttari]

Hello,
I had a question about the implementation of Cholesky decomposition in C using LAPACK. The LAPACK function for Cholesky decomposition that I'm using is spotrf.

The value of INFO is set to zero when the input matrix is symmetric positive definite and there are no errors. Suppose that the matrix for which I want to find the Cholesky decomposition is:

cov = [1 3 4;
3 13 16;
4 16 21]

This is symmetric and positive definite though when I pass it to spotrf, the value of INFO gets set to 1, signifying that the matrix cov is *not* positive definite. The situation is the same even the matrix that I pass (i.e., cov) is the identity matrix.

The relevant code that I'm using is:

int INFO, size;
char UPLO;
UPLO = 'L';
size = 3;
spotrf_(&UPLO, &size, cov, &size, &INFO);


I'm not sure why INFO is being set to 1 on return. Any help would be greatly appreciated.

Thanks,
Rahul.

Rahul,
can you please attach the piece of code you wrote? how is cov declared/filled?

alfredo

[woglinde]

Hi Alfredo,
Here's the code:

void main()
{
int INFO,size;
char UPLO;
double *cov;

size = 3;
cov = (double *)calloc(size * size, sizeof(double));

cov[0] = 1;
cov[1] = 3;
cov[2] = 4;
cov[3] = 3;
cov[4] = 13;
cov[5] = 16;
cov[6] = 4;
cov[7] = 16;
cov[8] = 21;

UPLO = 'L';
spotrf_(&UPLO,&size,cov,&size,&INFO);
}


Thanks,
Rahul.

[buttari]

Hi Alfredo,
Here's the code:

void main()
{
int INFO,size;
char UPLO;
double *cov;

size = 3;
cov = (double *)calloc(size * size, sizeof(double));

cov[0] = 1;
cov[1] = 3;
cov[2] = 4;
cov[3] = 3;
cov[4] = 13;
cov[5] = 16;
cov[6] = 4;
cov[7] = 16;
cov[8] = 21;

UPLO = 'L';
spotrf_(&UPLO,&size,cov,&size,&INFO);
}


Thanks,
Rahul.

So,
the answer is easy. if cov is a double precision matrix, then you have to use the double precision routine that is DPOTRF and not SPOTRF. The LAPACK routines interface follows the convention that whatever begins by D is for double precision data and whatever begins by S is for single precision. Read the LAPACK user's guide for further info about the interfaces.

Regards,
alfredo