MKL中getrf省略参数ipiv的问题

FLY

如果函数getrf参数中没有ipiv作为返回值放入getri中的ipiv里，计算出的结果是错误的，即使getrf中的ipiv是可选的。这是为什么？

FLY

vvt 发表于 2014-7-23 10:55
1. F77 语法本身，对接口的检查就不严格。建议使用 F95 接口
2. 是的，只有第一个参数是必选的。可选参数的 ...

      我输出了getrf返回的ipiv的值，发现第一个数据由3变成了2。而在getri的fortran77接口中明确说明ipiv是getrf的返回值，而在fortran95接口中，则没有强调。不知道这有没有什么联系。

fcode

[Fortran] 纯文本查看 复制代码

call getrf(a,ipiv,info)!矩阵分解成LU
!call getrf(a)!矩阵分解成LU
call getri(a,ipiv,info)!求逆阵

我这里两行代码计算结果是一样的。
  0.1764706     -5.8823526E-02  0.2941177
  0.2941176      0.2352941     -0.1764706
-0.2352941      0.4117647     -5.8823526E-02
你再检查一下看是不是别的问题吧

FLY

fcode 发表于 2014-7-23 15:11

[Fortran] 纯文本查看 复制代码

call getrf(a,ipiv,info)!矩阵分解成LU
!call getrf(a)!矩阵分解成LU
call getr ...

我的源代码如下：

[Fortran] 纯文本查看 复制代码

Subroutine Test_GeTrfTri!矩阵分解成LU,并求逆阵
   use lapack95
   implicit none
   real::a(3,3)=(/1,3,1,1,2,3,3,2,1/)
   integer::ipiv(3),info=1
   !integer::work(1)=3
   write(*,*) a(1,:)
   write(*,*) a(2,:)
   write(*,*) a(3,:)
   ipiv=3
   !call sgetrf(3,3,a,3,ipiv,info)!Fortran77接口矩阵分解成LU
   !call sgetri(3,a,3,ipiv,work(1),3,info)!求逆阵
   call getrf(a,ipiv=ipiv)!Fortran95简化接口
   !call getrf(a)
   call getri(a,ipiv=ipiv,info=info)
   write(*,*) info
   write(*,*) ipiv
   write(*,*) a(1,:)
   write(*,*) a(2,:)
   write(*,*) a(3,:)
End subroutine

你帮我运行一下看看，数组a来源于彭国伦的书第455页。

FLY

fcode 发表于 2014-7-23 15:11
[mw_shl_code=fortran,true]call getrf(a,ipiv,info)!矩阵分解成LU
!call getrf(a)!矩阵分解成LU
call getr ...

我继续试了几个，如果不在分解时用把ipiv的数据取出来，那么求逆阵都是不正确的，即使info的值为0。你帮我检查一下我的代码吧，谢啦

fcode

getrf 函数进行LU分解，使用了部分选主元。也就是会进行部分的行变换。变换结果会返回给 ipiv，在 getri 中使用。

如果你不返回，那么 getri 中的就是错误的。自然结果就不对了。（其实只是进行了行变换（交换）而已）

[Fortran] 纯文本查看 复制代码

ipiv=3
ipiv(1)=2 !// 加上这行就对了
call getrf(a)
call getri(a,ipiv=ipiv,info=info)

从上面的代码可以看出，ipiv不一定非要从getrf里返回。如果你事先知道，你自己生成也可以。所以与 getrf 无关。

FLY

fcode 发表于 2014-7-23 23:47
getrf 函数进行LU分解，使用了部分选主元。也就是会进行部分的行变换。变换结果会返回给 ipiv，在 getri 中 ...

明白了，给力。也就是说，很多情况下计算逆阵时，还是保留getrf计算出的ipiv值比较保险，因为如果自己知道ipiv值也不用来编程算了。

fcode

是的，就是如此。如果你为了求逆，那么最好保留。如果只是做LU分解，则可以不必保留。

PS：对管理员和版主，你不必评分，会从你的分数中扣除。管理员和版主不缺分。