在Fortran中使用r2c和c2r FFTW時，前向和後向維度是否相同？

-1

吹是一個主文件在Fortran中使用r2c和c2r FFTW時，前向和後向維度是否相同？

PROGRAM SPHEROID 

USE nrtype 
USE SUB_INFO 
INCLUDE "/usr/local/include/fftw3.f" 

INTEGER(I8B) :: plan_forward, plan_backward 
INTEGER(I4B) :: i, t, int_N 

REAL(DP) :: cth_i, sth_i, real_i, perturbation 
REAL(DP) :: PolarEffect, dummy, x1, x2, x3 

REAL(DP), DIMENSION(4096) :: dummy1, dummy2, gam, th, ph 
REAL(DP), DIMENSION(4096) :: k1, k2, k3, k4, l1, l2, l3, l4, f_in 

COMPLEX(DPC), DIMENSION(2049) :: output1, output2, f_out 

CHARACTER(1024) :: baseOutputFilename 
CHARACTER(1024) :: outputFile, format_string 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

int_N = 4096 

! File Open Section 

format_string = '(I5.5)' 

! Write the coodinates at t = 0 

do i = 1, N 
    real_i = real(i) 
    gam(i) = 2d0*pi/real_N 
    perturbation = 0.01d0*dsin(2d0*pi*real_i/real_N) 
    ph(i) = 2d0*pi*real_i/real_N + perturbation 
    th(i) = pi/3d0 + perturbation 
end do 

! Initialization Section for FFTW PLANS 

call dfftw_plan_dft_r2c_1d(plan_forward, int_N, f_in, f_out, FFTW_ESTIMATE) 
call dfftw_plan_dft_c2r_1d(plan_backward, int_N, f_out, f_in, FFTW_ESTIMATE) 

! Runge-Kutta 4th Order Method Section 

do t = 1, Iter_N 

    call integration(th, ph, gam, k1, l1) 

    do i = 1, N 
     dummy1(i) = th(i) + 0.5d0*dt*k1(i) 
    end do 
    do i = 1, N 
     dummy2(i) = ph(i) + 0.5d0*dt*l1(i) 
    end do 

    call integration(dummy1, dummy2, gam, k2, l2) 

    do i = 1, N 
     dummy1(i) = th(i) + 0.5d0*dt*k2(i) 
    end do 
    do i = 1, N 
     dummy2(i) = ph(i) + 0.5d0*dt*l2(i) 
    end do 

    call integration(dummy1, dummy2, gam, k3, l3) 

    do i = 1, N 
     dummy1(i) = th(i) + dt*k3(i) 
    end do 
    do i = 1, N 
     dummy2(i) = ph(i) + dt*l3(i) 
    end do 

    call integration(dummy1, dummy2, gam, k4, l4) 

    do i = 1, N 

     cth_i = dcos(th(i)) 
     sth_i = dsin(th(i)) 
     PolarEffect = (nv-sv)*dsqrt(1d0+a*sth_i**2) + (nv+sv)*cth_i 
     PolarEffect = PolarEffect/(sth_i**2) 
     th(i) = th(i) + dt*(k1(i) + 2d0*k2(i) + 2d0*k3(i) + k4(i))/6d0 
     ph(i) = ph(i) + dt*(l1(i) + 2d0*l2(i) + 2d0*l3(i) + l4(i))/6d0 
     ph(i) = ph(i) + dt*0.25d0*PolarEffect/pi 

    end do 

!! Fourier Filtering Section 

    call dfftw_execute_dft_r2c(plan_forward, th, output1) 

    do i = 1, N/2+1 
     dummy = abs(output1(i)) 
     if (dummy.lt.threshhold) then 
      output1(i) = dcmplx(0.0d0) 
     end if 
    end do 

    call dfftw_execute_dft_c2r(plan_backward, output1, th) 

    do i = 1, N 
     th(i) = th(i)/real_N 
    end do 

    call dfftw_execute_dft_r2c(plan_forward, ph, output2) 

    do i = 1, N/2+1 
     dummy = abs(output2(i)) 
     if (dummy.lt.threshhold) then 
      output2(i) = dcmplx(0.0d0) 
     end if 
    end do 

    call dfftw_execute_dft_c2r(plan_backward, output2, ph) 

    do i = 1, N 
     ph(i) = ph(i)/real_N 
    end do 

!! Data Writing Section 

    write(baseOutputFilename, format_string) t 
    outputFile = "xyz" // baseOutputFilename 
    open(unit=7, file=outputFile) 
    outputFile = "Fsptrm" // baseOutputFilename 
    open(unit=8, file=outputFile) 

    do i = 1, N 
     x1 = dsin(th(i))*dcos(ph(i)) 
     x2 = dsin(th(i))*dsin(ph(i)) 
     x3 = dsqrt(1d0+a)*dcos(th(i)) 
     write(7,*) x1, x2, x3 
    end do 

    do i = 1, N/2+1 
     write(8,*) abs(output1(i)), abs(output2(i)) 
    end do 

    close(7) 
    close(8) 

    do i = 1, N/2+1 
     output1(i) = dcmplx(0.0d0) 
    end do 

    do i = 1, N/2+1 
     output2(i) = dcmplx(0.0d0) 
    end do 

