Fortran Coder

查看: 9641|回复: 5
打印 上一主题 下一主题

[通用算法] MKL函数库中ZGETRF和ZGETRI函数的使用

[复制链接]

68

帖子

25

主题

0

精华

专家

F 币
321 元
贡献
197 点
跳转到指定楼层
楼主
发表于 2019-11-8 03:12:43 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
大家好,

我想向大家请教一个关于求逆矩阵的问题。我的矩阵是52行,52列的复数矩阵。我想求它的逆矩阵,我想用MKL函数库里的ZGETRFZGEFTI函数来求。为此我先做了一个测试,分别各自求了一个22列和一个33列的复数矩阵的逆矩阵。我把我的代码贴在下面了。
[Fortran] 纯文本查看 复制代码
PROGRAM TE
IMPLICIT NONE
INTEGER :: i, j, dime                                      !Ingeter for loop and order for complex matrix
DOUBLE COMPLEX, ALLOCATABLE :: a(:,:)       !Original complex matrix
DOUBLE COMPLEX, ALLOCATABLE :: b(:,:)       !Inversed complex matrix
DOUBLE COMPLEX, ALLOCATABLE :: c(:,:)       !Multiplication of the original andinversed complex matrices

!Variables used in LU decomposition subroutine (ZGETRF)
INTEGER                                     :: nu_r        !Number ofrow in target matrix
INTEGER                                     :: nu_c        !Number ofcolumn in target matrix
INTEGER                                     :: lu_lda      !Leadingdimension of target matrix
INTEGER, ALLOCATABLE               :: lu_ipiv(:)  !Arraywith dimension of .GE. MIN(nu_r,nu_c)
INTEGER                                     :: lu_info     !Judgement(0-successful;<0-error and infor = number, then
                                                                    !'1-number'th arguement hasillegal value;>0 and info =
                                                                    !number - LU matrix(number,number) is zero)

!Variables used in inverse subroutine (ZGETRI)
INTEGER                                     :: nu_o        !Order ofLU matrix
DOUBLE COMPLEX, ALLOCATABLE :: in_work(:) !Workspacearray with dimension of .GE. MAX(1,in_lwork)
INTEGER                      :: in_lda                     !Leadingdimension of LU decomposed matrix .GE. MAX(1,nu_o)
INTEGER                      :: in_lwork                  !Size ofin_work array with value .GE. nu_o
INTEGER                      :: in_info                    !Judgement(0-successful;<0-error and infor = number, then
                                                                    !'1-number'th arguement hasillegal value;>0 and info =
                                                                    !number; then, LU matrix(number,number) is zero and inversion
                                                                    !cannot be accomplished)

dime = 3                                                       !Thisparameter is used to control the size of the complex matrix
ALLOCATE (a(dime,dime))
ALLOCATE (b(dime,dime))
ALLOCATE (c(dime,dime))

DO i = 1, dime, 1
   DO j = 1, dime, 1
      a(i,j) =CMPLX(DBLE(i), DBLE(j))
   END DO
END DO

b = a

OPEN (UNIT=3, FILE='te_in.dat', STATUS='UNKNOWN')
WRITE (UNIT=3, FMT='(A1)') 'a'
DO i = 1, dime, 1
   DO j = 1, dime, 1
      WRITE (UNIT=3,FMT=*) a(i,j)
   END DO
END DO
WRITE (UNIT=3, FMT=*)

!Initialising parameters for LU decomposition subroutine(ZGETRF)
nu_r = dime
nu_c = dime
lu_lda = dime
ALLOCATE (lu_ipiv(dime))
lu_ipiv = 0
!

!Initialising parameters for inverse subroutine (ZGETRI)
nu_o = dime
in_lwork = dime
ALLOCATE (in_work(in_lwork))
in_work = (0.0d0, 0.0d0)
in_lda = dime
!

CALL ZGETRF(nu_r,nu_c,b,lu_lda,lu_ipiv,lu_info)
CALL ZGETRI(nu_o,b,in_lda,lu_ipiv,in_work,in_lwork,in_info)

WRITE (UNIT=3, FMT='(A1)') 'b'
DO i = 1, dime, 1
   DO j = 1, dime, 1
      WRITE (UNIT=3,FMT=*) b(i,j)
   END DO
END DO
WRITE (UNIT=3, FMT=*)

c = MATMUL(a,b)

WRITE (UNIT=3, FMT='(A1)') 'c'
DO i = 1, dime, 1
   DO j = 1, dime, 1
      WRITE (UNIT=3,FMT=*) c(i,j)
   END DO
END DO
WRITE (UNIT=3, FMT=*)

DEALLOCATE (lu_ipiv)
DEALLOCATE (a)
DEALLOCATE (b)
DEALLOCATE (c)

CLOSE (UNIT=3)
STOP
END PROGRAM TE
当测试22列的复数矩阵时,结果是正确的,因为原矩阵和它的逆矩阵乘积是一个单位矩阵。结果如下。
a
(1.00000000000000,1.00000000000000)
(1.00000000000000,2.00000000000000)
(2.00000000000000,1.00000000000000)
(2.00000000000000,2.00000000000000)

b
(-2.00000000000000,2.00000000000000)
(2.00000000000000,-1.00000000000000)
(1.00000000000000,-2.00000000000000)
(-1.00000000000000,1.00000000000000)

c
(1.00000000000000,-2.220446049250313E-016)
(0.000000000000000E+000,2.220446049250313E-016)
(0.000000000000000E+000,4.440892098500626E-016)
(1.00000000000000,0.000000000000000E+000)

