Bei Verwendung von r2c und c2r FFTW in Fortran sind die Vorwärts- und Rückwärtsabmessungen gleich?

Question

Schlag ist eine Hauptdatei

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 

Unten finden Sie eine Unterprogramm-Datei für die Integration ist

! 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 

Unten finden Sie eine Moduldatei

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 

Im Folgenden ein Unterprogramm Info

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 

Datei wird

ich sie mit dem folgenden Befehl kompiliert

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 &

Beachten Sie, dass 2049 = 4096/2 + 1

Wenn plan_backward machen, ist es nicht richtig, dass wir verwenden 2049 statt 4096, da die Dimension der Ausgabe 2049 ist?

Aber wenn ich das mache, explodiert es. (Blowing up bedeutet NAN Fehler)

Wenn ich 4096 bei der Erstellung von plan_backward verwende, ist alles in Ordnung, außer dass einige Fourier-Koeffizienten ungewöhnlich groß sind, was nicht passieren sollte.

Bitte helfen Sie mir, FFTW in Fortran richtig zu verwenden. Dieses Thema hat mich für eine lange Zeit entmutigt.

Answer 1

Ein Problem kann sein, dass dfftw_execute_dft_c2r den Inhalt des Eingangsarrays zerstören kann, wie in diesem page beschrieben. Der Schlüsselauszug lautet

FFTW_PRESERVE_INPUT gibt an, dass eine Out-of-Place-Transformation das Eingabearray nicht ändern darf. Dies ist normalerweise die Standardeinstellung, außer für c2r- und hc2r-Transformationen (d. H. Komplex zu reell), für die FFTW_DESTROY_INPUT die Standardeinstellung ist.

Wir dies überprüfen können, beispielsweise durch den Beispielcode durch @VladimirF Modifizieren, so daß es zu data_outdata_save direkt nach dem ersten FFT (R2C) -Aufruf speichert, und dann die Differenz nach der zweiten FFT-Berechnung (C2R) Anruf. Im Fall des OP-Codes scheint es also sicherer zu sein, output1 und output2 in verschiedenen Arrays zu speichern, bevor die zweite FFT (c2r) eingegeben wird.

Answer 2

Erstens, obwohl Sie behaupten, Ihr Beispiel ist minimal, es ist immer noch ziemlich groß, ich habe keine Zeit, es zu studieren.

Aber ich aktualisierte meine Gist-Code https://gist.github.com/LadaF/73eb430682ef527eea9972ceb96116c5, um auch die Rückwärts-Transformation zu zeigen und die Titelfrage über die Transformationsdimensionen zu beantworten.

Die logische Größe der Transformation ist die Größe des realen Arrays (Real-data DFT Array Format), aber der komplexe Teil ist aufgrund der inhärenten Symmetrien kleiner.

Aber wenn Sie zuerst r2c Transformation von realen Array der Größe n zu komplexen Array der Größe n/2+1 machen. und dann eine umgekehrte Rücktransformation, die reale Matrix sollte wieder die Größe n haben.

Das ist mein minimal Beispiel vom Kern:

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 

Bitte beachte, dass ich die moderne Fortran-Schnittstelle verwenden. Ich mag es nicht, das alte zu benutzen.