数据溢出

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

大家好，

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

我的出问题部分的代码。

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

!External variables
INTEGER                       :: vd
DOUBLE PRECISION, ALLOCATABLE :: kp(:,:)
INTEGER                       :: nu_wa
INTEGER                       :: nu_ws
INTEGER, ALLOCATABLE          :: nd(:)
DOUBLE PRECISION, ALLOCATABLE :: hr(:,:)
DOUBLE COMPLEX, ALLOCATABLE   :: tb(:,:)
DOUBLE COMPLEX, ALLOCATABLE   :: ec(:,:)
DOUBLE PRECISION, ALLOCATABLE :: ev(:,:)
DOUBLE COMPLEX, ALLOCATABLE   :: vh(:,:)
DOUBLE PRECISION              :: fermi
DOUBLE PRECISION              :: wf
DOUBLE PRECISION              :: bv
INTEGER                       :: nv
DOUBLE PRECISION              :: mv
INTEGER                       :: nn
DOUBLE PRECISION              :: dk(2)
DOUBLE COMPLEX, ALLOCATABLE   :: sr(:,:)
DOUBLE COMPLEX                :: sr1(3)
INTEGER                       :: nu_kp

INTEGER                       :: i, j, k, l, m
DOUBLE COMPLEX                :: dc(3)
DOUBLE COMPLEX, ALLOCATABLE   :: ex(:,:)
DOUBLE PRECISION              :: de

DOUBLE PRECISION              :: ey
DOUBLE COMPLEX                :: gr
DOUBLE COMPLEX                :: ga
DOUBLE COMPLEX                :: gd(2)
DOUBLE COMPLEX                :: su
DOUBLE COMPLEX, ALLOCATABLE   :: mt1(:,:)
DOUBLE COMPLEX, ALLOCATABLE   :: mt2(:,:)
DOUBLE COMPLEX, ALLOCATABLE   :: mt3(:,:)
DOUBLE COMPLEX, ALLOCATABLE   :: mt4(:,:)
DOUBLE COMPLEX, ALLOCATABLE   :: mt5(:,:)
DOUBLE COMPLEX                :: mt6(8)
DOUBLE COMPLEX                :: sp(6)
DOUBLE PRECISION, ALLOCATABLE :: nm(:)
DOUBLE COMPLEX                :: oc(3)

INTEGER                       :: omp_get_num_threads
INTEGER                       :: num_thread
INTEGER                       :: omp_get_thread_num
INTEGER                       :: id_thread
DOUBLE COMPLEX, ALLOCATABLE   :: sr2(:,:,:)
!$OMP PARALLEL DO PRIVATE(i, j, k, l, ey, sr1, mt1, mt2, mt3, mt4, mt5, mt6, nm, oc, gr, ga, gd, su)&

!$OMP             SHARED(ev, ec, ex, nv, nn, mv, wf, de, nu_kp, nu_wa, vd, vh, bv, dk, sr2)

ALLOCATE (sr2(num_thread,nv+nn,3))
sr2 = CMPLX(0.0d0, 0.0d0)

DO i = 1, nv + nn - 2, 1
   id_thread = omp_get_thread_num() + 1

   ey = mv - wf + DBLE(i) * de
   sr1 = CMPLX(0.0d0, 0.0d0)

   DO j = 1, nu_kp, 1
      !Normalized factor part
      DO k = 1, nu_wa, 1
         mt1(1,:) = DCONJG(ec(1+(j-1)*nu_wa:j*nu_wa,k))
         mt2(:,1) = ec(1+(j-1)*nu_wa:j*nu_wa,k)
         nm(k) = REAL(SUM(MATMUL(mt1,mt2)))
         nm(k) = 1.0d0 / DSQRT(nm(k))
      END DO
      IF (vd .EQ. 0) THEN
          mt5 = vh(1+(j-1)*nu_wa:j*nu_wa,1:nu_wa)
      ELSE
          mt5 = vh(1+(j-1)*nu_wa:j*nu_wa,1+nu_wa:2*nu_wa)
      END IF

      oc = CMPLX(0.0d0, 0.0d0)
      DO k = 1, nu_wa, 1
         mt1(1,:) = DCONJG(ec(1+(j-1)*nu_wa:j*nu_wa,k))
         mt2(:,1) = ec(1+(j-1)*nu_wa:j*nu_wa,k)
         gr = CMPLX (1.0d0, 0.0d0) / CMPLX(ey - ev(j,k), bv)
         ga = CMPLX (1.0d0, 0.0d0) / CMPLX(ey - ev(j,k), 0.0d0 - bv)
         gd(2) = (gr * gr + ga * ga)
         DO l = 1, nu_wa, 1
            mt3(1,:) = DCONJG(ec(1+(j-1)*nu_wa:j*nu_wa,l))
            mt4(:,1) = ec(1+(j-1)*nu_wa:j*nu_wa,l)
            gr = CMPLX (1.0d0, 0.0d0) / CMPLX(ey - ev(j,l), bv)
            ga = CMPLX (1.0d0, 0.0d0) / CMPLX(ey - ev(j,l), 0.0d0 - bv)
            gd(1) = (gr - ga)
            sp = CMPLX(0.0d0, 0.0d0)
            DO m = 1, nu_wa / 2, 1
               mt6(1) = mt2(m,1)
               mt6(2) = mt2(m+nu_wa/2,1)
               mt6(3) = mt3(1,m)
               mt6(4) = mt3(1,m+nu_wa/2)
               mt6(5) = mt1(1,m)
               mt6(6) = mt1(1,m+nu_wa/2)
               mt6(7) = mt4(m,1)
               mt6(8) = mt4(m+nu_wa/2,1)
               sp(1) = sp(1) + mt6(4) * mt6(1) + mt6(3) * mt6(2)
               sp(2) = sp(2) + dc(1) * (mt6(4) * mt6(1) - mt6(3) * mt6(2))
               sp(3) = sp(3) + mt6(3) * mt6(1) - mt6(4) * mt6(2)
               sp(4) = sp(4) + mt6(6) * mt6(7) + mt6(5) * mt6(8)
               sp(5) = sp(5) + dc(1) * (mt6(6) * mt6(7) - mt6(5) * mt6(8))
               sp(6) = sp(6) + mt6(5) * mt6(7) - mt6(6) * mt6(8)
            END DO
            su = (SUM(MATMUL(mt1,MATMUL(mt5,mt4)))*sp(1)-sp(4)*SUM(MATMUL(mt3,MATMUL(mt5,mt2))))/CMPLX(2.0d0, 0.0d0)
            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)
            su = (SUM(MATMUL(mt1,MATMUL(mt5,mt4)))*sp(2)-sp(5)*SUM(MATMUL(mt3,MATMUL(mt5,mt2))))/CMPLX(2.0d0, 0.0d0)
            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)
            su = (SUM(MATMUL(mt1,MATMUL(mt5,mt4)))*sp(3)-sp(6)*SUM(MATMUL(mt3,MATMUL(mt5,mt2))))/CMPLX(2.0d0, 0.0d0)
            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)
         END DO
      END DO
      sr1(1) = sr1(1) + ex(i,2) * oc(3) - ex(i,3) * oc(2)
      sr1(2) = sr1(2) + ex(i,3) * oc(1) - ex(i,1) * oc(3)
      sr1(3) = sr1(3) + ex(i,1) * oc(2) - ex(i,2) * oc(1)
