我在Fortran中计算一些数据的fft时遇到了问题。我不知道是否算法有问题，四舍五入，缺乏精度或什么。
代码
module fft_mod
public :: fft1D
integer,parameter :: cp = selected_real_kind(14)
real(cp),parameter :: PI = real(3.14159265358979,cp)
contains
subroutine fft1D(x)
complex(cp), dimension(:), intent(inout) :: x
complex(cp), dimension(:), allocatable :: temp
complex(cp) :: S
integer :: N
complex(cp) :: j ! sqrt(-1)
integer :: i,k
N=size(x)
allocate(temp(N))
j = cmplx(0.0_cp,1.0_cp,cp)
S = cmplx(0.0_cp,0.0_cp,cp)
temp = cmplx(0.0_cp,0.0_cp,cp)
do i = 1,N
do k = 1,N
S = S + x(k)*exp(real(-2.0,cp)*PI*j*real(k-1,cp)*real(i-1,cp)/real(N,cp))
enddo
temp(i) = S
S = cmplx(0.0_cp,0.0_cp,cp)
enddo
x = temp
deallocate(temp)
end subroutine
end module
program test
use fft_mod
implicit none
complex(cp), dimension(10) :: data = (/1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 4.0, 3.0, 2.0, 1.0/)
integer :: i
call fft1D(data)
do i=1,10
write(*,*) data(i)
end do
end program test
在fortran中运行
C:UsersCharlieDesktopfft>gmake
gfortran -J".\mod" -fimplicit-none -Wuninitialized -g -Wall -Wextra -fbacktrace
-fcheck=all -O0 -fopenmp -D_QUAD_PRECISION_ -cpp -c -o .\objtestFFT.o testFFT
.f90
gfortran -o .\test -J".\mod" -fimplicit-none -Wuninitialized -g -Wall -Wextra
-fbacktrace -fcheck=all -O0 -fopenmp -D_QUAD_PRECISION_ -cpp .\objtestFFT.o
C:UsersCharlieDesktopfft>test.exe
( 30.000000000000000 , 0.0000000000000000 )
( -9.4721355260035178 , -3.0776825738331275 )
( 1.2032715918097736E-007, 8.7422769579070803E-008)
(-0.52786408204828272 ,-0.72654221835813126 )
( 5.6810824045072650E-008, 1.7484555003832725E-007)
( 1.0325074129013956E-013, 2.6226834001115759E-007)
( -8.5216018574918451E-008, 2.6226836247200680E-007)
(-0.52786395855490920 , 0.72654325051559143 )
( -4.8130813040669906E-007, 3.4969128892559098E-007)
( -9.4721398159606647 , 3.0776922072585111 )
但是在matlab中运行相同的数据集得到
format long ; g = [1:5 5:-1:1]; fft(g)'
ans =
30.000000000000000
-9.472135954999580 + 3.077683537175253i
0
-0.527864045000420 + 0.726542528005361i
0
0
0
-0.527864045000420 - 0.726542528005361i
0
-9.472135954999580 - 3.077683537175253i
我相信我使用双精度使用selected_real_kind(14)，但看起来结果充其量只是单精度。我敢肯定，一些洒实(，cp)'的是不必要的，但我不知道在哪里，如何或为什么它看起来像结果是单一精度与matlab相比。
任何帮助都非常感谢!
更新:使用被接受的答案的建议，这里唯一需要更改的是:
real(cp),parameter :: PI = real(3.14159265358979,cp)
real(cp),parameter :: PI = 3.14159265358979_cp
问题在于如何定义实数，特别是pi。当定义
real(cp),parameter :: PI = real(3.14159265358979,cp)
将参数3.14159265358979传递给函数real。但是实数有默认的单精度，所以你的实数在进入函数时被转换成单精度。考虑下面的例子:
program main
integer,parameter :: cp = selected_real_kind(14)
real(cp),parameter :: pi = real(3.14159265358979,cp)
real(cp),parameter :: pj = 3.14159265358979_cp
write(*,*) pi
write(*,*) pj
end program main
用pgfortran编译，没有选项，这给了我:
3.141592741012573
3.141592653589790
当定义任何实数时，应使用[]_cp而不是real([],cp)来分配kind。
编辑:这个问题也会影响您如何定义0.0, 1.0和2.0，但是这些数字可以精确地转换为二进制，并且不会遭受相同的舍入误差。
如果你只需要一些"接近"的东西，上面这个公认的答案是一个合理的解决方案。双精度(Real(8))，因为您已经显式地将Pi定义为接近Real(8)的位数。Pi的实值(8)为
3.1415926535897931出口。
的参数 3.14159265358979如果您希望具有与_cp"一致的更一般的"全精度"。因此，您可能希望使用诸如
之类的内容Real(cP), Parameter :: Pi = 4_cp*ATan(1_cp)
***编辑:感谢francescalus指出拼写错误，这应该是ATan(1._cp)，但实际上，正如下面所解释的，参数'ized ';应该使用数字，避免使用"_cp"
现在，如果_cp是16，你的Pi将自动为:
3.14159265358979323846264338327950280现在，无论何时更改_cp，您的代码将具有与_cp"一致的"全精度"。自动.
顺便说一句，你可能还希望"参数"某些基本数字，例如:
Real(cp), Parameter :: Zero = 0_cp
Real(cp), Parameter :: One = 1_cp
:
…等
***编辑:感谢francescalus发现拼写错误，但在这些特定的表达式中，虽然0.0_cp和1.0_cp可能更好，但没关系，因为它们的Params声明会照顾到>
然后，在你的代码中，你可以这样写，例如:...
Real(cp), Parameter :: Pi = Four*ATan(One)
...
S = Cmplx(Zero, Zero)*Exp(-Two) ! etc etc
并且不必担心到处添加_cp等，并且更容易阅读/调试等。
这个策略还有其他优点，特别是在遗留代码中，比如
If( a == 1.0 ) Then
与
If( a == One ) Then
…但那是另一回事了。
顺便说一句，在某种程度上也是一个"风格"问题，在Fortran中，算术默认为"最大"。表达式中的精度/类型，因此
S = S + x(k)*exp(real(-2.0,cp)*PI*j*real(k-1,cp)*real(i-1,cp)/real(N,cp))
应该等同于
S = S + x(k)*exp(-2*PI*j*(k-1)*(i-1)/N)
或者更好的
S = S + x(k)*exp(-Two*PI*j*(k-1)*(i-1)/N)
就我个人而言，我发现越容易阅读，就越容易把它写对。
…尽管如上所述，如果你的算术只包含整数，那么可能需要额外的考虑。但是，如果您"参数"` ` ize ` `;只有Real，那么你就不应该有使用"2"的风险(比如div/0)。在真正需要2.0版本的地方，如果你使用"Two"
最后，当你有"Parameter'ized"Pi，为什么不是2*Pi，例如:
...
Real(cp), Parameter :: TwoPi = Two*Pi
现在，表达式
S = S + x(k)*exp(-Two*PI*j*(k-1)*(i-1)/N)
可以
S = S + x(k)*exp(-TwoPI*j*(k-1)*(i-1)/N)
可以节省一点CPU。事实上，扩展这个逻辑可以进一步提高代码的性能。
…另一个完全不同的问题是，您是否知道var Temp()可以是"自动"的，例如
complex(cp) :: temp(Size(x, Dim = 1))
则不需要Allocatable等。尽管从我的角度来看，看起来完全没有"S"也可以编写代码。和"Temp(:)"，提供更简单和更快的代码。
*****附录:根据对我的回答的一些评论，我认为它可能有助于显示对OP代码的可能更改以提高性能，并在某种程度上进一步提高精度。
然而，在此之前，OP没有说明为什么需要特定程度的精度。这可能是许多工业/现实世界设置中的一个重要问题。例如，在某些情况下，速度比精度重要得多，只要精度足够好，可以提供特定于上下文的可靠性/决策。例如，可以想象，在某些情况下，计算应该都是Real(4)。这件事需要单独讨论。
尽管如此，Ross的答案在"内存表示"方面纠正了精度问题，而我的原始答案在"你实际上需要提供足够的符号来开始"方面提高了精度。(即使声明是正确的)，减少FLOPS的数量不仅可以提高速度，还可以帮助提高精度，因为每个FLOP都会引入截断/舍入的可能性，特别是对于单个表达式中的长FLOPS序列。
此外，这个OP问题可以很好地解决与"标识"配方(4 *(1)每股)。在许多情况下，这要么是不可能的，要么是不实际的，然后需要一些不那么优雅的方法，但我把这些留到另一天。
以下只是一个/两个可能的选择。
我还没有测试过。
所以它可能需要调整;请随时提出改进意见。
如果运气好的话，它可能会减少大约80%的FLOPs。
…如果方便的话，我将感谢OP对不同的实现进行基准测试并报告结果。您可能还希望对_cp = 4,8,16进行基准测试，以演示速度和精度之间的权衡。
这个替代版本显然需要对调用s/r进行更新。
module fft_mod
!
public :: fft1D
!
Integer,parameter :: cp = selected_real_kind(14) ! just out of curiosity, why "14", see my comments on "reality vs. precision"
!>>real(cp),parameter :: PI = real(3.14159265358979,cp)
!
Real(cp), Parameter :: Zero = 0_cp ! these types of declarations can be created in a separate Module,
! say "MyDeclarationModule", and then just access with "Use"
!
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! NOTE: With respect to francescalus comments on issues with "1_cp" syntax, the usage here works fine as shown in the previous result,
! though francescalus' comments are fair in suggesting that other approaches make for better coding.
! However, and to be clear, I don't actually use "_xx" style myslef. I have used it here ONLY to be consistent with the
! the OP's and earlier answers. The usage of it here is NOT a recommendation. A discussion of the alternatives
! is much too long and strays from the immediate issues too far.
! ... Perhaps francescalus will take up the mantle and re-write the OP's code for alternatives
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!
Real(cp), Parameter :: One = 1_cp
Real(cp), Parameter :: Two = 2_cp
Real(cp), Parameter :: Four = 4_cp
!
! ... after this point, there is no need to see "_cp" again, in this example
!
Real(cp), Parameter :: Pi = Four*ATan(One) ! this guarrantees maximal precision for Pi, up to "_cp"
!
Real(cp), Parameter :: TwoPi = Two*Pi ! Vladimir: now, and only now (that I've talem most everything out of the loop),
! this declaration has less value ... unlike previously, when it was
! still in the inner NxN, and when it saved approx 15 - 20% of FLOPs
! Crucially, if you do a lot of computational mathematics, TwoPi, PiBy2, Root2Pi, etc etc
! arise with considerable frequency, and creating these as Parameters makes for much
! improvement, and I would think it a practice to be encouraged.
!
Complex(cp), Parameter :: SqrtM1 = Cmplx(Zero, One) ! sqrt(-1), this is your "j",
! sorry but "j" sounds just too much like an "integer" to me
! (for us "old timers"), so I changed the name to something more meaningful to me
!
!
Contains
!
Pure Subroutine fft1D(x, n, Temp) ! OPTIONAL, return the results explicitly to save vector assignment, and preserve original data
! OPTIONAL, send the size of x() explicitly as n, much cleaner
! you could leave it as s/r fft1D(x), and most of the savings below would still apply
! ... also, my practice is to make everything Pure/Elemental where reasonable
!
Integer, Intent(In) :: n ! Optional, but cleaner
! complex(cp), Intent(InOut) :: x(n) ! consider as just Intent(In), and return Temp() as Intent(Out)
complex(cp), Intent(In) :: x(n) ! consider as just Intent(In), and return Temp() as Intent(Out)
!
! Locals
! ------
!
! Complex(cp), Allocatable :: temp(:) ! no need for Allocatable
! Complex(cp) :: temp(Size(x), Dim = 1) ! I prefer to pass n explicitly
Complex(cp), Intent(Out) :: temp(n) ! make Intent(Out) to preserve data and save a vector assignment
!
! complex(cp) :: S
! integer :: N
! complex(cp) :: j ! sqrt(-1)
!
Integer :: i, k
!
Complex(cp) :: TPSdivN ! new reduce 4 nxn mults to 1
!
Complex(cp) :: TPSdivNiM1 ! new, to reduce 5 nxn mults to 1 nx1
!
Real(cp) :: rKM1(n) ! new, to obviate nxn k-1's in inner loop
!
! N=size(x) ! my preference is to pass N explicitly,
! but can be done this way if prefered (need to make appropriate changes above)
!
! allocate(temp(N)) ! Temp() can be either automatic or an Arg, no need for Allocate
! j = cmplx(0.0_cp,1.0_cp,cp) ! make Parameter, and rename to SqrtM1 (personal pref)
! S = cmplx(0.0_cp,0.0_cp,cp) ! Not required with "improved" inner loop
!
!
temp(1:n) = Cmplx(Zero, Zero) ! S is not needed, just start with Temp() directly
!
TPSdivN = -TwoPi*SqrtM1/n ! new, move everything out all loops that can be
!
ForAll( k=1:n) rKM1(k) = k - 1 ! new, this allows the elimination of NxN "k-1's" in the inner
! my preference is for ForAll's, but you could use Do,
!
do i = 1,N
!
TPSdivNiM1 = TPSdivN*(i-1) ! new. everything out of inner loop that can be
!
! improved, inner do, but could use "no-Do" alternative as illustrated below
!
! do k = 1,N
!>> S = S + x(k)*exp(real(-2.0,cp)*PI*j*real(k-1,cp)*real(i-1,cp)/real(N,cp)) ! original, many unneccessary nxn FLOPs, type conversions, etc
! temp(i) = temp(i) + x(k)*Exp(TPSdivNiM1*rKM1(k))) ! use array of k-1's, then the inner k-loop is no longer required,
! and can use Product()/Sum() or Dot_Product() intrinsics, see example below
! End Do
!
! alternative "direct array" approach to inner loop, or "no-DO" version ... there are other possibilities.
!
Temp(i) = Dot_Product( x(1:n), Exp(TPSdivNiM1*rKM1(1:n)) ) ! this approach can have certain advantages depending on circumstances
!
! temp(i) = S ! not needed
! S = cmplx(0.0_cp,0.0_cp,cp) ! not needed
!
End Do
!
! x(1:n) = temp(1:n) ! not need if Temp() is Intent(Out) that save this vector assignment, and the original data
! deallocate(temp) ! not needed
!
End Subroutine fft1D
!
! end module
End Module fft_mod