Changeset 1210 for palm/trunk/SOURCE/fft_xy.f90

Timestamp:

Aug 14, 2013 10:58:20 AM (11 years ago)

Author:

raasch

Message:

fftw support added; object file list in Makefile replaced by a short one line statement

File:

• palm/trunk/SOURCE/fft_xy.f90 (modified) (6 diffs)

palm/trunk/SOURCE/fft_xy.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! fftw added
 !
 ! Former revisions:
@@ -82 +82 @@
 #if defined( __cuda_fft )
     USE ISO_C_BINDING
+#elif defined( __fftw )
+    USE, INTRINSIC ::  ISO_C_BINDING
 #endif
     USE precision_kind
@@ -113 +115 @@
     INTEGER(C_INT), SAVE ::  plan_xf, plan_xi, plan_yf, plan_yi
     INTEGER, SAVE        ::  total_points_x_transpo, total_points_y_transpo
+#elif defined( __fftw )
+    INCLUDE  'fftw3.f03'
+    INTEGER(KIND=C_INT) ::  nx_c, ny_c
+    COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::  x_out, y_out
+    REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE            ::  x_in, y_in
+    TYPE(C_PTR), SAVE ::  plan_xf, plan_xi, plan_yf, plan_yi
 #endif
 
@@ -244 +252 @@
           CALL set99( trigs_y, ifax_y, ny+1 )
 
+       ELSEIF ( fft_method == 'fftw' )  THEN
+!
+!--       FFTW
+#if defined( __fftw )
+          nx_c = nx+1
+          ny_c = ny+1
+          ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2), &
+                    y_out(0:(ny+1)/2) )
+          plan_xf = FFTW_PLAN_DFT_R2C_1D( nx_c, x_in, x_out, FFTW_ESTIMATE )
+          plan_xi = FFTW_PLAN_DFT_C2R_1D( nx_c, x_out, x_in, FFTW_ESTIMATE )
+          plan_yf = FFTW_PLAN_DFT_R2C_1D( ny_c, y_in, y_out, FFTW_ESTIMATE )
+          plan_yi = FFTW_PLAN_DFT_C2R_1D( ny_c, y_out, y_in, FFTW_ESTIMATE )
+#else
+          message_string = 'preprocessor switch for fftw is missing'
+          CALL message( 'fft_init', 'PA0080', 1, 2, 0, 6, 0 )
+#endif
+
        ELSEIF ( fft_method == 'singleton-algorithm' )  THEN
 
@@ -413 +438 @@
           ENDIF
 
+       ELSEIF ( fft_method == 'fftw' )  THEN
+
+#if defined( __fftw )
+          IF ( forward_fft )  THEN
+
+             !$OMP PARALLEL PRIVATE ( work, i, j, k )
+             !$OMP DO
+             DO  k = nzb_x, nzt_x
+                DO  j = nys_x, nyn_x
+
+                   x_in(0:nx) = ar(0:nx,j,k)
+                   CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
+
+                   DO  i = 0, (nx+1)/2
+                      ar(i,j,k) = REAL( x_out(i) ) / ( nx+1 )
+                   ENDDO
+                   DO  i = 1, (nx+1)/2 - 1
+                      ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 )
+                   ENDDO
+
+                ENDDO
+             ENDDO
+             !$OMP END PARALLEL
+
+         ELSE
+             !$OMP PARALLEL PRIVATE ( work, i, j, k )
+             !$OMP DO
+             DO  k = nzb_x, nzt_x
+                DO  j = nys_x, nyn_x
+
+                   x_out(0) = CMPLX( ar(0,j,k), 0.0 )
+                   DO  i = 1, (nx+1)/2 - 1
+                      x_out(i)   = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
+                   ENDDO
+                   x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
+
+                   CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
+                   ar(0:nx,j,k) = x_in(0:nx)
+
+                ENDDO
+             ENDDO
+             !$OMP END PARALLEL
+
+         ENDIF
+#endif
+
        ELSEIF ( fft_method == 'system-specific' )  THEN
 
@@ -924 +995 @@
           ENDIF
 
+       ELSEIF ( fft_method == 'fftw' )  THEN
+
+#if defined( __fftw )
+          IF ( forward_fft )  THEN
+
+             !$OMP PARALLEL PRIVATE ( work, i, j, k )
+             !$OMP DO
+             DO  k = nzb_y, nzt_y
+                DO  i = nxl_y, nxr_y
+
+                   y_in(0:ny) = ar(0:ny,i,k)
+                   CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
+
+                   DO  j = 0, (ny+1)/2
+                      ar(j,i,k) = REAL( y_out(j) )  /(ny+1)
+                   ENDDO
+                   DO  j = 1, (ny+1)/2 - 1
+                      ar(ny+1-j,i,k) = AIMAG( y_out(j) )  /(ny+1)
+                   ENDDO
+
+                ENDDO
+             ENDDO
+             !$OMP END PARALLEL
+
+          ELSE
+
+             !$OMP PARALLEL PRIVATE ( work, i, j, k )
+             !$OMP DO
+             DO  k = nzb_y, nzt_y
+                DO  i = nxl_y, nxr_y
+
+                   y_out(0) = CMPLX( ar(0,i,k), 0.0 )
+                   DO  j = 1, (ny+1)/2 - 1
+                      y_out(j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
+                   ENDDO
+                   y_out((ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
+
+                   CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
+                   ar(0:ny,i,k) = y_in(0:ny)
+
+                ENDDO
+             ENDDO
+             !$OMP END PARALLEL
+
+          ENDIF
+#endif
+
        ELSEIF ( fft_method == 'system-specific' )  THEN