数据溢出

Kieran · 发表于 2020-9-25 12:54:16

大家好，

我写了一段代码。运行后，发现输出结果中出现了NaN类型的数据。我猜测可能是因为数据溢出，就是数据过小，超出了双精度复数的表示范围。但我并不确定是否是这个原因，所以我把我出问题部分的代码粘贴在下面了，麻烦大家帮我看下是否真的是因为数据溢出的错误。如果不是，麻烦大家给些建议是哪里出错了；如果是的话，麻烦大家给些建议如何判断哪些过小的双精度复数型数据，并把它们设置为零吧。多谢大家。

我的出问题部分的代码。

[Fortran] 纯文本查看 复制代码

001!External variables
002INTEGER                       :: vd
003DOUBLE PRECISION, ALLOCATABLE :: kp(:,:)
004INTEGER                       :: nu_wa
005INTEGER                       :: nu_ws
006INTEGER, ALLOCATABLE          :: nd(:)
007DOUBLE PRECISION, ALLOCATABLE :: hr(:,:)
008DOUBLE COMPLEX, ALLOCATABLE   :: tb(:,:)
009DOUBLE COMPLEX, ALLOCATABLE   :: ec(:,:)
010DOUBLE PRECISION, ALLOCATABLE :: ev(:,:)
011DOUBLE COMPLEX, ALLOCATABLE   :: vh(:,:)
012DOUBLE PRECISION              :: fermi
013DOUBLE PRECISION              :: wf
014DOUBLE PRECISION              :: bv
015INTEGER                       :: nv
016DOUBLE PRECISION              :: mv
017INTEGER                       :: nn
018DOUBLE PRECISION              :: dk(2)
019DOUBLE COMPLEX, ALLOCATABLE   :: sr(:,:)
020DOUBLE COMPLEX                :: sr1(3)
021INTEGER                       :: nu_kp
022 
023INTEGER                       :: i, j, k, l, m
024DOUBLE COMPLEX                :: dc(3)
025DOUBLE COMPLEX, ALLOCATABLE   :: ex(:,:)
026DOUBLE PRECISION              :: de
027 
028DOUBLE PRECISION              :: ey
029DOUBLE COMPLEX                :: gr
030DOUBLE COMPLEX                :: ga
031DOUBLE COMPLEX                :: gd(2)
032DOUBLE COMPLEX                :: su
033DOUBLE COMPLEX, ALLOCATABLE   :: mt1(:,:)
034DOUBLE COMPLEX, ALLOCATABLE   :: mt2(:,:)
035DOUBLE COMPLEX, ALLOCATABLE   :: mt3(:,:)
036DOUBLE COMPLEX, ALLOCATABLE   :: mt4(:,:)
037DOUBLE COMPLEX, ALLOCATABLE   :: mt5(:,:)
038DOUBLE COMPLEX                :: mt6(8)
039DOUBLE COMPLEX                :: sp(6)
040DOUBLE PRECISION, ALLOCATABLE :: nm(:)
041DOUBLE COMPLEX                :: oc(3)
042 
043INTEGER                       :: omp_get_num_threads
044INTEGER                       :: num_thread
045INTEGER                       :: omp_get_thread_num
046INTEGER                       :: id_thread
047DOUBLE COMPLEX, ALLOCATABLE   :: sr2(:,:,:)
048!$OMP PARALLEL DO PRIVATE(i, j, k, l, ey, sr1, mt1, mt2, mt3, mt4, mt5, mt6, nm, oc, gr, ga, gd, su)&
049 
050!$OMP             SHARED(ev, ec, ex, nv, nn, mv, wf, de, nu_kp, nu_wa, vd, vh, bv, dk, sr2)
051 
052ALLOCATE (sr2(num_thread,nv+nn,3))
053sr2 = CMPLX(0.0d0, 0.0d0)
054 
055DO i = 1, nv + nn - 2, 1
056   id_thread = omp_get_thread_num() + 1
057 
058   ey = mv - wf + DBLE(i) * de
059   sr1 = CMPLX(0.0d0, 0.0d0)
060 
061   DO j = 1, nu_kp, 1
062      !Normalized factor part
063      DO k = 1, nu_wa, 1
064         mt1(1,:) = DCONJG(ec(1+(j-1)*nu_wa:j*nu_wa,k))
065         mt2(:,1) = ec(1+(j-1)*nu_wa:j*nu_wa,k)
066         nm(k) = REAL(SUM(MATMUL(mt1,mt2)))
067         nm(k) = 1.0d0 / DSQRT(nm(k))
068      END DO
069      IF (vd .EQ. 0) THEN
070          mt5 = vh(1+(j-1)*nu_wa:j*nu_wa,1:nu_wa)
071      ELSE
072          mt5 = vh(1+(j-1)*nu_wa:j*nu_wa,1+nu_wa:2*nu_wa)
073      END IF
074 
075      oc = CMPLX(0.0d0, 0.0d0)
076      DO k = 1, nu_wa, 1
077         mt1(1,:) = DCONJG(ec(1+(j-1)*nu_wa:j*nu_wa,k))
078         mt2(:,1) = ec(1+(j-1)*nu_wa:j*nu_wa,k)
079         gr = CMPLX (1.0d0, 0.0d0) / CMPLX(ey - ev(j,k), bv)
080         ga = CMPLX (1.0d0, 0.0d0) / CMPLX(ey - ev(j,k), 0.0d0 - bv)
081         gd(2) = (gr * gr + ga * ga)
082         DO l = 1, nu_wa, 1
083            mt3(1,:) = DCONJG(ec(1+(j-1)*nu_wa:j*nu_wa,l))
084            mt4(:,1) = ec(1+(j-1)*nu_wa:j*nu_wa,l)
085            gr = CMPLX (1.0d0, 0.0d0) / CMPLX(ey - ev(j,l), bv)
086            ga = CMPLX (1.0d0, 0.0d0) / CMPLX(ey - ev(j,l), 0.0d0 - bv)
087            gd(1) = (gr - ga)
088            sp = CMPLX(0.0d0, 0.0d0)
089            DO m = 1, nu_wa / 2, 1
090               mt6(1) = mt2(m,1)
091               mt6(2) = mt2(m+nu_wa/2,1)
092               mt6(3) = mt3(1,m)
093               mt6(4) = mt3(1,m+nu_wa/2)
094               mt6(5) = mt1(1,m)
095               mt6(6) = mt1(1,m+nu_wa/2)
096               mt6(7) = mt4(m,1)
097               mt6(8) = mt4(m+nu_wa/2,1)
098               sp(1) = sp(1) + mt6(4) * mt6(1) + mt6(3) * mt6(2)
099               sp(2) = sp(2) + dc(1) * (mt6(4) * mt6(1) - mt6(3) * mt6(2))
100               sp(3) = sp(3) + mt6(3) * mt6(1) - mt6(4) * mt6(2)
101               sp(4) = sp(4) + mt6(6) * mt6(7) + mt6(5) * mt6(8)
102               sp(5) = sp(5) + dc(1) * (mt6(6) * mt6(7) - mt6(5) * mt6(8))
103               sp(6) = sp(6) + mt6(5) * mt6(7) - mt6(6) * mt6(8)
104            END DO
105            su = (SUM(MATMUL(mt1,MATMUL(mt5,mt4)))*sp(1)-sp(4)*SUM(MATMUL(mt3,MATMUL(mt5,mt2))))/CMPLX(2.0d0, 0.0d0)
106            oc(1) = oc(1)+su*gd(1)*gd(2)*CMPLX(de, 0.0d0)*nm(k)*nm(k)*nm(l)*nm(l)*dk(1)*dk(2)
107            su = (SUM(MATMUL(mt1,MATMUL(mt5,mt4)))*sp(2)-sp(5)*SUM(MATMUL(mt3,MATMUL(mt5,mt2))))/CMPLX(2.0d0, 0.0d0)
108            oc(2) = oc(2)+su*gd(1)*gd(2)*CMPLX(de, 0.0d0)*nm(k)*nm(k)*nm(l)*nm(l)*dk(1)*dk(2)
109            su = (SUM(MATMUL(mt1,MATMUL(mt5,mt4)))*sp(3)-sp(6)*SUM(MATMUL(mt3,MATMUL(mt5,mt2))))/CMPLX(2.0d0, 0.0d0)
110            oc(3) = oc(3)+su*gd(1)*gd(2)*CMPLX(de, 0.0d0)*nm(k)*nm(k)*nm(l)*nm(l)*dk(1)*dk(2)
111         END DO
112      END DO
113      sr1(1) = sr1(1) + ex(i,2) * oc(3) - ex(i,3) * oc(2)
114      sr1(2) = sr1(2) + ex(i,3) * oc(1) - ex(i,1) * oc(3)
115      sr1(3) = sr1(3) + ex(i,1) * oc(2) - ex(i,2) * oc(1)
116write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)!就是在这里输出的双精度复数型数据，有些数据显示为NaN
117write(unit=*,fmt=*)'T_y_interpolating',ex(i,3) * oc(1) - ex(i,1) * oc(3)!就是在这里输出的双精度复数型数据，有些数据显示为NaN
118write(unit=*,fmt=*)'T_z_interpolating',ex(i,1) * oc(2) - ex(i,2) * oc(1)!就是在这里输出的双精度复数型数据，有些数据显示为NaN
119   END DO
120 
121   sr2(id_thread,1+i,1) = sr1(1)
122   sr2(id_thread,1+i,2) = sr1(2)
123   sr2(id_thread,1+i,3) = sr1(3)
124END DO
125!$OMP END PARALLEL DO

