论坛 › 编程工具交流 › mkl矩阵求逆计算结果与MATLAB计算结果不符，求助！！

jkl1102925008 发表于 2022-6-23 17:30:44

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

本帖最后由 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-0096.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+0151.095275433704068E+0179.187343291859882E+015
-1.306043948102471E+0159.187343293930360E+015   585467958061173.

但是 a*c =
   1.00000000000000       0.000000000000000E+0000.000000000000000E+000
0.000000000000000E+0000.000000000000000E+0000.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-0089.313225746154785E-010
0.000000000000000E+0000.999999992549419       0.000000000000000E+000
1.490116119384766E-0081.192092895507812E-0070.999999992549419

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],)
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+0151.095275433704068E+0179.187343291859882E+015
-1.306043948102471E+0159.187343293930360E+015   585467958061173.

a*c=
   1.00000000000000       0.000000000000000E+0000.000000000000000E+000
0.000000000000000E+0000.999999992549419       0.000000000000000E+000
0.000000000000000E+0004.656612873077393E-010   1.00000000000000

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

输出逆矩阵c=
-7.566046942054234E+015 -1.441152316827714E+015 -1.306043946963708E+015
-1.441152330078778E+0151.095275433704068E+0179.187343291859882E+015
-1.306043948102471E+0159.187343293930360E+015   585467958061173.

a*c=
   1.00000000000000       0.000000000000000E+0000.000000000000000E+000
0.000000000000000E+0000.999999992549419       0.000000000000000E+000
0.000000000000000E+0004.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-095.2916055E-101.1994306E-08
5.2916055E-101.3229014E-102.9985765E-09
1.1994306E-082.9985765E-096.7967733E-08

输出逆矩阵c=
1.199430599285733E-0082.998576498214334E-0096.796773277528700E-008
0.176470584968382      -2.793420560336152E-0264.310277579995551E-016
4.411764624209553E-0020.250000000000000       0.000000000000000E+000

a*c=
8.422001670951175E-0162.105500417737794E-0164.772467535878760E-015
3.735250862938224E-0109.338127157345560E-0112.116642194849335E-009
2.256714087515359E-0105.641785218788397E-0111.278804673265910E-009



第三轮：
输入a=
1.4698904E-091.4698905E-105.8795617E-09
1.4698905E-101.4698904E-115.8795618E-10
5.8795617E-095.8795618E-102.3518247E-08

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

a*c=
0.000000000000000E+0000.000000000000000E+0000.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计算结果不符，求助！！