1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1999 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ################################################################# 9c ## ## 10c ## subroutine fftsetup -- setup 3-D Fast Fourier transform ## 11c ## ## 12c ################################################################# 13c 14c 15c "fftsetup" does initialization for a 3-D FFT to be computed 16c via either the FFTPACK or FFTW libraries 17c 18c 19 subroutine fftsetup 20 use fft 21 use openmp 22 use pme 23 implicit none 24 integer maxtable 25!$ integer ifront,iback 26!$ integer error,iguess 27c 28c 29c perform dynamic allocation of some global arrays 30c 31 allocate (qgrid(2,nfft1,nfft2,nfft3)) 32c 33c initialization of Fast Fourier transform using FFTW; 34c comment "dfftw_init_threads" and "dfftw_plan_with_threads" 35c calls if serial FFTW is used in place of OpenMP FFTW 36c 37!$ if (ffttyp .eq. 'FFTW') then 38!$ ifront = -1 39!$ iback = 1 40!$ error = 0 41!$ iguess = 0 42!$ call dfftw_init_threads (error) 43!$ call dfftw_plan_with_nthreads (nthread) 44!$ call dfftw_plan_dft_3d (planf,nfft1,nfft2,nfft3,qgrid, 45!$ & qgrid,ifront,iguess) 46!$ call dfftw_plan_dft_3d (planb,nfft1,nfft2,nfft3,qgrid, 47!$ & qgrid,iback,iguess) 48c 49c initialization of Fast Fourier transform using FFTPACK 50c 51!$ else 52 maxtable = 4 * max(nfft1,nfft2,nfft3) 53 if (allocated(ffttable)) then 54 if (size(ffttable) .ne. maxtable) deallocate (ffttable) 55 end if 56 if (.not. allocated(ffttable)) allocate (ffttable(maxtable,3)) 57 call cffti (nfft1,ffttable(1,1),iprime(1,1)) 58 call cffti (nfft2,ffttable(1,2),iprime(1,2)) 59 call cffti (nfft3,ffttable(1,3),iprime(1,3)) 60!$ end if 61 return 62 end 63c 64c 65c ############################################################### 66c ## ## 67c ## subroutine fftfront -- forward Fast Fourier transform ## 68c ## ## 69c ############################################################### 70c 71c 72c "fftfront" performs a 3-D FFT forward transform via a single 73c 3-D transform or three separate 1-D transforms 74c 75c 76 subroutine fftfront 77 use fft 78 use pme 79 implicit none 80 integer i,j,k 81 real*8, allocatable :: work(:,:) 82c 83c 84c perform a single 3-D forward transform using FFTW 85c 86!$ if (ffttyp .eq. 'FFTW') then 87!$ call dfftw_execute_dft (planf,qgrid,qgrid) 88!$ else 89c 90c perform three 1-D forward transforms using FFTPACK 91c 92 allocate (work(2,max(nfft1,nfft2,nfft3))) 93 do k = 1, nfft3 94 do j = 1, nfft2 95 do i = 1, nfft1 96 work(1,i) = qgrid(1,i,j,k) 97 work(2,i) = qgrid(2,i,j,k) 98 end do 99 call cfftf (nfft1,work,ffttable(1,1),iprime(1,1)) 100 do i = 1, nfft1 101 qgrid(1,i,j,k) = work(1,i) 102 qgrid(2,i,j,k) = work(2,i) 103 end do 104 end do 105 end do 106 do k = 1, nfft3 107 do i = 1, nfft1 108 do j = 1, nfft2 109 work(1,j) = qgrid(1,i,j,k) 110 work(2,j) = qgrid(2,i,j,k) 111 end do 112 call cfftf (nfft2,work,ffttable(1,2),iprime(1,2)) 113 do j = 1, nfft2 114 qgrid(1,i,j,k) = work(1,j) 115 qgrid(2,i,j,k) = work(2,j) 116 end do 117 end do 118 end do 119 do i = 1, nfft1 120 do j = 1, nfft2 121 do k = 1, nfft3 122 work(1,k) = qgrid(1,i,j,k) 123 work(2,k) = qgrid(2,i,j,k) 124 end do 125 call cfftf (nfft3,work,ffttable(1,3),iprime(1,3)) 126 do k = 1, nfft3 127 qgrid(1,i,j,k) = work(1,k) 128 qgrid(2,i,j,k) = work(2,k) 129 end do 130 end do 131 end do 132 deallocate (work) 133!$ end if 134 return 135 end 136c 137c 138c ############################################################### 139c ## ## 140c ## subroutine fftback -- backward Fast Fourier transform ## 141c ## ## 142c ############################################################### 143c 144c 145c "fftback" performs a 3-D FFT backward transform via a single 146c 3-D transform or three separate 1-D transforms 147c 148c 149 subroutine fftback 150 use fft 151 use pme 152 implicit none 153 integer i,j,k 154 real*8, allocatable :: work(:,:) 155c 156c 157c perform a single 3-D backward transform using FFTW 158c 159!$ if (ffttyp .eq. 'FFTW') then 160!$ call dfftw_execute_dft (planb,qgrid,qgrid) 161!$ else 162c 163c perform three 1-D backward transforms using FFTPACK 164c 165 allocate (work(2,max(nfft1,nfft2,nfft3))) 166 do k = 1, nfft3 167 do j = 1, nfft2 168 do i = 1, nfft1 169 work(1,i) = qgrid(1,i,j,k) 170 work(2,i) = qgrid(2,i,j,k) 171 end do 172 call cfftb (nfft1,work,ffttable(1,1),iprime(1,1)) 173 do i = 1, nfft1 174 qgrid(1,i,j,k) = work(1,i) 175 qgrid(2,i,j,k) = work(2,i) 176 end do 177 end do 178 end do 179 do k = 1, nfft3 180 do i = 1, nfft1 181 do j = 1, nfft2 182 work(1,j) = qgrid(1,i,j,k) 183 work(2,j) = qgrid(2,i,j,k) 184 end do 185 call cfftb (nfft2,work,ffttable(1,2),iprime(1,2)) 186 do j = 1, nfft2 187 qgrid(1,i,j,k) = work(1,j) 188 qgrid(2,i,j,k) = work(2,j) 189 end do 190 end do 191 end do 192 do i = 1, nfft1 193 do j = 1, nfft2 194 do k = 1, nfft3 195 work(1,k) = qgrid(1,i,j,k) 196 work(2,k) = qgrid(2,i,j,k) 197 end do 198 call cfftb (nfft3,work,ffttable(1,3),iprime(1,3)) 199 do k = 1, nfft3 200 qgrid(1,i,j,k) = work(1,k) 201 qgrid(2,i,j,k) = work(2,k) 202 end do 203 end do 204 end do 205 deallocate (work) 206!$ end if 207 return 208 end 209c 210c 211c ################################################################# 212c ## ## 213c ## subroutine fftclose -- close 3-D Fast Fourier transform ## 214c ## ## 215c ################################################################# 216c 217c 218c "fftclose" does cleanup after performing a 3-D FFT by destroying 219c the FFTW plans for the forward and backward transforms 220c 221c 222 subroutine fftclose 223 use fft 224 use pme 225 implicit none 226c 227c 228c remove the FFTW plans to avoid a cumulative memory leak 229c 230!$ if (ffttyp .eq. 'FFTW') then 231!$ call dfftw_destroy_plan (planf) 232!$ call dfftw_destroy_plan (planb) 233!$ end if 234c 235c perform deallocation of some global arrays 236c 237 if (allocated(qgrid)) deallocate (qgrid) 238 return 239 end 240