计算化学公社

标题: Fortran书中关于复数的坑 [打印本页]

作者
Author: beefly 时间: 2021-3-6 17:03
标题: Fortran书中关于复数的坑
本帖最后由 beefly 于 2021-3-6 17:11 编辑

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

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

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

接下来是白海波的《Fortran程序设计权威指南》
(, 下载次数 Times of downloads: 89)
这里讲得很明确，三种形式完全等价。

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

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

为了验证书中的说法是否正确，准备了以下小程序（test.f90），分别在浮点数和复数情况下打印不同精度的pi。为了保证精度，pi是用四精度浮点数计算的。

call test_real()
call test_cplx()
end
subroutine test_real()
real(kind=16) :: a16
real(kind=4) :: a4
real(kind=8) :: a8
real(16) b16
real(8) b8
real(4) b4
real*16 c16
real*8 c8
real*4 c4
! pi
a16=acos(-1.0q0)
b16=a16
c16=a16
a8=a16
b8=a16
c8=a16
a4=a16
b4=a16
c4=a16
write(*,"(f38.34)")a16,b16,c16,a8,b8,c8,a4,b4,c4
return
end
subroutine test_cplx()
complex(kind=16) :: a16
complex(kind=4) :: a4
complex(kind=8) :: a8
complex(16) b16
complex(8) b8
complex(4) b4
complex*16 c16
complex*8 c8
!complex*4 c4
! pi + i * pi
a16=(1.q0,1.q0)*acos(-1.0q0)
b16=a16
c16=a16
a8=a16
b8=a16
c8=a16
a4=a16
b4=a16
!c4=a16
write(*,"(2f38.34)")a16,b16,c16,a8,b8,c8,a4,b4 !,c4
return
end

复制代码

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

总结：

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

作者
Author: wxhwbh 时间: 2021-3-6 19:10
之前吃过这个亏，后来同组的同学就测试过这个，结论是用real*n，complex(n)，complex*2n才是匹配的

作者
Author: Warm_Cloud 时间: 2021-3-6 22:36
统一用kind最保险

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

欢迎光临计算化学公社 (http://bbs.keinsci.com/)