Fortran Coder

查看: 14658|回复: 7

[数学库] MKL中getrf省略参数ipiv的问题

[复制链接]

31

帖子

5

主题

0

精华

入门

F 币
105 元
贡献
54 点
发表于 2014-7-23 11:07:56 | 显示全部楼层 |阅读模式
如果函数getrf参数中没有ipiv作为返回值放入getri中的ipiv里,计算出的结果是错误的,即使getrf中的ipiv是可选的。这是为什么?

31

帖子

5

主题

0

精华

入门

F 币
105 元
贡献
54 点
 楼主| 发表于 2014-7-23 11:16:41 | 显示全部楼层
vvt 发表于 2014-7-23 10:55
1. F77 语法本身,对接口的检查就不严格。建议使用 F95 接口
2. 是的,只有第一个参数是必选的。可选参数的 ...

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

1947

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1295 元
贡献
545 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2014-7-23 15:11:08 | 显示全部楼层
[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
你再检查一下看是不是别的问题吧

31

帖子

5

主题

0

精华

入门

F 币
105 元
贡献
54 点
 楼主| 发表于 2014-7-23 20:57:34 | 显示全部楼层
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页。

31

帖子

5

主题

0

精华

入门

F 币
105 元
贡献
54 点
 楼主| 发表于 2014-7-23 21:10:27 | 显示全部楼层
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。你帮我检查一下我的代码吧,谢啦

1947

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1295 元
贡献
545 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2014-7-23 23:47:37 | 显示全部楼层
getrf 函数进行LU分解,使用了部分选主元。也就是会进行部分的行变换。变换结果会返回给 ipiv,在 getri 中使用。

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

[Fortran] 纯文本查看 复制代码
ipiv=3
ipiv(1)=2 !// 加上这行就对了
call getrf(a)
call getri(a,ipiv=ipiv,info=info)


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

评分

参与人数 1F 币 +5 贡献 +5 收起 理由
FLY + 5 + 5 很给力!

查看全部评分

31

帖子

5

主题

0

精华

入门

F 币
105 元
贡献
54 点
 楼主| 发表于 2014-7-23 23:56:52 | 显示全部楼层
fcode 发表于 2014-7-23 23:47
getrf 函数进行LU分解,使用了部分选主元。也就是会进行部分的行变换。变换结果会返回给 ipiv,在 getri 中 ...

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

1947

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1295 元
贡献
545 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2014-7-24 00:42:12 | 显示全部楼层
是的,就是如此。如果你为了求逆,那么最好保留。如果只是做LU分解,则可以不必保留。

PS:对管理员和版主,你不必评分,会从你的分数中扣除。管理员和版主不缺分。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-3-28 21:55

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表