时间戳

首页 > 解决方案 > 如果十进制更长，Fortran-逆矩阵结果不同

matrix - 如果十进制更长，Fortran-逆矩阵结果不同

问题描述

我的真实数据是第一次输入，但结果的倒数太大了。当您与第一个和第二个输入进行比较时，它们是相同的数据。只有十进制大小不同。为什么会有不同的结果？因为它们是相同的数据。他们怎么会有不同的结果？您可以看到结果和输入。太奇怪了。

program test

Implicit none

double precision,allocatable,dimension(:,:)         :: A       
double precision,allocatable,dimension(:)           :: WORK
integer ,allocatable,dimension(:)       :: ipiv
integer                                 :: n,info,M
 external     DGETRF,DGETRI
M=8
allocate(A(M,M),WORK(M),IPIV(M))
!!! First Input !!!!
A(1,:)=(/3.740486048842566D-4, 0.0D0, 0.0D0, 4.987315029057229D-5, 0.0D0, 0.0D0, 0.0D0, 0.0D0/)
A(2,:)=(/0.0D0 , 3.740486048842566D-4, 0.0D0, 0.0D0, 4.987315029057229D-5 ,0.0D0 ,0.0D0 ,0.0D0 /)
A(3,:)=(/0.0D0 , 0.0D0 ,3.740486048842566D-4, 0.0D0 ,0.0D0, 4.987315029057229D-5, 0.0D0 ,0.0D0/)
A(4,:)=(/4.987315029057229D-5 ,0.0D0 ,0.0D0 ,6.649753768432517D-6, 0.0D0 ,0.0D0, 0.0D0, 0.0D0 /)
A(5,:)=(/0.0D0 , 4.987315029057229D-5, 0.0D0, 0.0D0 ,6.649753768432517D-6 ,0.0D0 ,0.0D0 ,0.0D0 /)
A(6,:)=(/0.0D0, 0.0D0, 4.987315029057229D-5, 0.0D0 ,0.0D0, 6.649753768432517D-6, 0.0D0 ,0.0D0 /)
A(7,:)=(/0.0D0, 0.0D0 ,0.0D0, 0.0D0 ,0.0D0 ,0.0D0 ,1.499999910593033D-11, 0.0D0 /)
A(8,:)=(/0.0D0 ,0.0D0 ,0.0D0 ,0.0D0 ,0.0D0 ,0.0D0, 0.0D0 ,1.499999910593033D-11 /)
 !!!! Second Input !!!! 
!A(1,:)=(/3.74D-4, 0.0D0, 0.0D0, 4.98D-5, 0.0D0, 0.0D0, 0.0D0, 0.0D0/)
!A(2,:)=(/0.0D0 , 3.74D-4, 0.0D0, 0.0D0, 4.98D-5 ,0.0D0 ,0.0D0 ,0.0D0 /)
!A(3,:)=(/0.0D0 , 0.0D0 ,3.74D-4, 0.0D0 ,0.0D0, 4.98D-5, 0.0D0 ,0.0D0/)
!A(4,:)=(/4.98D-5 ,0.0D0 ,0.0D0 ,6.64D-6, 0.0D0 ,0.0D0, 0.0D0, 0.0D0 /)
!A(5,:)=(/0.0D0 , 4.98D-5, 0.0D0, 0.0D0 ,6.64D-6 ,0.0D0 ,0.0D0 ,0.0D0 /)
!A(6,:)=(/0.0D0, 0.0D0, 4.98D-5, 0.0D0 ,0.0D0, 6.64D-6, 0.0D0 ,0.0D0 /)
!A(7,:)=(/0.0D0, 0.0D0 ,0.0D0, 0.0D0 ,0.0D0 ,0.0D0 ,1.49D-11, 0.0D0 /)
!A(8,:)=(/0.0D0 ,0.0D0 ,0.0D0 ,0.0D0 ,0.0D0 ,0.0D0, 0.0D0 ,1.49D-11 /)


call DGETRF(M,M,A,M,IPIV,info)
if(info .eq. 0) then
Print *,'succeded'
else
Print *,'failed'
end if

call DGETRI(M,A,M,IPIV,WORK,M,info)
if(info .eq. 0) then
 Print *,'succeded'
else
Print *,'failed'
end if
Print *,A

deallocate(A,IPIV,WORK)

end 
!!!!! Second Input Result
!1.0e+10 *
! 0.0002     0       0   -0.0015       0      0        0   0
!     0      0.0002  0       0       -0.0015  0        0   0
!     0      0    0.0002     0         0     -0.0015   0   0
! -0.0015    0       0     0.0113      0      0        0   0
!     0     -0.0015  0       0       0.0113   0        0   0
!     0      0   -0.0015     0         0    0.0113     0   0
!     0      0       0       0         0      0     6.7114 0
!     0      0       0       0         0      0        0   6.7114

!!! First Input Result
!   1.0e+21 *

!-0.0238         0         0    0.1783         0         0         0         0
!     0   -0.0238         0         0    0.1783         0         0         0
!     0         0    0.0000         0         0   -0.0000         0         0
! 0.1783         0         0   -1.3375         0         0         0         0
!     0    0.1783         0         0   -1.3375         0         0         0
!     0         0   -0.0000         0         0    0.0000         0         0
!     0         0         0         0         0         0    0.0000         0
!     0         0         0         0         0         0         0    0.0000

标签： matrixfortranfortran90matrix-inversefortran95

解决方案

创建矩阵逆不是一个难题。我将您之前的示例转换为使用一种简单的方法，该方法基于带有阴影单位矩阵的高斯消元法，适用于大多数情况。附加的程序会反转您之前的对称矩阵，而不需要对行进行旋转。它不需要“黑匣子”。

用不同的系数得到不同的结果并不奇怪。由于输入值的明显微小变化导致结果显着变化，显示您正在使用的方程关系的灵敏度和可能较差的条件。

https://www.dropbox.com/s/ssotjx45yrz5sf9/dgetri.f90?dl=0

关于“第一个输入”的附加响应

https://www.dropbox.com/s/hximfoin977rmov/dgetri_piv4.f90?dl=0

这个最新的链接 (16-6) 包含两个数据集。在“第一个输入”中，您的方程式基本上是 4:6 行，即 1:3 / 7.5 + small_noise 行。

这个最新的代码示例在矩阵求逆期间和之后都进行了准确性检查。测试期间检查行更改是否正确，而后检查是“AA^-1 - I”和“A - (A^-1)^-1”，这更好地表明准确性差。

有趣的是，“第二输入”（有更多噪音）报告了相当准确的结果。无法获得 8 字节实数的逆需要一个相当人为的矩阵！类似地，随机数导出的系数示例显示出良好的准确性。

这些例子表明，我提出的准确性测试并不总是能识别定义不明确的方程关系。您检查逆以识别值的大变化也很有用。

鉴于方程的定义方式，我不确定您想要的结果是什么。

推荐阅读