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

jkl1102925008

存在待求逆矩阵 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计算的结果是错误的。。。

jkl1102925008

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

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

zjk0112

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

胡文刚

你都用双精度就好了。

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

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

jkl1102925008

一样的代码，今天在我重新启动电脑后
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函数呢

胡文刚

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

  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

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，计算结果仍有错误欸

Transpose

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