计算化学公社

 找回密码 Forget password
 注册 Register
Views: 5656|回复 Reply: 3

[Fortran] Fortran书中关于复数的坑

[复制链接 Copy URL]

597

帖子

20

威望

4153

eV
积分
5150

Level 6 (一方通行)

发表于 Post on 2021-3-6 17:03:58 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 beefly 于 2021-3-6 17:11 编辑

在Fortran语言中,复数的声明有以下三种形式:
第一种:complex(kind=n) a
第二种:complex(n) a
第三种:complex*n a
其中n=4或8。

但是在不同的Fortran书中,对这三种形式的解释是不一样的。

首先是彭国伦的《Fortran 95程序设计》
01.jpg
这里讲得不够明确:从FORTRAN 77继承下来的第二、三种形式的精度在Fortran后续版本中是否做了修改?给人的感觉,在n相同的情况下,三种声明复数的形式应当具有相同的精度。

接下来是白海波的《Fortran程序设计权威指南》

白海波fortran

白海波fortran

这里讲得很明确,三种形式完全等价。

此外还有一些书只介绍第一种声明复数的形式。这以Chapman的《Fortran程序设计》为代表。

按照以上两本书中的说法,在实际的复数运算中会遇到各种问题。例如,LAPACK库(严格遵守FORTRAN 77标准)中的复数声明全部为“complex*16 a”这种形式。难道LAPACK库用了四精度复数?当然不可能!

为了验证书中的说法是否正确,准备了以下小程序(test.f90),分别在浮点数和复数情况下打印不同精度的pi。为了保证精度,pi是用四精度浮点数计算的。
  1. call test_real()
  2. call test_cplx()
  3. end

  4. subroutine test_real()
  5. real(kind=16) :: a16
  6. real(kind=4)  :: a4
  7. real(kind=8)  :: a8
  8. real(16) b16
  9. real(8) b8
  10. real(4) b4
  11. real*16 c16
  12. real*8 c8
  13. real*4 c4

  14. ! pi
  15. a16=acos(-1.0q0)
  16. b16=a16
  17. c16=a16
  18. a8=a16
  19. b8=a16
  20. c8=a16
  21. a4=a16
  22. b4=a16
  23. c4=a16

  24. write(*,"(f38.34)")a16,b16,c16,a8,b8,c8,a4,b4,c4

  25. return
  26. end

  27. subroutine test_cplx()
  28. complex(kind=16) :: a16
  29. complex(kind=4)  :: a4
  30. complex(kind=8)  :: a8
  31. complex(16) b16
  32. complex(8) b8
  33. complex(4) b4
  34. complex*16 c16
  35. complex*8 c8
  36. !complex*4 c4

  37. ! pi + i * pi
  38. a16=(1.q0,1.q0)*acos(-1.0q0)
  39. b16=a16
  40. c16=a16
  41. a8=a16
  42. b8=a16
  43. c8=a16
  44. a4=a16
  45. b4=a16
  46. !c4=a16

  47. write(*,"(2f38.34)")a16,b16,c16,a8,b8,c8,a4,b4  !,c4

  48. return
  49. end
复制代码
用gfortran或ifort(不支持pgf90)编译后运行。

总结:
  • 对于浮点数,“real(kind=n) a”,“real(n) a”,“real*n a”三种声明的形式等价。
  • 对于复数,第一种和第二种声明形式等价。
  • 复数的第三种声明形式的精度低了一个级别:n=16是双精度,n=8是单精度,至于n=4则根本不支持。书的作者显然是受浮点数声明的影响,写错了。
  • 通过前两种方式,当n=16时可以定义书中没有讲的四精度复数,不过计算非常慢。


评分 Rate

参与人数
Participants 4
eV +18 收起 理由
Reason
伞阳 + 3 谢谢
朙天儿 + 5 正解
sobereva + 8
wxhwbh + 2

查看全部评分 View all ratings

249

帖子

4

威望

3692

eV
积分
4021

Level 6 (一方通行)

发表于 Post on 2021-3-6 19:10:03 | 显示全部楼层 Show all
之前吃过这个亏,后来同组的同学就测试过这个,结论是用real*n,complex(n),complex*2n才是匹配的

评分 Rate

参与人数
Participants 1
eV +5 收起 理由
Reason
北大-陶豫 + 5

查看全部评分 View all ratings

贫困U 退学与疯子工程学院

268

帖子

3

威望

4335

eV
积分
4663

Level 6 (一方通行)

发表于 Post on 2021-3-6 22:36:52 | 显示全部楼层 Show all
统一用kind最保险

评分 Rate

参与人数
Participants 3
eV +9 收起 理由
Reason
ene + 5
卡开发发 + 2 大师大法好
zjxitcc + 2 我很赞同

查看全部评分 View all ratings

142

帖子

0

威望

612

eV
积分
754

Level 4 (黑子)

发表于 Post on 2022-6-17 09:35:17 | 显示全部楼层 Show all
关于这个complex(kind=8),看有些代码里也用double complex表示,同时也有用double precision表示real(kind=8)的。作为C语言习惯用户,感觉还是double这个词看着比较亲切。另外这个double好像只是在主流PC平台上默认是kind=8,而且也可以通过编译器选项比如dbl, -r8等指定。

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 03:26 , Processed in 0.208753 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list