|
OpenMCTDHB v2.3
|
|
OpenMCTDHB v2.3
|
00001 ! vim:fdm=marker: 00002 !------------------------------------------------------------------------ 00014 00015 #if defined(FFT_FFTW3) 00016 MODULE modulefft 00017 IMPLICIT NONE 00018 00019 CONTAINS 00020 00021 !------------------------------------------------------------------------ 00026 SUBROUTINE get_FFT1D(instruction,length1,vec,INFO) 00027 !---{{{ 00028 IMPLICIT NONE 00029 INCLUDE 'fftw3.f' 00030 00031 INTEGER*8 :: plan 00032 INTEGER,INTENT(INOUT) :: INFO 00033 INTEGER,INTENT(IN) :: length1,instruction 00034 COMPLEX*16,INTENT(INOUT) :: vec(length1) 00035 00036 00037 do_what: SELECT CASE(instruction) 00038 CASE (1) 00039 CALL dfftw_plan_dft_1d(plan,length1,vec,vec,FFTW_BACKWARD,FFTW_ESTIMATE) 00040 CALL dfftw_execute_dft(plan,vec,vec) 00041 CASE (-1) 00042 CALL dfftw_plan_dft_1d(plan,length1,vec,vec,FFTW_FORWARD,FFTW_ESTIMATE) 00043 CALL dfftw_execute_dft(plan,vec,vec) 00044 CASE DEFAULT 00045 WRITE(*,*)"WRONG DIRECTION IN get_FFT1D (FFTW3), instruction",instruction 00046 STOP 00047 END SELECT do_what 00048 00049 CALL dfftw_destroy_plan(plan) 00050 INFO = 0 !< there is no INFO flag in dfftw_execute_dft 00051 !---}}} 00052 END SUBROUTINE get_FFT1D 00053 00054 !------------------------------------------------------------------------ 00059 SUBROUTINE get_FFT2D(instruction,length1,length2,vec,INFO) 00060 !---{{{ 00061 IMPLICIT NONE 00062 INCLUDE 'fftw3.f' 00063 00064 INTEGER*8 :: plan 00065 INTEGER,INTENT(INOUT) :: INFO 00066 INTEGER,INTENT(IN) :: length1,length2,instruction 00067 COMPLEX*16,INTENT(INOUT) :: vec(length1*length2) 00068 COMPLEX*16 :: vec_FFTW(length1,length2) 00069 00070 vec_FFTW = RESHAPE(vec,(/length1,length2/)) 00071 00072 do_what: SELECT CASE(instruction) 00073 CASE (1) 00074 CALL dfftw_plan_dft_2d(plan,length1,length2,vec_FFTW,vec_FFTW,FFTW_BACKWARD,FFTW_ESTIMATE) 00075 CALL dfftw_execute_dft(plan,vec_FFTW,vec_FFTW) 00076 CASE (-1) 00077 CALL dfftw_plan_dft_2d(plan,length1,length2,vec_FFTW,vec_FFTW,FFTW_FORWARD,FFTW_ESTIMATE) 00078 CALL dfftw_execute_dft(plan,vec_FFTW,vec_FFTW) 00079 CASE DEFAULT 00080 WRITE(*,*)"WRONG DIRECTION IN get_FFT1D (FFTW3), instruction",instruction 00081 STOP 00082 END SELECT do_what 00083 00084 CALL dfftw_destroy_plan(plan) 00085 INFO = 0 !< there is no INFO flag in dfftw_execute_dft 00086 00087 vec = RESHAPE(vec_FFTW,(/length1*length2/)) 00088 00089 !---}}} 00090 END SUBROUTINE get_FFT2D 00091 00092 !------------------------------------------------------------------------ 00097 SUBROUTINE get_FFT3D(instruction,length1,length2,length3,vec,INFO) 00098 !---{{{ 00099 IMPLICIT NONE 00100 INCLUDE 'fftw3.f' 00101 00102 INTEGER*8 :: plan 00103 INTEGER,INTENT(INOUT) :: INFO 00104 INTEGER,INTENT(IN) :: length1,length2,length3,instruction 00105 COMPLEX*16,INTENT(INOUT) :: vec(length1*length2*length3) 00106 COMPLEX*16 :: vec_FFTW(length1,length2,length3) 00107 00108 vec_FFTW = RESHAPE(vec,(/length1,length2,length3/)) 00109 00110 do_what: SELECT CASE(instruction) 00111 CASE (1) 00112 CALL dfftw_plan_dft_3d(plan,length1,length2,length3,vec_FFTW,vec_FFTW,FFTW_BACKWARD,FFTW_ESTIMATE) 00113 CALL dfftw_execute_dft(plan,vec_FFTW,vec_FFTW) 00114 CASE (-1) 00115 CALL dfftw_plan_dft_3d(plan,length1,length2,length3,vec_FFTW,vec_FFTW,FFTW_FORWARD,FFTW_ESTIMATE) 00116 CALL dfftw_execute_dft(plan,vec_FFTW,vec_FFTW) 00117 CASE DEFAULT 00118 WRITE(*,*)"WRONG DIRECTION IN get_FFT1D (FFTW3), instruction",instruction 00119 STOP 00120 END SELECT do_what 00121 00122 CALL dfftw_destroy_plan(plan) 00123 INFO = 0 !< there is no INFO flag in dfftw_execute_dft 00124 00125 vec = RESHAPE(vec_FFTW,(/length1*length2*length3/)) 00126 !---}}} 00127 END SUBROUTINE get_FFT3D 00128 00129 00130 END MODULE modulefft 00131 00132 #elif defined (FFT_ACML) 00133 00134 MODULE modulefft 00135 IMPLICIT NONE 00136 00137 CONTAINS 00138 00139 !------------------------------------------------------------------------ 00144 SUBROUTINE get_FFT1D(instruction,length1,vec,INFO) 00145 !---{{{ 00146 IMPLICIT NONE 00147 00148 INTEGER,INTENT(INOUT) :: INFO 00149 INTEGER,INTENT(IN) :: length1,instruction 00150 COMPLEX*16,INTENT(INOUT) :: vec(length1) 00151 COMPLEX*16 :: DUM(length1) 00152 00153 REAL*8 :: SCALEFAC 00154 COMPLEX*16 :: COMM(3*length1+100) 00155 LOGICAL :: INPL 00156 00157 00158 SCALEFAC = 1.D0 00159 INPL = .TRUE. 00160 00161 do_what: SELECT CASE(instruction) 00162 CASE (1) 00163 CALL ZFFT1DX(2,SCALEFAC,INPL,length1,vec,1,DUM,1,COMM,INFO) !< default initialization + bwd transform 00164 CASE (-1) 00165 CALL ZFFT1DX(-2,SCALEFAC,INPL,length1,vec,1,DUM,1,COMM,INFO) !< default initialization + fwd transform 00166 CASE DEFAULT 00167 WRITE(*,*)"ILLEGAL INSTRUCTION IN get_FFT1D (acml), instruction",instruction 00168 STOP 00169 END SELECT do_what 00170 00171 IF (INFO /= 0) THEN 00172 WRITE(*,*)"1D ZFFT1DX FAILED, INFO",INFO 00173 STOP 00174 END IF 00175 00176 00177 !---}}} 00178 END SUBROUTINE get_FFT1D 00179 00180 !------------------------------------------------------------------------ 00185 SUBROUTINE get_FFT2D(instruction,length1,length2,vec,INFO) 00186 !---{{{ 00187 IMPLICIT NONE 00188 00189 INTEGER,INTENT(INOUT) :: INFO 00190 INTEGER,INTENT(IN) :: length1,length2,instruction 00191 COMPLEX*16,INTENT(INOUT) :: vec(length1*length2) 00192 COMPLEX*16 :: DUM(length1*length2) 00193 00194 REAL*8 :: SCALEFAC 00195 COMPLEX*16 :: COMM(length1*length2+3*(length1+length2)+200) 00196 LOGICAL :: INPL,LTRANS 00197 INTEGER :: INCX1,INCX2,INCY1,INCY2 00198 00199 SCALEFAC = 1.D0 00200 INPL = .TRUE. 00201 LTRANS = .TRUE. 00202 INCX1 = 1 00203 INCX2 = length1 00204 INCY1 = 1 00205 INCY2 = 1 00206 00207 do_what: SELECT CASE(instruction) 00208 CASE (1) 00209 CALL ZFFT2DX(2,SCALEFAC,LTRANS,INPL,length1,length2,vec,INCX1,INCX2,DUM,INCY1,INCY2,COMM,INFO) !< default initialization + BWD transform 00210 CASE (-1) 00211 CALL ZFFT2DX(-2,SCALEFAC,LTRANS,INPL,length1,length2,vec,INCX1,INCX2,DUM,INCY1,INCY2,COMM,INFO) !< default initialization + FWD transform 00212 CASE DEFAULT 00213 WRITE(*,*)"ILLEGAL INSTRUCTION IN get_FFT2D (acml), instruction",instruction 00214 STOP 00215 END SELECT do_what 00216 00217 IF (INFO /= 0) THEN 00218 WRITE(*,*)"2D ZFFT2DX FAILED, INFO",INFO 00219 STOP 00220 END IF 00221 !---}}} 00222 END SUBROUTINE get_FFT2D 00223 00224 !------------------------------------------------------------------------ 00229 SUBROUTINE get_FFT3D(instruction,length1,length2,length3,vec,INFO) 00230 !---{{{ 00231 IMPLICIT NONE 00232 00233 INTEGER,INTENT(INOUT) :: INFO 00234 INTEGER,INTENT(IN) :: length1,length2,length3,instruction 00235 COMPLEX*16,INTENT(INOUT) :: vec(length1*length2*length3) 00236 COMPLEX*16 :: DUM(length1*length2*length3) 00237 00238 REAL*8 :: SCALEFAC 00239 COMPLEX*16 :: COMM(length1*length2*length3+4*(length1+length2+length3)+300) 00240 LOGICAL :: INPL 00241 INTEGER :: LCOMM,INCX1,INCX2,INCX3,INCY1,INCY2,INCY3 00242 00243 LCOMM = length1*length2*length3+4*(length1+length2+length3)+300 00244 SCALEFAC = 1.D0 00245 INPL = .TRUE. 00246 INCX1 = 1 00247 INCX2 = length1 00248 INCX3 = length1*length2 00249 INCY1 = 1 00250 INCY2 = 1 00251 INCY3 = 1 00252 00253 do_what: SELECT CASE(instruction) 00254 CASE (1) 00255 CALL ZFFT3DY(2,SCALEFAC,INPL,length1,length2,length3,vec,INCX1,INCX2,INCX3,DUM,INCY1,INCY2,INCY3,& 00256 COMM,LCOMM,INFO) !< default initialization + BWD transform 00257 CASE (-1) 00258 CALL ZFFT3DY(-2,SCALEFAC,INPL,length1,length2,length3,vec,INCX1,INCX2,INCX3,DUM,INCY1,INCY2,INCY3,& 00259 COMM,LCOMM,INFO) !< default initialization + FWD transform 00260 CASE DEFAULT 00261 WRITE(*,*)"ILLEGAL INSTRUCTION IN get_FFT3D (acml), instruction",instruction 00262 STOP 00263 END SELECT do_what 00264 00265 IF (INFO /= 0) THEN 00266 WRITE(*,*)"3D ZFFT3DY FAILED, INFO",INFO 00267 STOP 00268 END IF 00269 !---}}} 00270 END SUBROUTINE get_FFT3D 00271 00272 00273 END MODULE modulefft 00274 00275 #endif
1.7.3