输出的部分结果如下。
T_z_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_y_interpolating (-7.844408716089390E-318,-1.627452237401066E-320)
T_z_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_x_interpolating (NaN,NaN)
T_y_interpolating (NaN,NaN)
T_z_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_x_interpolating (-4.748344416916300E-003,-2.184128388170520E-010)
T_y_interpolating (3.377560741471222E-008,4.832304342329601E-011)
T_y_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_z_interpolating (7.890054335804031E-008,-2.098265883494228E-010)
T_x_interpolating (NaN,NaN)
T_y_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_z_interpolating (NaN,NaN)
T_x_interpolating (3.764965159798266E-007,7.032141314956383E-022)
T_x_interpolating (1.435727593204141E-317,-1.963519148161803E-317)
T_y_interpolating (9.405725326138147E-319,9.193968790330589E-319)
T_z_interpolating (3.839154229161444E-006,2.624325922453405E-010)
T_x_interpolating (NaN,NaN)
T_y_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_z_interpolating (NaN,NaN)
T_x_interpolating (-0.166424638577942,7.705225622188310E-013)
T_y_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_z_interpolating (5.576638611415991E-002,1.704192603919102E-011)
T_x_interpolating (NaN,NaN)
T_y_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_z_interpolating (NaN,NaN)
T_x_interpolating (NaN,NaN)
T_y_interpolating (0.000000000000000E+000,0.000000000000000E+000)
T_z_interpolating (NaN,NaN)
T_x_interpolating (9.167435739615091E-005,-3.946845693068038E-020)
T_y_interpolating (0.000000000000000E+000,0.000000000000000E+000)