但是当测试33列的矩阵时,结果却是错误的,原矩阵和它的逆矩阵乘积不是单位矩阵。结果如下。
a
(1.00000000000000,1.00000000000000)
(1.00000000000000,2.00000000000000)
(1.00000000000000,3.00000000000000)
(2.00000000000000,1.00000000000000)
(2.00000000000000,2.00000000000000)
(2.00000000000000,3.00000000000000)
(3.00000000000000,1.00000000000000)
(3.00000000000000,2.00000000000000)
(3.00000000000000,3.00000000000000)

b
(723205779577744.,262983919846453.)
(-1.446411559155488E+015,-525967839692903.)
(723205779577744.,262983919846451.)
(-1.446411559155489E+015,-525967839692904.)
(2.892823118310976E+015,1.051935679385808E+015)
(-1.446411559155487E+015,-525967839692903.)
(723205779577745.,262983919846452.)
(-1.446411559155488E+015,-525967839692905.)
(723205779577743.,262983919846452.)

c
(0.437500000000000,0.000000000000000E+000)
(1.25000000000000,0.000000000000000E+000)
(-0.250000000000000,0.000000000000000E+000)
(0.125000000000000,0.000000000000000E+000)
(0.500000000000000,0.000000000000000E+000)
(0.750000000000000,0.000000000000000E+000)
(-0.500000000000000,0.000000000000000E+000)
(0.000000000000000E+000,-1.00000000000000)
(1.50000000000000,0.500000000000000)

我有测试了更高阶的复数矩阵,发现结果也是错误的(原矩阵和它的逆矩阵乘积不是单位矩阵)。我用来编译FORTRAN程序的脚本文件如下。
#!/bin/bash
module switch PrgEnv-cray PrgEnv-intel
#module load intel
#----------------------------------------------------------#
ifort -O3 te.f90 \
  -Wl,--start-group${MKLROOT}/lib/intel64/libmkl_intel_lp64.a \
                   ${MKLROOT}/lib/intel64/libmkl_core.a \
                   ${MKLROOT}/lib/intel64/libmkl_sequential.a \
  -Wl,--end-group-lpthread -lm –ldl

我没有找出其中的错误。能否麻烦大家帮我看下,我的代码中哪里出了问题呢?要如何改正呢?

另外,有没有专门计算复数矩阵的逆矩阵的程序,或算法推荐给我呢?多谢大家的帮助啦。
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

834

帖子

2

主题

0

精华

大宗师

F 币
3922 元
贡献
2332 点
沙发
发表于 2019-11-10 14:36:24 | 只看该作者
本帖最后由 li913 于 2019-11-10 14:38 编辑

放宽心,代码没问题,是误差问题。A不是对角占优的,导致逆矩阵B某些值特别大(15次方),B的误差可能在0.1-1之间,最终导致C非单位阵。你试试让对角元素大一些,比如 if(i==j) a(i,j)=2*a(i,j),结果就ok。
另外,推荐使用f95接口,更加简单好使。CALL GETRF(b,lu_ipiv,lu_info)
call getri( b, lu_ipiv,in_info)

68

帖子

25

主题

0

精华

专家

F 币
321 元
贡献
197 点
板凳
 楼主| 发表于 2019-11-11 03:05:01 | 只看该作者
li913 发表于 2019-11-10 14:36
放宽心,代码没问题,是误差问题。A不是对角占优的,导致逆矩阵B某些值特别大(15次方),B的误差可能在0.1-1 ...

谢谢你的回复和建议。

我用的那个服务器上没有f95版本的程序。我又没有管理员权限,无法安装。

我在看徐士良的《FORTRAN算法常用程序集》,里面有介绍计算复数矩阵的逆矩阵的算法。我打算用那个程序算下逆矩阵。

非常感谢你的建议。

834

帖子

2

主题

0

精华

大宗师

F 币
3922 元
贡献
2332 点
地板
发表于 2019-11-12 13:27:15 | 只看该作者
1、mkl库自带f95接口,无需额外安装;
2、徐士良程序集,你自己写代码,先不说正确性,效率肯定不如mkl。仅可作为验证手段。

68

帖子

25

主题

0

精华

专家

F 币
321 元
贡献
197 点
5#
 楼主| 发表于 2019-11-13 05:52:13 | 只看该作者
li913 发表于 2019-11-12 13:27
1、mkl库自带f95接口,无需额外安装;
2、徐士良程序集,你自己写代码,先不说正确性,效率肯定不如mkl。仅 ...

非常感谢你的解释。

我手算了我程序里的复数矩阵的行列式的值,发现3乘3的复数矩阵的行列式的值是(0,0),也就是一个奇异矩阵,所以我现在可以理解3乘3的复数矩阵没有逆矩阵了。

可我不理解的是ZGETRF和ZGETRI两个函数里的lu_info,in_info的参数值都是0,就是说这两个函数并没有报错,认为原3乘3的复数矩阵没有问题。

我用这两个函数计算了我的代码里的52乘52的复数矩阵的逆矩阵,发现逆矩阵与原矩阵相乘,结果也不是单位矩阵,所以逆矩阵出了问题。我想这可能与编译器或者服务器系统的精度有关。就是原复数矩阵的行列式的值虽不为零(理论上应该有逆矩阵),但却是个非常小的值。这也许会使得函数求得的逆矩阵是一个不合理的结果。所以逆矩阵与原矩阵相乘,结果不是单位矩阵。

不知道你是否也遇到过这类的问题呢?又是如何解决的呢?非常感谢。

834

帖子

2

主题

0

精华

大宗师

F 币
3922 元
贡献
2332 点
6#
发表于 2019-11-14 14:55:46 | 只看该作者
没有详细研究过。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-11-1 08:07

Powered by Tencent X3.4

© 2013-2024 Tencent

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