Fortran Coder

查看: 97|回复: 7

[数学库] mkl矩阵求逆计算结果与MATLAB计算结果不符,求助!!

[复制链接]

5

帖子

2

主题

0

精华

新人

F 币
27 元
贡献
12 点
发表于 2022-6-23 17:30:44 | 显示全部楼层 |阅读模式
本帖最后由 jkl1102925008 于 2022-6-23 21:13 编辑

存在待求逆矩阵 a=
1.46989040000000e-09        8.08439760000000e-10        -9.40729890000000e-09
8.08439760000000e-10        4.44641850000000e-10        -5.17401410000000e-09
-9.40729890000000e-09        -5.17401410000000e-09        6.02067090000000e-08



subroutine Kalman_Gain_calc
    use lapack95


    real(kind=8) :: a(3,3),b(3)
    integer :: ipiv(3)
    integer :: n,info1,info2
    b=0.0
    ipiv=0
    a=0.0
    n=3

    a=st_sim_cov
    call sgetrf(n,n,a,n,ipiv,info1)
    call sgetri(n,a,n,ipiv,b,n,info2)
end subroutine





计算结果 a=
   54127.1604943661       2.562316454118586E-315 -9.407298939834163E-009
   22340.6659669089       4.446417947376324E-010 -5.174014106046343E-009
  0.000000000000000E+000 -5.174014106046343E-009  6.020670895168223E-008

而MATLAB INV(a) 计算结果为
-6.61095895614289e+15        2.60655716871615e+15        -808961396190492
2.60655716871615e+15        1.05500250710316e+17        9.47370245921866e+15
-808961396190492        9.47370245921866e+15        687746102360067

大佬们帮我看看,是我哪里搞错了吗,问题出在哪里,求求了,嘤嘤嘤!

-----------------------------------------------------------------------------------------
使用sgetrf+sgetri计算的结果是错误的。。。



70

帖子

2

主题

1

精华

专家

Vim

F 币
594 元
贡献
301 点

规矩勋章

发表于 2022-6-23 21:38:50 | 显示全部楼层
精度不够,使用双精度 dgetrf,dgetri

5

帖子

2

主题

0

精华

新人

F 币
27 元
贡献
12 点
 楼主| 发表于 2022-6-23 22:06:07 | 显示全部楼层
Transpose 发表于 2022-6-23 21:38
精度不够,使用双精度 dgetrf,dgetri

使用 dgetrf+dgetri后 获得 c=
-7.566046942054234E+015 -1.441152316827714E+015 -1.306043946963708E+015
-1.441152330078778E+015  1.095275433704068E+017  9.187343291859882E+015
-1.306043948102471E+015  9.187343293930360E+015   585467958061173.

但是 a*c =
   1.00000000000000       0.000000000000000E+000  0.000000000000000E+000
  0.000000000000000E+000  0.000000000000000E+000  0.250000000000000
  0.000000000000000E+000   64.0000000000000       0.000000000000000E+000

而不是E,计算结果仍有错误欸

92

帖子

0

主题

0

精华

大师

F 币
633 元
贡献
281 点

规矩勋章元老勋章新人勋章

发表于 2022-6-24 08:22:51 | 显示全部楼层
本帖最后由 胡文刚 于 2022-6-24 08:23 编辑

我用以下代码测试的结果是正确的。

  0.999999997206032      -1.490116119384766E-008  9.313225746154785E-010
  0.000000000000000E+000  0.999999992549419       0.000000000000000E+000
  1.490116119384766E-008  1.192092895507812E-007  0.999999992549419

[Fortran] 纯文本查看 复制代码
Program Main
  use lapack95
  implicit none
  integer , parameter :: N = 3
  real(kind=8) :: a(N,N),b(N),t(N,N)
  integer :: ipiv(N) = 0,info1,info2
  a=reshape([1.46989040000000d-09, 8.08439760000000d-10,-9.40729890000000d-09,&
             8.08439760000000d-10, 4.44641850000000d-10,-5.17401410000000d-09,&
            -9.40729890000000d-09,-5.17401410000000d-09, 6.02067090000000d-08],[N,N])
  t=a
  call dgetrf(n,n,a,n,ipiv,info1)
  call dgetri(n,a,n,ipiv,b,n,info2)
  write(*,*) matmul(a,t)
End Program Main

天之道,损有余而补不足

5

帖子

2

主题

0

精华

新人

F 币
27 元
贡献
12 点
 楼主| 发表于 2022-6-24 09:43:39 | 显示全部楼层
本帖最后由 jkl1102925008 于 2022-6-24 10:22 编辑

