mkl矩阵求逆计算结果与MATLAB计算结果不符，求助！！

jkl1102925008 · 发表于 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计算的结果是错误的。。。

Transpose · 发表于 2022-6-23 21:38:50

精度不够，使用双精度 dgetrf,dgetri

jkl1102925008 · 发表于 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，计算结果仍有错误欸

胡文刚 · 发表于 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

jkl1102925008 · 发表于 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函数呢

胡文刚 · 发表于 2022-6-24 11:27:39

你都用双精度就好了。

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

除非fortran中很难实现，没有现成的东西，否则不要在fortran中调用matlab。麻烦，效率低，依赖度高。

zjk0112 · 发表于 2022-6-24 17:13:29

就是精度的问题，你用real(kind=8)表示的就是双精度实型，你下面函数却用sgetrf，sgetri。这两个函数是对单精度的用的。你用法有瑕疵

jkl1102925008 · 发表于 2022-6-30 09:29:47

zjk0112 发表于 2022-6-24 17:13
就是精度的问题，你用real(kind=8)表示的就是双精度实型，你下面函数却用sgetrf，sgetri。这两个函数是对单 ...

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

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

规矩勋章

元老勋章

新人勋章