麻烦大家给些建议吧。谢谢啦。

necrohan · 发表于 2020-9-26 07:46:12

write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)
这里用的变量没有初始化

Kieran · 发表于 2020-9-26 14:10:55

necrohan 发表于 2020-9-26 07:46
write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)
这里用的变量没有初始化 ...

谢谢你的回复。

我仅仅只是粘贴了出现问题部分的代码，ex变量在前面的代码部分已经赋值过了，他们也是复数型的数组。

我有尝试在含有注释的代码（write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)!就是在这里输出的双精度复数型数据，有些数据显示为NaN）前面加入判断数据溢出的IF语句，如下面所示。

[Fortran] 纯文本查看 复制代码

01      IF (ABS(ex(i,2)*oc(3)-ex(i,3)*oc(2)) <= 1.17549435d-38) THEN
02         sr1(1) = sr1(1) + CMPLX(0.0d0, 0.0d0)
03      ELSE
04         sr1(1) = sr1(1) + ex(i,2)*oc(3)-ex(i,3)*oc(2)
05      END IF
06      IF (ABS(ex(i,3)*oc(1)-ex(i,1)*oc(3)) <= 1.17549435d-38) THEN
07         sr1(2) = sr1(2) + CMPLX(0.0d0, 0.0d0)
08      ELSE
09         sr1(2) = sr1(2) + ex(i,3)*oc(1)-ex(i,1)*oc(3)
10      END IF
11      IF (ABS(ex(i,1)*oc(2)-ex(i,2)*oc(1)) <= 1.17549435d-38) THEN
12         sr1(3) = sr1(3) + CMPLX(0.0d0, 0.0d0)
13      ELSE
14         sr1(3) = sr1(3) + ex(i,1)*oc(2)-ex(i,2)*oc(1)
15      END IF
16write(unit=*,fmt=*)'T_x_interpolating',sr1(1)
17write(unit=*,fmt=*)'T_y_interpolating',sr1(2)
18write(unit=*,fmt=*)'T_z_interpolating',sr1(3)

但输出的结果还是出现了NaN类型的数据。

麻烦你给些建议要如何修改代码吧。多谢啦。

Kieran · 发表于 2020-9-26 14:11:46

necrohan 发表于 2020-9-26 07:46
write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)
这里用的变量没有初始化 ...