end do 

! Destroying Process for FFTW PLANS 

call dfftw_destroy_plan(plan_forward) 
call dfftw_destroy_plan(plan_backward) 

END PROGRAM

下面是集成

! We implemented Shelly's spectrally accurate convergence method 

SUBROUTINE integration(in1,in2,in3,out1,out2) 

USE nrtype 
USE SUB_INFO 

INTEGER(I4B) :: i, j 
REAL(DP) :: th_i, th_j, gi, ph_i, ph_j, gam_j, v1, v2 
REAL(DP), DIMENSION(N), INTENT(INOUT) :: in1, in2, in3, out1, out2 
REAL(DP) :: ui, uj, part1, part2, gj, cph, sph 
REAL(DP) :: denom, numer, temp 

do i = 1, N 
    out1(i) = 0d0 
end do 
do i = 1, N 
    out2(i) = 0d0 
end do 

do i = 1, N 

    th_i = in1(i) 
    ph_i = in2(i) 
    ui = dcos(th_i) 
    part1 = dsqrt(1d0+a)/(dsqrt(-a)*ui+dsqrt(1d0+a-a*ui*ui)) 
    part1 = part1**(dsqrt(-a)) 
    part2 = (dsqrt(1d0+a-a*ui*ui)+ui)/(dsqrt(1d0+a-a*ui*ui)-ui) 
    part2 = dsqrt(part2) 

    gi = dsqrt(1d0-ui*ui)*part1*part2 

    do j = 1, N 

     if (mod(i+j,2).eq.1) then 

      th_j = in1(j) 
      ph_j = in2(j) 
      gam_j = in3(j) 

      uj = dcos(th_j) 
      part1 = dsqrt(1d0+a)/(dsqrt(-a)*uj+dsqrt(1d0+a-a*uj*uj)) 
      part1 = part1**(dsqrt(-a)) 
      part2 = (dsqrt(1d0+a-a*uj*uj)+uj)/(dsqrt(1d0+a-a*uj*uj)-uj) 
      part2 = dsqrt(part2) 

      gj = dsqrt(1d0-ui*ui)*part1*part2 

      cph = dcos(ph_i-ph_j) 
      sph = dsin(ph_i-ph_j) 

      numer = dsqrt(1d0-uj*uj)*sph 
      denom = (gj/gi*(1d0-ui*ui) + gi/gj*(1d0-uj*uj))*0.5d0 
      denom = denom - dsqrt((1d0-ui*ui)*(1d0-uj*uj))*cph 
      denom = denom + krasny_delta 
      v1 = -0.25d0*gam_j*numer/denom/pi 

      temp = dsqrt(1d0+(1d0-ui*ui)*a) 
      numer = -(gj/gi)*(temp+ui) 
      numer = numer + (gi/gj)*((1d0-uj*uj)/(1d0-ui*ui))*(temp-ui) 
      numer = numer + 2d0*ui*dsqrt((1d0-uj*uj)/(1d0-ui*ui))*cph 
      numer = 0.5d0*numer 
      v2 = -0.25d0*gam_j*numer/denom/pi 

      out1(i) = out1(i) + 2d0*v1 
      out2(i) = out2(i) + 2d0*v2 

     end if 

    end do 

end do 

END

下面的子程序文件模塊文件

module nrtype 
Implicit none 
!integer 
INTEGER, PARAMETER :: I8B = SELECTED_INT_KIND(20) 
INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9) 
INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4) 
INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2) 
!real 
INTEGER, PARAMETER :: SP = KIND(1.0) 
INTEGER, PARAMETER :: DP = KIND(1.0D0) 
!complex 
INTEGER, PARAMETER :: SPC = KIND((1.0,1.0)) 
INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0)) 
!defualt logical 
INTEGER, PARAMETER :: LGT = KIND(.true.) 
!mathematical constants 
REAL(DP), PARAMETER :: pi = 3.141592653589793238462643383279502884197_dp 
!derived data type s for sparse matrices,single and double precision 
!User-Defined Constants 
INTEGER(I4B), PARAMETER :: N = 4096, Iter_N = 20000 
REAL(DP), PARAMETER :: real_N = 4096d0 
REAL(DP), PARAMETER :: a = -0.1d0, dt = 0.001d0, krasny_delta = 0.01d0 
REAL(DP), PARAMETER :: nv = 0d0, sv = 0d0, threshhold = 0.00000000001d0 
!N : The Number of Point Vortices, Iter_N * dt = Total time, dt : Time Step 
!krasny_delta : Smoothing Parameter introduced by R.Krasny 
!nv : Northern Vortex Strength, sv : Southern Vortex Strength 
!a : The Eccentricity in the direction of z , threshhold : Filtering Threshhold 
end module nrtype

