Fortran Coder

查看: 224|回复: 6

[数值问题] 数据溢出

[复制链接]

32

帖子

10

主题

0

精华

熟手

F 币
166 元
贡献
102 点
发表于 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)

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

使用道具 举报

138

帖子

2

主题

0

精华

大师

F 币
910 元
贡献
475 点

规矩勋章

发表于 2020-9-26 07:46:12 | 显示全部楼层
write(unit=*,fmt=*)'T_x_interpolating',ex(i,2) * oc(3) - ex(i,3) * oc(2)
这里用的变量没有初始化

32

帖子

10

主题

0

精华

熟手

F 币
166 元
贡献
102 点
 楼主| 发表于 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类型的数据。

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

32

帖子

10

主题

0

精华

熟手

F 币
166 元
贡献
102 点
 楼主| 发表于 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类型的数据。

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

32

帖子

10

主题

0

精华

熟手

F 币
166 元
贡献
102 点
 楼主| 发表于 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类型的数据。麻烦你给些建议要如何修改代码吧。多谢啦。

66

帖子

2

主题

0

精华

专家

F 币
476 元
贡献
217 点

规矩勋章

发表于 2020-9-26 22:31:04 | 显示全部楼层
Kieran 发表于 2020-9-26 14:13
谢谢你的回复。

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

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

32

帖子

10

主题

0

精华

熟手

F 币
166 元
贡献
102 点
 楼主| 发表于 2020-9-27 09:23:30 | 显示全部楼层
风平老涡 发表于 2020-9-26 22:31
先边译成非并行程序,是否有同样问题。

好的,我先试下。谢谢你的建议。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2020-10-28 01:49

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表