LAPACKE_sgetrf是如何工作的？

Question

我刚开始使用lapack。我想让我们从LU分解开始，分解如下矩阵A：

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

我写了下面的程序，

#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);
}

这应该会给出

L=

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

和U=

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

但在文档中我读到L en U在A中返回。如何返回？也许可以跳过L的对角线上的1，然后合并这两个。这是真的吗？

我知道

$ gcc -Wall  sgetrf.c -llapacke -llapack

我看到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

我对我来说没有任何意义。

Answer 1

您的程序是完全正确的，LAPACKE_sgetrf实际上返回了带有部分旋转的矩阵的LU分解。它返回：

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

L矩阵是：

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

U矩阵是：

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

除了行的排列之外，乘积类似于A：

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

您提供的LU分解是正确的(除了L的最后一行，它必须是-1 -3 0 1而不是-1 3 0 1)。但这样的LU分解并不总是可能的(参见wikipedia)。这就是为什么LAPACK计算一个PLU分解的原因，其中P是一个置换矩阵，在ipiv中返回。事实上，这样的因式分解总是可能的。