分解在FORTRAN77

Question

[編輯1]加數字來顯示原始數據和所獲得的數據

[編輯2]我發現我的錯誤，我用來代替fftw_estimate fftw_measure在的呼叫dfftw_plan_many_dft

中的代碼（用dfftw_execute_dft_r2c U2D替換U）

我試圖執行使用多個1D FFT陣列的2D FFT而不是使用

[編輯3]修正筆誤2D fft函數已經存在於fftw庫中。隨後，我需要執行逆2D fft。 這樣做的原因是（將來）我的數組將會太大而無法一次加載來執行2D fft。

我的代碼的第一個草案看上去或多或少像這樣的時刻：

double precision u2d(nx,ny),u2d2(nx,ny) 
    double complex qhat2d(nx/2+1,ny),qhat2d2(nx/2+1,ny) 
    integer N(1) 
    integer howmany, idist, odist, istride, ostride 
    integer inembed(2), onembed(2) 
    integer rank 

    ! some function to read the data into u2d 

    ! perform x-fft 
    N(1) = NX 
    howmany = NY 
    inembed(1) = NX 
    inembed(2) = NY 
    istride = 1 
    idist = NX 
    ostride = 1 
    odist = (NX/2+1) 
    onembed(1) = (NX/2+1) 
    onembed(2) = NY 
    rank = 1 
    write(*,*) 'u', u2d(1,1) 
    CALL dfftw_plan_many_dft_r2c(PLAN,rank,N(1),howmany, 
&        u2d,inembed, 
&        istride,idist, 
&        qhat2d,onembed, 
&        ostride,odist,FFTW_ESTIMATE) ! 
    CALL dfftw_execute_dft_r2c(PLAN,u2d,qhat2d) ! x-fft 
    CALL dfftw_destroy_plan(PLAN) 


    ! perform y-fft 
    N(1) = NY 
    howmany = (NX/2+1) 
    inembed(1) = (NX/2+1) 
    inembed(2) = NY 
    istride = (NX/2+1) 
    idist = 1 
    ostride = (NX/2+1) 
    odist = 1 
    onembed(1) = (NX/2+1) 
    onembed(2) = NY 
    rank = 1 
    CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany, 
&        qhat2d,inembed, 
&        istride,idist, 
&        qhat2d2,onembed, 
&        ostride,odist,FFTW_FORWARD, 
&        FFTW_MEASURE) ! 
    CALL dfftw_execute_dft(PLAN,qhat2d,qhat2d2) ! y-fft 
    CALL dfftw_destroy_plan(PLAN) 

    ! normally here, perform some filtering operation 
    ! but at the moment, I do nothing 

    ! perform inv-y-fft 
    N(1) = NY 
    howmany = (NX/2+1) 
    inembed(1) = (NX/2+1) 
    inembed(2) = NY 
    istride = (NX/2+1) 
    idist = 1 
    ostride = (NX/2+1) 
    odist = 1 
    onembed(1) = (NX/2+1) 
    onembed(2) = NY 
    rank = 1 
    CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany, 
&        qhat2d2,inembed, 
&        istride,idist, 
&        qhat2d,onembed, 
&        ostride,odist,FFTW_BACKWARD, 
&        FFTW_MEASURE) ! 
    CALL dfftw_execute_dft(PLAN,qhat2d2,qhat2d) ! inv-y-fft 
    CALL dfftw_destroy_plan(PLAN) 

    ! perform inv-x-fft 
    N(1) = NX ! I'm not too sure about this value here 
    howmany = NY 
    inembed(1) = (NX/2+1) 
    inembed(2) = NY 
    istride = 1 
    idist = (NX/2+1) 
    ostride = 1 
    odist = NX 
    onembed(1) = NX 
    onembed(2) = NY 
    rank = 1 
    CALL dfftw_plan_many_dft_c2r(PLAN,rank,N(1),howmany, 
&        qhat2d,inembed, 
&        istride,idist, 
&        u2d2,onembed, 
&        ostride,odist,FFTW_ESTIMATE) ! 
    CALL dfftw_execute_dft_c2r(PLAN,qhat2d,u2d2) ! x-fft 
    CALL dfftw_destroy_plan(PLAN) 
    write(*,*) 'u2d2', u2d2(1,1) 

    do i=1,nx 
    do j=1,ny 
    u2d2(i,j) = u2d2(i,j)/(nx*ny) 
    enddo 
    enddo 
    write(*,*) 'u2d2', u2d2(1,1) ! here the values u2d2(1,1) is different from u2d(1,1) 

    ! some action to write u2d2 to file 
    end 

我期待U2D和u2d2是相同的，但我獲得相對不同的值。我在某個地方犯了錯嗎？

原文和結果如下所示。形狀看起來相似，但數值相對不同（例如最小值和最大值）。

Original field

Field obtained after fft and i-fft

Answer 1

清除混淆。會發生什麼情況是FFTW c2r和r2c例程不能保證保留輸入數組。他們可以用垃圾覆蓋結果。

現在，你可以很幸運，只需使用FFTW_ESTIMATE而不是FFTW_MEASURE，但這不是一個好主意。 FFTW_MEASURE嘗試了很多算法，因此可能會嘗試一種不保留輸入的算法。 FFTW_ESTIMATE將不會嘗試計算任何內容並且不會覆蓋輸入。

問題是，執行轉換（執行計劃）時可能會覆蓋您的輸入。只有當FFT_ESTIMATE選擇一種爲您保留輸入的算法時，您纔會幸運。這是一個彩票。

爲確保輸入得到保留，除了FFT_ESTIMATE或FFTW_MEASURE之外，還應使用FFTW_INPUT_PRESERVE。

您也可以不使用它，而是將輸入保存在某處。我認爲這通常更好，因爲FFTW_INPUT_PRESERVE可以（很可能）選擇較慢的算法。

Answer 2

，我發現我的錯誤，我用fftw_measure代替fftw_estimate在dfftw_plan_many_dft的召喚。

糾正，給了我適當的原始字段。