write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)!就是在这里输出的双精度复数型数据，有些数据显示为NaN
write(unit=*,fmt=*)'T_y_interpolating',ex(i,3) * oc(1) - ex(i,1) * oc(3)!就是在这里输出的双精度复数型数据，有些数据显示为NaN
write(unit=*,fmt=*)'T_z_interpolating',ex(i,1) * oc(2) - ex(i,2) * oc(1)!就是在这里输出的双精度复数型数据，有些数据显示为NaN
   END DO

   sr2(id_thread,1+i,1) = sr1(1)
   sr2(id_thread,1+i,2) = sr1(2)
   sr2(id_thread,1+i,3) = sr1(3)
END DO
!$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] 纯文本查看 复制代码

      IF (ABS(ex(i,2)*oc(3)-ex(i,3)*oc(2)) <= 1.17549435d-38) THEN
         sr1(1) = sr1(1) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(1) = sr1(1) + ex(i,2)*oc(3)-ex(i,3)*oc(2)
      END IF
      IF (ABS(ex(i,3)*oc(1)-ex(i,1)*oc(3)) <= 1.17549435d-38) THEN
         sr1(2) = sr1(2) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(2) = sr1(2) + ex(i,3)*oc(1)-ex(i,1)*oc(3)
      END IF
      IF (ABS(ex(i,1)*oc(2)-ex(i,2)*oc(1)) <= 1.17549435d-38) THEN
         sr1(3) = sr1(3) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(3) = sr1(3) + ex(i,1)*oc(2)-ex(i,2)*oc(1)
      END IF
write(unit=*,fmt=*)'T_x_interpolating',sr1(1)
write(unit=*,fmt=*)'T_y_interpolating',sr1(2)
write(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] 纯文本查看 复制代码

     IF (ABS(ex(i,2)*oc(3)-ex(i,3)*oc(2)) <= 1.17549435d-38) THEN
         sr1(1) = sr1(1) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(1) = sr1(1) + ex(i,2)*oc(3)-ex(i,3)*oc(2)
      END IF
      IF (ABS(ex(i,3)*oc(1)-ex(i,1)*oc(3)) <= 1.17549435d-38) THEN
         sr1(2) = sr1(2) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(2) = sr1(2) + ex(i,3)*oc(1)-ex(i,1)*oc(3)
      END IF
      IF (ABS(ex(i,1)*oc(2)-ex(i,2)*oc(1)) <= 1.17549435d-38) THEN
         sr1(3) = sr1(3) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(3) = sr1(3) + ex(i,1)*oc(2)-ex(i,2)*oc(1)
      END IF
write(unit=*,fmt=*)'T_x_interpolating',sr1(1)
write(unit=*,fmt=*)'T_y_interpolating',sr1(2)
write(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] 纯文本查看 复制代码

     IF (ABS(ex(i,2)*oc(3)-ex(i,3)*oc(2)) <= 1.17549435d-38) THEN
         sr1(1) = sr1(1) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(1) = sr1(1) + ex(i,2)*oc(3)-ex(i,3)*oc(2)
      END IF
      IF (ABS(ex(i,3)*oc(1)-ex(i,1)*oc(3)) <= 1.17549435d-38) THEN
         sr1(2) = sr1(2) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(2) = sr1(2) + ex(i,3)*oc(1)-ex(i,1)*oc(3)
      END IF
      IF (ABS(ex(i,1)*oc(2)-ex(i,2)*oc(1)) <= 1.17549435d-38) THEN
         sr1(3) = sr1(3) + CMPLX(0.0d0, 0.0d0)
      ELSE
         sr1(3) = sr1(3) + ex(i,1)*oc(2)-ex(i,2)*oc(1)
      END IF
write(unit=*,fmt=*)'T_x_interpolating',sr1(1)
write(unit=*,fmt=*)'T_y_interpolating',sr1(2)
write(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

数据溢出怎么解决的呀

[数值问题] 数据溢出

规矩勋章