下面是一個子程序信息文件

MODULE SUB_INFO 

INTERFACE 
    SUBROUTINE integration(in1,in2,in3,out1,out2) 
    USE nrtype 
    INTEGER(I4B) :: i, j 
    REAL(DP) :: th_i, th_j, gi, ph_i, ph_j, gam_j, v1, v2 
    REAL(DP), DIMENSION(N), INTENT(INOUT) :: in1, in2, in3, out1, out2 
    REAL(DP) :: ui, uj, part1, part2, gj, cph, sph 
    REAL(DP) :: denom, numer, temp 
    END SUBROUTINE 
END INTERFACE 

END MODULE

我使用下面的命令編譯他們

gfortran -o p0 -fbounds-check nrtype.f90 spheroid_sub_info.f90 spheroid_sub_integration.f90 spheroid_main.f90 -lfftw3 -lm -Wall -pedantic -pg

nohup的./p0 &

注意，2049 = 4096/2 + 1

當製作plan_backward，是不是正確，我們使用2049而不是4096，因爲輸出的維數是2049？

但是當我這樣做的時候，它就爆炸了。（爆炸意味着NAN錯誤）

如果我使用4096製作plan_backward，除了一些傅立葉係數異常大而不應該發生之外，一切都很好。

請幫助我正確使用Fortran中的FFTW。這個問題已經讓我感到很沮喪了。

來源

2016-04-27 Sah-moo Kim

我相信我的代碼非常小。並且我證實了這個代碼可以很好地通過繪製的結果xyz座標進行驗證。但傅立葉係數是不正常的。所以我懷疑使用FFTW –

有一個問題可能是dfftw_execute_dft_c2r可能會破壞輸入數組的內容，如page中所述。關鍵部分是

FFTW_PRESERVE_INPUT指定不在場轉換不得改變其輸入數組。這通常是默認值，除了c2r和hc2r（即複數到實數）變換，其中FFTW_DESTROY_INPUT是默認值。

我們可以驗證這一點，例如，由第二FFT之後通過修改@VladimirF示例代碼，使得其data_out到data_save第一FFT（R2C）調用之後保存，然後計算它們的差（C2R ）電話。因此，就OP代碼而言，在輸入第二個FFT（c2r）之前，將output1和output2保存到不同的陣列似乎更安全。

來源

2016-04-27 21:02:21 roygvib

即使我沒有這樣的知識解決它，我真的很感激你的答案。它真正幫助我理解這種現象。 –

@ Sah-mooKim嗨，我不知道這種特殊的行爲，因爲我以前從未使用過複雜到真實的FFT。我相信FFTW網站應該更加清楚地強調這一點...... – roygvib

首先，雖然你聲稱你的例子很少，但它仍然很大，我沒有時間研究它。

但我更新了我的代碼https://gist.github.com/LadaF/73eb430682ef527eea9972ceb96116c5以顯示反向變換並回答關於變換維度的標題問題。

變換的邏輯大小是實數組的大小（Real-data DFT Array Format），但複雜部分由於固有的對稱性而較小。

但是，當您將第一個r2c從大小爲n的真實數組轉換爲大小爲n/2+1的複雜數組時。然後進行相反的變換，真正的陣列應該再次大小爲n。

這是我從要點小例子：

module FFTW3 
    use, intrinsic :: iso_c_binding 
    include "fftw3.f03" 
end module 

    use FFTW3 
    implicit none 

    integer, parameter :: n = 100 

    real(c_double), allocatable :: data_in(:) 
    complex(c_double_complex), allocatable :: data_out(:) 
    type(c_ptr) :: planf, planb 

    allocate(data_in(n)) 
    allocate(data_out(n/2+1)) 

    call random_number(data_in) 

    planf = fftw_plan_dft_r2c_1d(size(data_in), data_in, data_out, FFTW_ESTIMATE+FFTW_UNALIGNED) 
    planb = fftw_plan_dft_c2r_1d(size(data_in), data_out, data_in, FFTW_ESTIMATE+FFTW_UNALIGNED) 

    print *, "real input:", real(data_in) 

    call fftw_execute_dft_r2c(planf, data_in, data_out) 

    print *, "result real part:", real(data_out) 

    print *, "result imaginary part:", aimag(data_out) 

    call fftw_execute_dft_c2r(planb, data_out, data_in) 

    print *, "real output:", real(data_in)/n 

    call fftw_destroy_plan(planf) 
    call fftw_destroy_plan(planb) 
end

注意，我用的是現代的Fortran接口。我不喜歡使用舊的。

來源

2016-04-27 17:45:39

時出現問題我解決了我的問題。經過逆FFT後，複數傅里葉變量發生了變化。它造成了我的問題。但我不知道爲什麼會發生這樣的錯誤... –

在Fortran中使用r2c和c2r FFTW時，前向和後向維度是否相同？

回答

相關問題