谢谢你的回复。

我仅仅只是粘贴了出现问题部分的代码，ex变量在前面的代码部分已经赋值过了，他们也是复数型的数组。

我有尝试在含有注释的代码（write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)!就是在这里输出的双精度复数型数据，有些数据显示为NaN）前面加入判断数据溢出的IF语句，如下面所示。

[Fortran] 纯文本查看 复制代码

01     IF (ABS(ex(i,2)*oc(3)-ex(i,3)*oc(2)) <= 1.17549435d-38) THEN
02         sr1(1) = sr1(1) + CMPLX(0.0d0, 0.0d0)
03      ELSE
04         sr1(1) = sr1(1) + ex(i,2)*oc(3)-ex(i,3)*oc(2)
05      END IF
06      IF (ABS(ex(i,3)*oc(1)-ex(i,1)*oc(3)) <= 1.17549435d-38) THEN
07         sr1(2) = sr1(2) + CMPLX(0.0d0, 0.0d0)
08      ELSE
09         sr1(2) = sr1(2) + ex(i,3)*oc(1)-ex(i,1)*oc(3)
10      END IF
11      IF (ABS(ex(i,1)*oc(2)-ex(i,2)*oc(1)) <= 1.17549435d-38) THEN
12         sr1(3) = sr1(3) + CMPLX(0.0d0, 0.0d0)
13      ELSE
14         sr1(3) = sr1(3) + ex(i,1)*oc(2)-ex(i,2)*oc(1)
15      END IF
16write(unit=*,fmt=*)'T_x_interpolating',sr1(1)
17write(unit=*,fmt=*)'T_y_interpolating',sr1(2)
18write(unit=*,fmt=*)'T_z_interpolating',sr1(3)

但输出的结果还是出现了NaN类型的数据。

麻烦你给些建议要如何修改代码吧。多谢啦。

Kieran · 发表于 2020-9-26 14:13:33

necrohan 发表于 2020-9-26 07:46
write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)
这里用的变量没有初始化 ...

谢谢你的回复。

我仅仅只是粘贴了出现问题部分的代码，ex变量在前面的代码部分已经赋值过了，它们也是复数型的数组。我有尝试在含有注释的代码（write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)!就是在这里输出的双精度复数型数据，有些数据显示为NaN）前面加入判断数据溢出的IF语句，如下面所示。

[Fortran] 纯文本查看 复制代码

01     IF (ABS(ex(i,2)*oc(3)-ex(i,3)*oc(2)) <= 1.17549435d-38) THEN
02         sr1(1) = sr1(1) + CMPLX(0.0d0, 0.0d0)
03      ELSE
04         sr1(1) = sr1(1) + ex(i,2)*oc(3)-ex(i,3)*oc(2)
05      END IF
06      IF (ABS(ex(i,3)*oc(1)-ex(i,1)*oc(3)) <= 1.17549435d-38) THEN
07         sr1(2) = sr1(2) + CMPLX(0.0d0, 0.0d0)
08      ELSE
09         sr1(2) = sr1(2) + ex(i,3)*oc(1)-ex(i,1)*oc(3)
10      END IF
11      IF (ABS(ex(i,1)*oc(2)-ex(i,2)*oc(1)) <= 1.17549435d-38) THEN
12         sr1(3) = sr1(3) + CMPLX(0.0d0, 0.0d0)
13      ELSE
14         sr1(3) = sr1(3) + ex(i,1)*oc(2)-ex(i,2)*oc(1)
15      END IF
16write(unit=*,fmt=*)'T_x_interpolating',sr1(1)
17write(unit=*,fmt=*)'T_y_interpolating',sr1(2)
18write(unit=*,fmt=*)'T_z_interpolating',sr1(3)

但输出的结果还是出现了NaN类型的数据。麻烦你给些建议要如何修改代码吧。多谢啦。

风平老涡 · 发表于 2020-9-26 22:31:04

Kieran 发表于 2020-9-26 14:13
谢谢你的回复。

我仅仅只是粘贴了出现问题部分的代码，ex变量在前面的代码部分已经赋值过了，它们也是复 ...

先边译成非并行程序，是否有同样问题。

Kieran · 发表于 2020-9-27 09:23:30

风平老涡发表于 2020-9-26 22:31
先边译成非并行程序，是否有同样问题。

好的，我先试下。谢谢你的建议。

wx2203526 · 发表于 2022-12-6 09:37:44

数据溢出怎么解决的呀

[数值问题] 数据溢出

规矩勋章