problem about SPOTRF

[nl2219]

Hi all,

I was trying to use SPOTRF to compute the Cholesky factorization of a real symmetric
positive definite matrix A.

A = U**T * U, if UPLO = 'U', or
A = L * L**T, if UPLO = 'L',

where U is an upper triangular matrix and L is lower triangular.

However, the output U I got is not uppertriangular, and I commputed U**T*U, it doesn't not equal to the original matrix A.

So why it is like that please?

Below is my code.

#include <stdio.h>
#include <time.h>
#include "f2c.h"
//#include "clapack.h"
extern "C" { int spotrf_(char *uplo,integer *n, real *a,integer *lda,integer* info);}
int main(void)
{
char uplo;
integer n,lda,info;
real a[2][2]={3.,2.,2.,3.};

n=2;
uplo='U';
lda=2;
spotrf_(&uplo,&n,a[0],&lda,&info);
printf("INFO=%d\n",info);
for(int i=0;i<2;i++)
{
for(int j=0;j<2;j++)
printf("%f ",a[i][j]);
printf("\n");
}
printf("\n");
for(int i=0;i<2;i++)
{
for(int j=0;j<2;j++)
{
float sum=0.;
for(int k=0;k<2;k++)
sum+=a[k][i]*a[k][j];
printf("%f ",sum);
}
printf("\n");
}

return 0;
}

And the result I got is
U:
1.732051 2.000000
1.154701 1.290995

U**T*U
4.333333 4.954813
4.954813 5.666667

Thanks in advance.

[Julien Langou]

[1] LAPACK works with column major (not row major). This is because there is a FORTRAN story behind. All array are 1D, then you
need to go by column to describe your 2D matrix. The output of LAPACK is not

Code: Select all

U:
1.732051 2.000000
1.154701 1.290995


as you have written but

Code: Select all

U:
1.732051 1.154701
2.000000 1.290995


[2] With the 'U' flag, LAPACK does not touch the lower part of the matrix (in your case only A(2,1)). Whatever is there in input will be there in output. Your input matrix is symmetric in input so whatever you put on the upper part is assumed to be the same in the lower part, so no need to look at it. The output matrix is triangular so you know that below the diagonal the entries are zeros. No need to put explicit zeros. Using this LAPACK xPOTRF performs less memory reference, moreover the user is free to store whatever he wants in the lower part of the matrix. So although LAPACK returns

Code: Select all

U:
1.732051 1.154701
2.000000 1.290995


it means

Code: Select all

U:
1.732051 1.154701
       0 1.290995


You can check that

Code: Select all

U:
( 1.732051        0 ) ( 1.732051 1.154701 )    (  3  2 )
( 1.154701 1.290995 ) (        0 1.290995 ) =  (  2  3 )