Fortran书中关于复数的坑

beefly

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

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

首先是彭国伦的《Fortran 95程序设计》

这里讲得不够明确：从FORTRAN 77继承下来的第二、三种形式的精度在Fortran后续版本中是否做了修改？给人的感觉，在n相同的情况下，三种声明复数的形式应当具有相同的精度。

接下来是白海波的《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. 
   
5. subroutine test_real()
   
6. real(kind=16) :: a16
   
7. real(kind=4)  :: a4
   
8. real(kind=8)  :: a8
   
9. real(16) b16
   
10. real(8) b8
    
11. real(4) b4
    
12. real*16 c16
    
13. real*8 c8
    
14. real*4 c4
    
15. 
    
16. ! pi
    
17. a16=acos(-1.0q0)
    
18. b16=a16
    
19. c16=a16
    
20. a8=a16
    
21. b8=a16
    
22. c8=a16
    
23. a4=a16
    
24. b4=a16
    
25. c4=a16
    
26. 
    
27. write(*,"(f38.34)")a16,b16,c16,a8,b8,c8,a4,b4,c4
    
28. 
    
29. return
    
30. end
    
31. 
    
32. subroutine test_cplx()
    
33. complex(kind=16) :: a16
    
34. complex(kind=4)  :: a4
    
35. complex(kind=8)  :: a8
    
36. complex(16) b16
    
37. complex(8) b8
    
38. complex(4) b4
    
39. complex*16 c16
    
40. complex*8 c8
    
41. !complex*4 c4
    
42. 
    
43. ! pi + i * pi
    
44. a16=(1.q0,1.q0)*acos(-1.0q0)
    
45. b16=a16
    
46. c16=a16
    
47. a8=a16
    
48. b8=a16
    
49. c8=a16
    
50. a4=a16
    
51. b4=a16
    
52. !c4=a16
    
53. 
    
54. write(*,"(2f38.34)")a16,b16,c16,a8,b8,c8,a4,b4  !,c4
    
55. 
    
56. return
    
57. end
    

复制代码

用gfortran或ifort（不支持pgf90）编译后运行。

总结：

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

wxhwbh

之前吃过这个亏，后来同组的同学就测试过这个，结论是用real*n，complex(n)，complex*2n才是匹配的

Warm_Cloud

统一用kind最保险

snljty2

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