Cholesky decomposition

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

Cholesky decomposition

Postby woglinde » Sat Aug 12, 2006 2:34 am

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

Postby buttari » Sat Aug 12, 2006 10:04 am

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

Postby woglinde » Sat Aug 12, 2006 2:29 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

Postby buttari » Sat Aug 12, 2006 2:48 pm

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


Return to User Discussion

Who is online

Users browsing this forum: No registered users and 4 guests