CODEWITHSUNDEEP
lapack
ericj
ericj

Reputation: 2301

How does LAPACKE_sgetrf work?

I just started with lapack. I thought let's start with LU decomposition, of the following matrix A:

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

I wrote the follwinng program,

#include <stdio.h>
#include <lapacke/lapacke.h>

int main (int argc, const char * argv[]) {
  float a[4][4]={
{1,1,0,3},
{2,1,-1,1},
{3,-1,-1,2},
{-1,2,3,-1}
};
  lapack_int m=4,n=4,lda=4,info;
  int i,j;
  lapack_int*ipiv=(lapack_int*)malloc(4*4*sizeof(lapack_int));

   m = 4;
   n = 4;
   lda = 4;

   info= LAPACKE_sgetrf(LAPACK_ROW_MAJOR,m,n,*a,lda,ipiv);

   for(i=0;i<m;i++)
   {
      for(j=0;j<n;j++)
      {
         printf("%lf ",a[i][j]);
      }
      printf("\n");
   }
   return(info);
}

This should give

L=

1 0 0 0
2 1 0 0
3 4 1 0
-1 3 0 1

and U=

1 1 0 3
0 -1 -1 -5
0 0 3 13
0 0 0 -13

But in the docs I read that L en U are returned in A. How? Maybe skip the 1's on the diagonal of L en merge the two. Is that true?

I do

$ gcc -Wall  sgetrf.c -llapacke -llapack

I see a= 3.000000 -1.000000 -1.000000 2.000000 0.666667 1.666667 -0.333333 -0.333333 -0.333333 1.000000 3.000000 0.000000 0.333333 0.800000 0.200000 2.600000

I does not make sense to me.

Upvotes: 2

Views: 564

Answers (1)

francis
francis

Reputation: 9817

Your program is perfectly correct and LAPACKE_sgetrf actually return an LU factorization of the matrix with partial pivoting. It returned :

3     -1   -1     2
2/3   5/3 -1/3  -1/3
-1/3   1    3     0
1/3   0.8  0.2   2.6

The L matrix is:

1       0    0    0
2/3    1   -1/3  -1/3
-1/3   1    1     0
1/3   0.8  0.2    1

The U matrix is:

3     -1   -1     2
0     5/3 -1/3  -1/3
0     0     3     0
0     0     0    2.6

The product is similar to A except for the permutation of the rows:

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

The LU factorization you provided is correct (except for the last line of L, which must be -1 -3 0 1 instead of -1 3 0 1). But such a LU factorization is not always possible (See wikipedia). This is the reason why LAPACK compute a PLU factorization, where P is a permutation matrix, returned in ipiv. Indeed, such a factorization is always possible.

Upvotes: 1

Related Questions