LAPACKE_sgetrf 是如何工作的？

Question

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

我写了下面的程序，

#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=

和U=

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

我愿意

$ gcc -Wall  sgetrf.c -llapacke -llapack

我看到一个= 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 实际上 return 矩阵的 LU 因式分解具有部分主元。它 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

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 是一个置换矩阵，returned in ipiv。事实上，这样的分解总是可能的。

LAPACKE_sgetrf 是如何工作的？

How does LAPACKE_sgetrf work?

lapack