一样的代码,今天在我重新启动电脑后
c=
-7.566046942054234E+015 -1.441152316827714E+015 -1.306043946963708E+015
-1.441152330078778E+015  1.095275433704068E+017  9.187343291859882E+015
-1.306043948102471E+015  9.187343293930360E+015   585467958061173.

a*c=
   1.00000000000000       0.000000000000000E+000  0.000000000000000E+000
  0.000000000000000E+000  0.999999992549419       0.000000000000000E+000
  0.000000000000000E+000  4.656612873077393E-010   1.00000000000000

然而在我进行循环的时候。。。。再一次出错。 单精度和双精度相乘是否。。。可能出现错误呢。。
第一轮:
输入a=
  1.4698904E-09  8.0843976E-10 -9.4072989E-09
  8.0843976E-10  4.4464185E-10 -5.1740141E-09
-9.4072989E-09 -5.1740141E-09  6.0206709E-08

输出逆矩阵c=
-7.566046942054234E+015 -1.441152316827714E+015 -1.306043946963708E+015
-1.441152330078778E+015  1.095275433704068E+017  9.187343291859882E+015
-1.306043948102471E+015  9.187343293930360E+015   585467958061173.

a*c=
   1.00000000000000       0.000000000000000E+000  0.000000000000000E+000
  0.000000000000000E+000  0.999999992549419       0.000000000000000E+000
  0.000000000000000E+000  4.656612873077393E-010   1.00000000000000
(我尝试在matlab中进行a*c 最终结果为)
a*c=
    1.1864    0.4574    0.0684
   -0.0095    0.9327   -0.0073
   -0.3732    1.0487    1.0318

第二轮:
输入a=
  2.1166422E-09  5.2916055E-10  1.1994306E-08
  5.2916055E-10  1.3229014E-10  2.9985765E-09
  1.1994306E-08  2.9985765E-09  6.7967733E-08

输出逆矩阵c=
  1.199430599285733E-008  2.998576498214334E-009  6.796773277528700E-008
  0.176470584968382      -2.793420560336152E-026  4.310277579995551E-016
  4.411764624209553E-002  0.250000000000000       0.000000000000000E+000

a*c=
  8.422001670951175E-016  2.105500417737794E-016  4.772467535878760E-015
  3.735250862938224E-010  9.338127157345560E-011  2.116642194849335E-009
  2.256714087515359E-010  5.641785218788397E-011  1.278804673265910E-009



第三轮:
输入a=
  1.4698904E-09  1.4698905E-10  5.8795617E-09
  1.4698905E-10  1.4698904E-11  5.8795618E-10
  5.8795617E-09  5.8795618E-10  2.3518247E-08

输出逆矩阵c=
-6.127803175762061E+024   16.0000000000000       1.531950793940515E+024
  1.157130282314151E+016 -7.585009868576445E+017  1.606969932372038E+016
  1.531950793651233E+024  1.896252502950575E+016 -3.829876988868712E+023

a*c=
  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
  0.000000000000000E+000   1.00000000000000       0.000000000000000E+000
  0.500000000000000       3.125000000000000E-002   2.00000000000000

是否能在fortran中使用matlab函数呢


92

帖子

0

主题

0

精华

大师

F 币
633 元
贡献
281 点

规矩勋章元老勋章新人勋章

发表于 2022-6-24 11:27:39 | 显示全部楼层
你都用双精度就好了。

如果你第一次运行是好的,循环中出了错。请优先怀疑是自己的问题。

除非fortran中很难实现,没有现成的东西,否则不要在fortran中调用matlab。麻烦,效率低,依赖度高。
天之道,损有余而补不足

19

帖子

2

主题

0

精华

入门

F 币
82 元
贡献
39 点
发表于 2022-6-24 17:13:29 | 显示全部楼层
就是精度的问题,你用real(kind=8)表示的就是双精度实型,你下面函数却用sgetrf,sgetri。这两个函数是对单精度的用的。你用法有瑕疵

5

帖子

2

主题

0

精华

新人

F 币
27 元
贡献
12 点
 楼主| 发表于 5 天前 | 显示全部楼层
zjk0112 发表于 2022-6-24 17:13
就是精度的问题,你用real(kind=8)表示的就是双精度实型,你下面函数却用sgetrf,sgetri。这两个函数是对单 ...

后修改后仍出现问题。后来我直接matlab打包exe后在fortran中调用,需要安装matlab runtime,且速度较慢,但可以接受,并且由于matlab中有许多现成代码可以直接使用,感觉还行。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2022-7-5 21:47

Powered by Tencent X3.4

© 2013-2022 Tencent

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