## Cholesky decomposition

Open discussion regarding features, bugs, issues, vendors, etc.

### Cholesky decomposition

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.
woglinde

Posts: 2
Joined: Sat Aug 12, 2006 2:04 am

### Re: Cholesky decomposition

woglinde wrote: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
buttari

Posts: 51
Joined: Tue Jul 11, 2006 2:11 pm

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.
woglinde

Posts: 2
Joined: Sat Aug 12, 2006 2:04 am

woglinde wrote: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
buttari

Posts: 51
Joined: Tue Jul 11, 2006 2:11 pm