00001
00002 module fft_3d
00003 integer,parameter,private:: nradix=6
00004 integer,dimension(nradix),parameter,private:: radix=(/6, 8, 3, 4, 2, 5/)
00005 integer,parameter,private:: nfmax = 30
00006 logical,parameter,private:: prefer4 = .false.
00007 private:: fact235,genw,fftv,bftv,fft_tp1,fft_tp,fft_tps1,fft_tps
00008 private:: fftrd2,bftrd2,fftrd3,bftrd3,fftrd4,bftrd4,fftrd5,bftrd5
00009 private:: fftrd6,bftrd6,fftrd8,bftrd8
00010
00011 type fft3_struct
00012 integer:: nx,ny,nz,nxyz,lx,ly,lz,lxyz
00013 integer:: nfx,nfy,nfz
00014 integer,dimension(nfmax):: facx,facy,facz
00015 real(kind=8),pointer,dimension(:,:):: wx,wy,wz
00016 end type fft3_struct
00017 contains
00018
00019 subroutine fact235(n,nf,fac)
00020 implicit none
00021 integer,intent(in):: n
00022 integer,intent(out):: nf,fac(nfmax)
00023
00024 integer:: p,q,r,j
00025 integer:: n2,n3,n4,n5,n6,n8
00026
00027 p = n
00028 fac(:) = 0
00029
00030 n2 = 0
00031 q = p/2
00032 r = p-2*q
00033 do while ( (q .gt. 0) .and. (r .eq. 0) )
00034 n2 = n2 + 1
00035 p = q
00036 q = p/2
00037 r = p-2*q
00038 end do
00039
00040 n3 = 0
00041 q = p/3
00042 r = p-3*q
00043 do while ( (q .gt. 0) .and. (r .eq. 0) )
00044 n3 = n3 + 1
00045 p = q
00046 q = p/3
00047 r = p-3*q
00048 end do
00049
00050 n5 = 0
00051 q = p/5
00052 r = p-5*q
00053 do while ( (q .gt. 0) .and. (r .eq. 0) )
00054 n5 = n5 + 1
00055 p = q
00056 q = p/5
00057 r = p-5*q
00058 end do
00059
00060 if (p .ne. 1) then
00061 write(6,*) n, ' can not be factorized with 2, 3, 5'
00062 stop
00063 end if
00064
00065 n6 = min(n2,n3)
00066 n2 = n2 - n6
00067 n3 = n3 - n6
00068 if (prefer4) then
00069 n4 = n2/2
00070 j = mod(n2,2)
00071 if ( (j .eq. 1) .and. (n4 .ge. 1) ) then
00072 n8 = 1
00073 n4 = n4 - 1
00074 else
00075 n8 = 0
00076 end if
00077 else
00078 n8 = n2/3
00079 j = mod(n2,3)
00080 if (j .eq. 2) then
00081 n4 = 1
00082 else if ( (j .eq. 1) .and. (n8 .ge. 1) ) then
00083 n8 = n8 - 1
00084 n4 = 2
00085 else
00086 n4 = 0
00087 end if
00088 end if
00089
00090 n2 = n2 - n4*2 - n8*3
00091 if (n2 .lt. 0) then
00092 write(6,*) n2, ' n2 become negative'
00093 stop
00094 end if
00095 if (n2+n3+n4+n5+n6+n8 .gt. nfmax) then
00096 write(6,*) 'too many factors'
00097 stop
00098 end if
00099
00100 p = 1
00101 nf = 0
00102 do j=1,n2
00103 nf = nf + 1
00104 fac(nf) = 2
00105 p = p*2
00106 end do
00107 do j=1,n3
00108 nf = nf + 1
00109 fac(nf) = 3
00110 p = p*3
00111 end do
00112 do j=1,n4
00113 nf = nf + 1
00114 fac(nf) = 4
00115 p = p*4
00116 end do
00117 do j=1,n5
00118 nf = nf + 1
00119 fac(nf) = 5
00120 p = p*5
00121 end do
00122 do j=1,n6
00123 nf = nf + 1
00124 fac(nf) = 6
00125 p = p*6
00126 end do
00127 do j=1,n8
00128 nf = nf + 1
00129 fac(nf) = 8
00130 p = p*8
00131 end do
00132
00133 if (p .ne. n) then
00134 write(6,*) p, ' is not equal to n = ', n
00135 stop
00136 end if
00137
00138 return
00139 end subroutine fact235
00140
00141
00142 subroutine genw(n,nf,fac,w)
00143 implicit none
00144 integer,intent(in):: n,nf
00145 integer:: fac(nfmax)
00146 real(kind=8),intent(out):: w(2,n)
00147
00148 real(kind=8),parameter:: tpi=2.0d0*3.14159265358979323846d0
00149 integer:: iw,darg,ir,rad,narg,k,p
00150 real(kind=8):: argh,argi,arg
00151
00152 iw = 1
00153 darg = n
00154 do ir = 1,nf
00155 argh = tpi/darg
00156 rad = fac(ir)
00157 narg = darg/rad
00158 do p = 0,narg-1
00159 do k = 1,rad-1
00160 argi = argh*k
00161 arg = argi*p
00162 w(1,iw) = cos(arg)
00163 w(2,iw) = -sin(arg)
00164 iw = iw + 1
00165 end do
00166 end do
00167 darg = darg/rad
00168 end do
00169 return
00170 end subroutine genw
00171
00172 subroutine fftrd2(nvk,np,w,xr,xi,yr,yi)
00173 implicit none
00174 integer,intent(in):: nvk,np
00175 real(kind=8),intent(in):: xr(nvk,np,2),xi(nvk,np,2)
00176 real(kind=8),intent(out):: yr(nvk,2,np),yi(nvk,2,np)
00177 real(kind=8),intent(in):: w(2,np)
00178
00179 real(kind=8):: wr1,wi1,tr,ti
00180 integer:: p,ivk
00181
00182 do p = 1,np
00183 wr1 = w(1,p)
00184 wi1 = w(2,p)
00185 do ivk = 1,nvk
00186 yr(ivk,1,p) = xr(ivk,p,1) + xr(ivk,p,2)
00187 yi(ivk,1,p) = xi(ivk,p,1) + xi(ivk,p,2)
00188
00189 tr = xr(ivk,p,1) - xr(ivk,p,2)
00190 ti = xi(ivk,p,1) - xi(ivk,p,2)
00191 yr(ivk,2,p) = wr1*tr - wi1*ti
00192 yi(ivk,2,p) = wr1*ti + wi1*tr
00193 end do
00194 end do
00195 return
00196 end subroutine fftrd2
00197
00198 subroutine bftrd2(nvk,np,w,xr,xi,yr,yi)
00199 implicit none
00200 integer,intent(in):: nvk,np
00201 real(kind=8),intent(in):: xr(nvk,np,2),xi(nvk,np,2)
00202 real(kind=8),intent(out):: yr(nvk,2,np),yi(nvk,2,np)
00203 real(kind=8),intent(in):: w(2,np)
00204
00205 real(kind=8):: wr1,wi1,tr,ti
00206 integer:: p,ivk
00207
00208 do p = 1,np
00209 wr1 = w(1,p)
00210 wi1 = w(2,p)
00211 do ivk = 1,nvk
00212 yr(ivk,1,p) = xr(ivk,p,1) + xr(ivk,p,2)
00213 yi(ivk,1,p) = xi(ivk,p,1) + xi(ivk,p,2)
00214
00215 tr = xr(ivk,p,1) - xr(ivk,p,2)
00216 ti = xi(ivk,p,1) - xi(ivk,p,2)
00217 yr(ivk,2,p) = wr1*tr + wi1*ti
00218 yi(ivk,2,p) = wr1*ti - wi1*tr
00219 end do
00220 end do
00221 return
00222 end subroutine bftrd2
00223
00224 subroutine fftrd3(nvk,np,w,xr,xi,yr,yi)
00225 implicit none
00226 integer,intent(in):: nvk,np
00227 real(kind=8),intent(in):: xr(nvk,np,3),xi(nvk,np,3)
00228 real(kind=8),intent(out):: yr(nvk,3,np),yi(nvk,3,np)
00229 real(kind=8),intent(in):: w(2,2*np)
00230
00231 real(kind=8),parameter:: sp3 = 0.866025403784438647d0
00232
00233 real(kind=8):: wr1,wi1,wr2,wi2,tr1,ti1,tr2,ti2,tr3,ti3,tr,ti
00234 integer:: p,ivk,iw
00235
00236 iw = 1
00237 do p = 1,np
00238 wr1 = w(1,iw)
00239 wi1 = w(2,iw)
00240 iw = iw+1
00241 wr2 = w(1,iw)
00242 wi2 = w(2,iw)
00243 iw = iw+1
00244 do ivk = 1,nvk
00245 tr1 = xr(ivk,p,2) + xr(ivk,p,3)
00246 ti1 = xi(ivk,p,2) + xi(ivk,p,3)
00247 tr2 = xr(ivk,p,1) - 0.5d0*tr1
00248 ti2 = xi(ivk,p,1) - 0.5d0*ti1
00249 tr3 = sp3*(xi(ivk,p,2) - xi(ivk,p,3))
00250 ti3 = -sp3*(xr(ivk,p,2) - xr(ivk,p,3))
00251 yr(ivk,1,p) = xr(ivk,p,1) + tr1
00252 yi(ivk,1,p) = xi(ivk,p,1) + ti1
00253 tr = tr2 + tr3
00254 ti = ti2 + ti3
00255 yr(ivk,2,p) = wr1*tr - wi1*ti
00256 yi(ivk,2,p) = wr1*ti + wi1*tr
00257 tr = tr2 - tr3
00258 ti = ti2 - ti3
00259 yr(ivk,3,p) = wr2*tr - wi2*ti
00260 yi(ivk,3,p) = wr2*ti + wi2*tr
00261 end do
00262 end do
00263 return
00264 end subroutine fftrd3
00265
00266 subroutine bftrd3(nvk,np,w,xr,xi,yr,yi)
00267 implicit none
00268 integer,intent(in):: nvk,np
00269 real(kind=8),intent(in):: xr(nvk,np,3),xi(nvk,np,3)
00270 real(kind=8),intent(out):: yr(nvk,3,np),yi(nvk,3,np)
00271 real(kind=8),intent(in):: w(2,2*np)
00272
00273 real(kind=8),parameter:: sp3 = 0.866025403784438647d0
00274
00275 real(kind=8):: wr1,wi1,wr2,wi2,tr1,ti1,tr2,ti2,tr3,ti3,tr,ti
00276 integer:: p,ivk,iw
00277
00278 iw = 1
00279 do p = 1,np
00280 wr1 = w(1,iw)
00281 wi1 = w(2,iw)
00282 iw = iw+1
00283 wr2 = w(1,iw)
00284 wi2 = w(2,iw)
00285 iw = iw+1
00286 do ivk = 1,nvk
00287 tr1 = xr(ivk,p,2) + xr(ivk,p,3)
00288 ti1 = xi(ivk,p,2) + xi(ivk,p,3)
00289 tr2 = xr(ivk,p,1) - 0.5d0*tr1
00290 ti2 = xi(ivk,p,1) - 0.5d0*ti1
00291 tr3 = -sp3*(xi(ivk,p,2) - xi(ivk,p,3))
00292 ti3 = sp3*(xr(ivk,p,2) - xr(ivk,p,3))
00293 yr(ivk,1,p) = xr(ivk,p,1) + tr1
00294 yi(ivk,1,p) = xi(ivk,p,1) + ti1
00295 tr = tr2 + tr3
00296 ti = ti2 + ti3
00297 yr(ivk,2,p) = wr1*tr + wi1*ti
00298 yi(ivk,2,p) = wr1*ti - wi1*tr
00299 tr = tr2 - tr3
00300 ti = ti2 - ti3
00301 yr(ivk,3,p) = wr2*tr + wi2*ti
00302 yi(ivk,3,p) = wr2*ti - wi2*tr
00303 end do
00304 end do
00305 return
00306 end subroutine bftrd3
00307
00308 subroutine fftrd4(nvk,np,w,xr,xi,yr,yi)
00309 implicit none
00310 integer,intent(in):: nvk,np
00311 real(kind=8),intent(in):: xr(nvk,np,4),xi(nvk,np,4)
00312 real(kind=8),intent(out):: yr(nvk,4,np),yi(nvk,4,np)
00313 real(kind=8),intent(in):: w(2,3*np)
00314
00315 real(kind=8):: wr1,wi1,wr2,wi2,wr3,wi3
00316 real(kind=8):: tr1,ti1,tr2,ti2,tr3,ti3,tr4,ti4,tr,ti
00317 integer:: p,ivk,iw
00318
00319 iw = 1
00320 do p = 1,np
00321 wr1 = w(1,iw)
00322 wi1 = w(2,iw)
00323 iw = iw+1
00324 wr2 = w(1,iw)
00325 wi2 = w(2,iw)
00326 iw = iw+1
00327 wr3 = w(1,iw)
00328 wi3 = w(2,iw)
00329 iw = iw+1
00330 do ivk = 1,nvk
00331 tr1 = xr(ivk,p,1) + xr(ivk,p,3)
00332 ti1 = xi(ivk,p,1) + xi(ivk,p,3)
00333 tr2 = xr(ivk,p,1) - xr(ivk,p,3)
00334 ti2 = xi(ivk,p,1) - xi(ivk,p,3)
00335 tr3 = xr(ivk,p,2) + xr(ivk,p,4)
00336 ti3 = xi(ivk,p,2) + xi(ivk,p,4)
00337 tr4 = xi(ivk,p,2) - xi(ivk,p,4)
00338 ti4 = xr(ivk,p,4) - xr(ivk,p,2)
00339 yr(ivk,1,p) = tr1 + tr3
00340 yi(ivk,1,p) = ti1 + ti3
00341 tr = tr2 + tr4
00342 ti = ti2 + ti4
00343 yr(ivk,2,p) = wr1*tr - wi1*ti
00344 yi(ivk,2,p) = wr1*ti + wi1*tr
00345 tr = tr1 - tr3
00346 ti = ti1 - ti3
00347 yr(ivk,3,p) = wr2*tr - wi2*ti
00348 yi(ivk,3,p) = wr2*ti + wi2*tr
00349 tr = tr2 - tr4
00350 ti = ti2 - ti4
00351 yr(ivk,4,p) = wr3*tr - wi3*ti
00352 yi(ivk,4,p) = wr3*ti + wi3*tr
00353 end do
00354 end do
00355 return
00356 end subroutine fftrd4
00357
00358 subroutine bftrd4(nvk,np,w,xr,xi,yr,yi)
00359 implicit none
00360 integer,intent(in):: nvk,np
00361 real(kind=8),intent(in):: xr(nvk,np,4),xi(nvk,np,4)
00362 real(kind=8),intent(out):: yr(nvk,4,np),yi(nvk,4,np)
00363 real(kind=8),intent(in):: w(2,3*np)
00364
00365 real(kind=8):: wr1,wi1,wr2,wi2,wr3,wi3
00366 real(kind=8):: tr1,ti1,tr2,ti2,tr3,ti3,tr4,ti4,tr,ti
00367 integer:: p,ivk,iw
00368
00369 iw = 1
00370 do p = 1,np
00371 wr1 = w(1,iw)
00372 wi1 = w(2,iw)
00373 iw = iw+1
00374 wr2 = w(1,iw)
00375 wi2 = w(2,iw)
00376 iw = iw+1
00377 wr3 = w(1,iw)
00378 wi3 = w(2,iw)
00379 iw = iw+1
00380 do ivk = 1,nvk
00381 tr1 = xr(ivk,p,1) + xr(ivk,p,3)
00382 ti1 = xi(ivk,p,1) + xi(ivk,p,3)
00383 tr2 = xr(ivk,p,1) - xr(ivk,p,3)
00384 ti2 = xi(ivk,p,1) - xi(ivk,p,3)
00385 tr3 = xr(ivk,p,2) + xr(ivk,p,4)
00386 ti3 = xi(ivk,p,2) + xi(ivk,p,4)
00387 tr4 = xi(ivk,p,4) - xi(ivk,p,2)
00388 ti4 = xr(ivk,p,2) - xr(ivk,p,4)
00389 yr(ivk,1,p) = tr1 + tr3
00390 yi(ivk,1,p) = ti1 + ti3
00391 tr = tr2 + tr4
00392 ti = ti2 + ti4
00393 yr(ivk,2,p) = wr1*tr + wi1*ti
00394 yi(ivk,2,p) = wr1*ti - wi1*tr
00395 tr = tr1 - tr3
00396 ti = ti1 - ti3
00397 yr(ivk,3,p) = wr2*tr + wi2*ti
00398 yi(ivk,3,p) = wr2*ti - wi2*tr
00399 tr = tr2 - tr4
00400 ti = ti2 - ti4
00401 yr(ivk,4,p) = wr3*tr + wi3*ti
00402 yi(ivk,4,p) = wr3*ti - wi3*tr
00403 end do
00404 end do
00405 return
00406 end subroutine bftrd4
00407
00408 subroutine fftrd5(nvk,np,w,xr,xi,yr,yi)
00409 implicit none
00410 integer,intent(in):: nvk,np
00411 real(kind=8),intent(in):: xr(nvk,np,5),xi(nvk,np,5)
00412 real(kind=8),intent(out):: yr(nvk,5,np),yi(nvk,5,np)
00413 real(kind=8),intent(in):: w(2,4*np)
00414
00415 real(kind=8),parameter:: sp25 = 0.951056516295153572d0
00416 real(kind=8),parameter:: sq54 = 0.559016994374947424d0
00417 real(kind=8),parameter:: ss = 0.618033988749894848d0
00418
00419 real(kind=8):: wr1,wi1,wr2,wi2,wr3,wi3,wr4,wi4
00420 real(kind=8):: tr1,ti1,tr2,ti2,tr3,ti3,tr4,ti4,tr5,ti5
00421 real(kind=8):: tr6,ti6,tr7,ti7,tr8,ti8,tr9,ti9,tr10,ti10,tr11,ti11
00422 real(kind=8):: tr,ti
00423 integer:: p,ivk,iw
00424
00425 do ivk = 1,nvk
00426 tr1 = xr(ivk,1,2) + xr(ivk,1,5)
00427 ti1 = xi(ivk,1,2) + xi(ivk,1,5)
00428 tr2 = xr(ivk,1,3) + xr(ivk,1,4)
00429 ti2 = xi(ivk,1,3) + xi(ivk,1,4)
00430 tr3 = sp25*(xr(ivk,1,2) - xr(ivk,1,5))
00431 ti3 = sp25*(xi(ivk,1,2) - xi(ivk,1,5))
00432 tr4 = sp25*(xr(ivk,1,3) - xr(ivk,1,4))
00433 ti4 = sp25*(xi(ivk,1,3) - xi(ivk,1,4))
00434 tr5 = tr1 + tr2
00435 ti5 = ti1 + ti2
00436 tr6 = sq54*(tr1 - tr2)
00437 ti6 = sq54*(ti1 - ti2)
00438 tr7 = xr(ivk,1,1) - 0.25d0*tr5
00439 ti7 = xi(ivk,1,1) - 0.25d0*ti5
00440 tr8 = tr7 + tr6
00441 ti8 = ti7 + ti6
00442 tr9 = tr7 - tr6
00443 ti9 = ti7 - ti6
00444 tr10 = ti3 + ss*ti4
00445 ti10 = - tr3 - ss*tr4
00446 tr11 = ss*ti3 - ti4
00447 ti11 = - ss*tr3 + tr4
00448 yr(ivk,1,1) = xr(ivk,1,1) + tr5
00449 yi(ivk,1,1) = xi(ivk,1,1) + ti5
00450 yr(ivk,2,1) = tr8 + tr10
00451 yi(ivk,2,1) = ti8 + ti10
00452 yr(ivk,3,1) = tr9 + tr11
00453 yi(ivk,3,1) = ti9 + ti11
00454 yr(ivk,4,1) = tr9 - tr11
00455 yi(ivk,4,1) = ti9 - ti11
00456 yr(ivk,5,1) = tr8 - tr10
00457 yi(ivk,5,1) = ti8 - ti10
00458 end do
00459
00460 iw = 5
00461 do p = 2,np
00462 wr1 = w(1,iw)
00463 wi1 = w(2,iw)
00464 iw = iw+1
00465 wr2 = w(1,iw)
00466 wi2 = w(2,iw)
00467 iw = iw+1
00468 wr3 = w(1,iw)
00469 wi3 = w(2,iw)
00470 iw = iw+1
00471 wr4 = w(1,iw)
00472 wi4 = w(2,iw)
00473 iw = iw+1
00474 do ivk = 1,nvk
00475 tr1 = xr(ivk,p,2) + xr(ivk,p,5)
00476 ti1 = xi(ivk,p,2) + xi(ivk,p,5)
00477 tr2 = xr(ivk,p,3) + xr(ivk,p,4)
00478 ti2 = xi(ivk,p,3) + xi(ivk,p,4)
00479 tr3 = sp25*(xr(ivk,p,2) - xr(ivk,p,5))
00480 ti3 = sp25*(xi(ivk,p,2) - xi(ivk,p,5))
00481 tr4 = sp25*(xr(ivk,p,3) - xr(ivk,p,4))
00482 ti4 = sp25*(xi(ivk,p,3) - xi(ivk,p,4))
00483 tr5 = tr1 + tr2
00484 ti5 = ti1 + ti2
00485 tr6 = sq54*(tr1 - tr2)
00486 ti6 = sq54*(ti1 - ti2)
00487 tr7 = xr(ivk,p,1) - 0.25d0*tr5
00488 ti7 = xi(ivk,p,1) - 0.25d0*ti5
00489 tr8 = tr7 + tr6
00490 ti8 = ti7 + ti6
00491 tr9 = tr7 - tr6
00492 ti9 = ti7 - ti6
00493 tr10 = ti3 + ss*ti4
00494 ti10 = - tr3 - ss*tr4
00495 tr11 = ss*ti3 - ti4
00496 ti11 = - ss*tr3 + tr4
00497 yr(ivk,1,p) = xr(ivk,p,1) + tr5
00498 yi(ivk,1,p) = xi(ivk,p,1) + ti5
00499 tr = tr8 + tr10
00500 ti = ti8 + ti10
00501 yr(ivk,2,p) = wr1*tr - wi1*ti
00502 yi(ivk,2,p) = wr1*ti + wi1*tr
00503 tr = tr9 + tr11
00504 ti = ti9 + ti11
00505 yr(ivk,3,p) = wr2*tr - wi2*ti
00506 yi(ivk,3,p) = wr2*ti + wi2*tr
00507 tr = tr9 - tr11
00508 ti = ti9 - ti11
00509 yr(ivk,4,p) = wr3*tr - wi3*ti
00510 yi(ivk,4,p) = wr3*ti + wi3*tr
00511 tr = tr8 - tr10
00512 ti = ti8 - ti10
00513 yr(ivk,5,p) = wr4*tr - wi4*ti
00514 yi(ivk,5,p) = wr4*ti + wi4*tr
00515 end do
00516 end do
00517 return
00518 end subroutine fftrd5
00519
00520 subroutine bftrd5(nvk,np,w,xr,xi,yr,yi)
00521 implicit none
00522 integer,intent(in):: nvk,np
00523 real(kind=8),intent(in):: xr(nvk,np,5),xi(nvk,np,5)
00524 real(kind=8),intent(out):: yr(nvk,5,np),yi(nvk,5,np)
00525 real(kind=8),intent(in):: w(2,4*np)
00526
00527 real(kind=8),parameter:: sp25 = 0.951056516295153572d0
00528 real(kind=8),parameter:: sq54 = 0.559016994374947424d0
00529 real(kind=8),parameter:: ss = 0.618033988749894848d0
00530
00531 real(kind=8):: wr1,wi1,wr2,wi2,wr3,wi3,wr4,wi4
00532 real(kind=8):: tr1,ti1,tr2,ti2,tr3,ti3,tr4,ti4,tr5,ti5
00533 real(kind=8):: tr6,ti6,tr7,ti7,tr8,ti8,tr9,ti9,tr10,ti10,tr11,ti11
00534 real(kind=8):: tr,ti
00535 integer:: p,ivk,iw
00536
00537 do ivk = 1,nvk
00538 tr1 = xr(ivk,1,2) + xr(ivk,1,5)
00539 ti1 = xi(ivk,1,2) + xi(ivk,1,5)
00540 tr2 = xr(ivk,1,3) + xr(ivk,1,4)
00541 ti2 = xi(ivk,1,3) + xi(ivk,1,4)
00542 tr3 = sp25*(xr(ivk,1,2) - xr(ivk,1,5))
00543 ti3 = sp25*(xi(ivk,1,2) - xi(ivk,1,5))
00544 tr4 = sp25*(xr(ivk,1,3) - xr(ivk,1,4))
00545 ti4 = sp25*(xi(ivk,1,3) - xi(ivk,1,4))
00546 tr5 = tr1 + tr2
00547 ti5 = ti1 + ti2
00548 tr6 = sq54*(tr1 - tr2)
00549 ti6 = sq54*(ti1 - ti2)
00550 tr7 = xr(ivk,1,1) - 0.25d0*tr5
00551 ti7 = xi(ivk,1,1) - 0.25d0*ti5
00552 tr8 = tr7 + tr6
00553 ti8 = ti7 + ti6
00554 tr9 = tr7 - tr6
00555 ti9 = ti7 - ti6
00556 tr10 = - ti3 - ss*ti4
00557 ti10 = tr3 + ss*tr4
00558 tr11 = - ss*ti3 + ti4
00559 ti11 = ss*tr3 - tr4
00560 yr(ivk,1,1) = xr(ivk,1,1) + tr5
00561 yi(ivk,1,1) = xi(ivk,1,1) + ti5
00562 yr(ivk,2,1) = tr8 + tr10
00563 yi(ivk,2,1) = ti8 + ti10
00564 yr(ivk,3,1) = tr9 + tr11
00565 yi(ivk,3,1) = ti9 + ti11
00566 yr(ivk,4,1) = tr9 - tr11
00567 yi(ivk,4,1) = ti9 - ti11
00568 yr(ivk,5,1) = tr8 - tr10
00569 yi(ivk,5,1) = ti8 - ti10
00570 end do
00571
00572 iw = 5
00573 do p = 2,np
00574 wr1 = w(1,iw)
00575 wi1 = w(2,iw)
00576 iw = iw+1
00577 wr2 = w(1,iw)
00578 wi2 = w(2,iw)
00579 iw = iw+1
00580 wr3 = w(1,iw)
00581 wi3 = w(2,iw)
00582 iw = iw+1
00583 wr4 = w(1,iw)
00584 wi4 = w(2,iw)
00585 iw = iw+1
00586 do ivk = 1,nvk
00587 tr1 = xr(ivk,p,2) + xr(ivk,p,5)
00588 ti1 = xi(ivk,p,2) + xi(ivk,p,5)
00589 tr2 = xr(ivk,p,3) + xr(ivk,p,4)
00590 ti2 = xi(ivk,p,3) + xi(ivk,p,4)
00591 tr3 = sp25*(xr(ivk,p,2) - xr(ivk,p,5))
00592 ti3 = sp25*(xi(ivk,p,2) - xi(ivk,p,5))
00593 tr4 = sp25*(xr(ivk,p,3) - xr(ivk,p,4))
00594 ti4 = sp25*(xi(ivk,p,3) - xi(ivk,p,4))
00595 tr5 = tr1 + tr2
00596 ti5 = ti1 + ti2
00597 tr6 = sq54*(tr1 - tr2)
00598 ti6 = sq54*(ti1 - ti2)
00599 tr7 = xr(ivk,p,1) - 0.25d0*tr5
00600 ti7 = xi(ivk,p,1) - 0.25d0*ti5
00601 tr8 = tr7 + tr6
00602 ti8 = ti7 + ti6
00603 tr9 = tr7 - tr6
00604 ti9 = ti7 - ti6
00605 tr10 = - ti3 - ss*ti4
00606 ti10 = tr3 + ss*tr4
00607 tr11 = - ss*ti3 + ti4
00608 ti11 = ss*tr3 - tr4
00609 yr(ivk,1,p) = xr(ivk,p,1) + tr5
00610 yi(ivk,1,p) = xi(ivk,p,1) + ti5
00611 tr = tr8 + tr10
00612 ti = ti8 + ti10
00613 yr(ivk,2,p) = wr1*tr + wi1*ti
00614 yi(ivk,2,p) = wr1*ti - wi1*tr
00615 tr = tr9 + tr11
00616 ti = ti9 + ti11
00617 yr(ivk,3,p) = wr2*tr + wi2*ti
00618 yi(ivk,3,p) = wr2*ti - wi2*tr
00619 tr = tr9 - tr11
00620 ti = ti9 - ti11
00621 yr(ivk,4,p) = wr3*tr + wi3*ti
00622 yi(ivk,4,p) = wr3*ti - wi3*tr
00623 tr = tr8 - tr10
00624 ti = ti8 - ti10
00625 yr(ivk,5,p) = wr4*tr + wi4*ti
00626 yi(ivk,5,p) = wr4*ti - wi4*tr
00627 end do
00628 end do
00629 return
00630 end subroutine bftrd5
00631
00632 subroutine fftrd6(nvk,np,w,xr,xi,yr,yi)
00633 implicit none
00634 integer,intent(in):: nvk,np
00635 real(kind=8),intent(in):: xr(nvk,np,6),xi(nvk,np,6)
00636 real(kind=8),intent(out):: yr(nvk,6,np),yi(nvk,6,np)
00637 real(kind=8),intent(in):: w(2,5*np)
00638
00639 real(kind=8),parameter:: sp3 = 0.866025403784438647d0
00640
00641 real(kind=8):: wr1,wi1,wr2,wi2,wr3,wi3,wr4,wi4,wr5,wi5
00642 real(kind=8):: tr1,ti1,tr2,ti2,tr3,ti3,tr,ti
00643 real(kind=8):: ur1,ui1,ur2,ui2,ur3,ui3
00644 real(kind=8):: vr1,vi1,vr2,vi2,vr3,vi3
00645 integer:: p,ivk,iw
00646
00647 do ivk = 1,nvk
00648 tr1 = xr(ivk,1,3) + xr(ivk,1,5)
00649 ti1 = xi(ivk,1,3) + xi(ivk,1,5)
00650 tr2 = xr(ivk,1,1) - 0.5d0*tr1
00651 ti2 = xi(ivk,1,1) - 0.5d0*ti1
00652 tr3 = -sp3*(xr(ivk,1,3) - xr(ivk,1,5))
00653 ti3 = -sp3*(xi(ivk,1,3) - xi(ivk,1,5))
00654 ur1 = xr(ivk,1,1) + tr1
00655 ui1 = xi(ivk,1,1) + ti1
00656 ur2 = tr2 - ti3
00657 ui2 = ti2 + tr3
00658 ur3 = tr2 + ti3
00659 ui3 = ti2 - tr3
00660
00661 tr1 = xr(ivk,1,6) + xr(ivk,1,2)
00662 ti1 = xi(ivk,1,6) + xi(ivk,1,2)
00663 tr2 = xr(ivk,1,4) - 0.5d0*tr1
00664 ti2 = xi(ivk,1,4) - 0.5d0*ti1
00665 tr3 = -sp3*(xr(ivk,1,6) - xr(ivk,1,2))
00666 ti3 = -sp3*(xi(ivk,1,6) - xi(ivk,1,2))
00667 vr1 = xr(ivk,1,4) + tr1
00668 vi1 = xi(ivk,1,4) + ti1
00669 vr2 = tr2 - ti3
00670 vi2 = ti2 + tr3
00671 vr3 = tr2 + ti3
00672 vi3 = ti2 - tr3
00673
00674 yr(ivk,1,1) = ur1 + vr1
00675 yi(ivk,1,1) = ui1 + vi1
00676 yr(ivk,5,1) = ur2 + vr2
00677 yi(ivk,5,1) = ui2 + vi2
00678 yr(ivk,3,1) = ur3 + vr3
00679 yi(ivk,3,1) = ui3 + vi3
00680 yr(ivk,4,1) = ur1 - vr1
00681 yi(ivk,4,1) = ui1 - vi1
00682 yr(ivk,2,1) = ur2 - vr2
00683 yi(ivk,2,1) = ui2 - vi2
00684 yr(ivk,6,1) = ur3 - vr3
00685 yi(ivk,6,1) = ui3 - vi3
00686 end do
00687
00688 iw = 6
00689 do p = 2,np
00690 wr1 = w(1,iw)
00691 wi1 = w(2,iw)
00692 iw = iw+1
00693 wr2 = w(1,iw)
00694 wi2 = w(2,iw)
00695 iw = iw+1
00696 wr3 = w(1,iw)
00697 wi3 = w(2,iw)
00698 iw = iw+1
00699 wr4 = w(1,iw)
00700 wi4 = w(2,iw)
00701 iw = iw+1
00702 wr5 = w(1,iw)
00703 wi5 = w(2,iw)
00704 iw = iw+1
00705 do ivk = 1,nvk
00706 tr1 = xr(ivk,p,3) + xr(ivk,p,5)
00707 ti1 = xi(ivk,p,3) + xi(ivk,p,5)
00708 tr2 = xr(ivk,p,1) - 0.5d0*tr1
00709 ti2 = xi(ivk,p,1) - 0.5d0*ti1
00710 tr3 = -sp3*(xr(ivk,p,3) - xr(ivk,p,5))
00711 ti3 = -sp3*(xi(ivk,p,3) - xi(ivk,p,5))
00712 ur1 = xr(ivk,p,1) + tr1
00713 ui1 = xi(ivk,p,1) + ti1
00714 ur2 = tr2 - ti3
00715 ui2 = ti2 + tr3
00716 ur3 = tr2 + ti3
00717 ui3 = ti2 - tr3
00718
00719 tr1 = xr(ivk,p,6) + xr(ivk,p,2)
00720 ti1 = xi(ivk,p,6) + xi(ivk,p,2)
00721 tr2 = xr(ivk,p,4) - 0.5d0*tr1
00722 ti2 = xi(ivk,p,4) - 0.5d0*ti1
00723 tr3 = -sp3*(xr(ivk,p,6) - xr(ivk,p,2))
00724 ti3 = -sp3*(xi(ivk,p,6) - xi(ivk,p,2))
00725 vr1 = xr(ivk,p,4) + tr1
00726 vi1 = xi(ivk,p,4) + ti1
00727 vr2 = tr2 - ti3
00728 vi2 = ti2 + tr3
00729 vr3 = tr2 + ti3
00730 vi3 = ti2 - tr3
00731
00732 yr(ivk,1,p) = ur1 + vr1
00733 yi(ivk,1,p) = ui1 + vi1
00734 tr = ur2 + vr2
00735 ti = ui2 + vi2
00736 yr(ivk,5,p) = wr4*tr - wi4*ti
00737 yi(ivk,5,p) = wr4*ti + wi4*tr
00738 tr = ur3 + vr3
00739 ti = ui3 + vi3
00740 yr(ivk,3,p) = wr2*tr - wi2*ti
00741 yi(ivk,3,p) = wr2*ti + wi2*tr
00742 tr = ur1 - vr1
00743 ti = ui1 - vi1
00744 yr(ivk,4,p) = wr3*tr - wi3*ti
00745 yi(ivk,4,p) = wr3*ti + wi3*tr
00746 tr = ur2 - vr2
00747 ti = ui2 - vi2
00748 yr(ivk,2,p) = wr1*tr - wi1*ti
00749 yi(ivk,2,p) = wr1*ti + wi1*tr
00750 tr = ur3 - vr3
00751 ti = ui3 - vi3
00752 yr(ivk,6,p) = wr5*tr - wi5*ti
00753 yi(ivk,6,p) = wr5*ti + wi5*tr
00754 end do
00755 end do
00756 return
00757 end subroutine fftrd6
00758
00759 subroutine bftrd6(nvk,np,w,xr,xi,yr,yi)
00760 implicit none
00761 integer,intent(in):: nvk,np
00762 real(kind=8),intent(in):: xr(nvk,np,6),xi(nvk,np,6)
00763 real(kind=8),intent(out):: yr(nvk,6,np),yi(nvk,6,np)
00764 real(kind=8),intent(in):: w(2,5*np)
00765
00766 real(kind=8),parameter:: sp3 = 0.866025403784438647d0
00767
00768 real(kind=8):: wr1,wi1,wr2,wi2,wr3,wi3,wr4,wi4,wr5,wi5
00769 real(kind=8):: tr1,ti1,tr2,ti2,tr3,ti3,tr,ti
00770 real(kind=8):: ur1,ui1,ur2,ui2,ur3,ui3
00771 real(kind=8):: vr1,vi1,vr2,vi2,vr3,vi3
00772 integer:: p,ivk,iw
00773
00774 do ivk = 1,nvk
00775 tr1 = xr(ivk,1,3) + xr(ivk,1,5)
00776 ti1 = xi(ivk,1,3) + xi(ivk,1,5)
00777 tr2 = xr(ivk,1,1) - 0.5d0*tr1
00778 ti2 = xi(ivk,1,1) - 0.5d0*ti1
00779 tr3 = sp3*(xr(ivk,1,3) - xr(ivk,1,5))
00780 ti3 = sp3*(xi(ivk,1,3) - xi(ivk,1,5))
00781 ur1 = xr(ivk,1,1) + tr1
00782 ui1 = xi(ivk,1,1) + ti1
00783 ur2 = tr2 - ti3
00784 ui2 = ti2 + tr3
00785 ur3 = tr2 + ti3
00786 ui3 = ti2 - tr3
00787
00788 tr1 = xr(ivk,1,6) + xr(ivk,1,2)
00789 ti1 = xi(ivk,1,6) + xi(ivk,1,2)
00790 tr2 = xr(ivk,1,4) - 0.5d0*tr1
00791 ti2 = xi(ivk,1,4) - 0.5d0*ti1
00792 tr3 = sp3*(xr(ivk,1,6) - xr(ivk,1,2))
00793 ti3 = sp3*(xi(ivk,1,6) - xi(ivk,1,2))
00794 vr1 = xr(ivk,1,4) + tr1
00795 vi1 = xi(ivk,1,4) + ti1
00796 vr2 = tr2 - ti3
00797 vi2 = ti2 + tr3
00798 vr3 = tr2 + ti3
00799 vi3 = ti2 - tr3
00800
00801 yr(ivk,1,1) = ur1 + vr1
00802 yi(ivk,1,1) = ui1 + vi1
00803 yr(ivk,5,1) = ur2 + vr2
00804 yi(ivk,5,1) = ui2 + vi2
00805 yr(ivk,3,1) = ur3 + vr3
00806 yi(ivk,3,1) = ui3 + vi3
00807 yr(ivk,4,1) = ur1 - vr1
00808 yi(ivk,4,1) = ui1 - vi1
00809 yr(ivk,2,1) = ur2 - vr2
00810 yi(ivk,2,1) = ui2 - vi2
00811 yr(ivk,6,1) = ur3 - vr3
00812 yi(ivk,6,1) = ui3 - vi3
00813 end do
00814
00815
00816 iw = 6
00817 do p = 2,np
00818 wr1 = w(1,iw)
00819 wi1 = w(2,iw)
00820 iw = iw+1
00821 wr2 = w(1,iw)
00822 wi2 = w(2,iw)
00823 iw = iw+1
00824 wr3 = w(1,iw)
00825 wi3 = w(2,iw)
00826 iw = iw+1
00827 wr4 = w(1,iw)
00828 wi4 = w(2,iw)
00829 iw = iw+1
00830 wr5 = w(1,iw)
00831 wi5 = w(2,iw)
00832 iw = iw+1
00833 do ivk = 1,nvk
00834 tr1 = xr(ivk,p,3) + xr(ivk,p,5)
00835 ti1 = xi(ivk,p,3) + xi(ivk,p,5)
00836 tr2 = xr(ivk,p,1) - 0.5d0*tr1
00837 ti2 = xi(ivk,p,1) - 0.5d0*ti1
00838 tr3 = sp3*(xr(ivk,p,3) - xr(ivk,p,5))
00839 ti3 = sp3*(xi(ivk,p,3) - xi(ivk,p,5))
00840 ur1 = xr(ivk,p,1) + tr1
00841 ui1 = xi(ivk,p,1) + ti1
00842 ur2 = tr2 - ti3
00843 ui2 = ti2 + tr3
00844 ur3 = tr2 + ti3
00845 ui3 = ti2 - tr3
00846
00847 tr1 = xr(ivk,p,6) + xr(ivk,p,2)
00848 ti1 = xi(ivk,p,6) + xi(ivk,p,2)
00849 tr2 = xr(ivk,p,4) - 0.5d0*tr1
00850 ti2 = xi(ivk,p,4) - 0.5d0*ti1
00851 tr3 = sp3*(xr(ivk,p,6) - xr(ivk,p,2))
00852 ti3 = sp3*(xi(ivk,p,6) - xi(ivk,p,2))
00853 vr1 = xr(ivk,p,4) + tr1
00854 vi1 = xi(ivk,p,4) + ti1
00855 vr2 = tr2 - ti3
00856 vi2 = ti2 + tr3
00857 vr3 = tr2 + ti3
00858 vi3 = ti2 - tr3
00859
00860 yr(ivk,1,p) = ur1 + vr1
00861 yi(ivk,1,p) = ui1 + vi1
00862 tr = ur2 + vr2
00863 ti = ui2 + vi2
00864 yr(ivk,5,p) = wr4*tr + wi4*ti
00865 yi(ivk,5,p) = wr4*ti - wi4*tr
00866 tr = ur3 + vr3
00867 ti = ui3 + vi3
00868 yr(ivk,3,p) = wr2*tr + wi2*ti
00869 yi(ivk,3,p) = wr2*ti - wi2*tr
00870 tr = ur1 - vr1
00871 ti = ui1 - vi1
00872 yr(ivk,4,p) = wr3*tr + wi3*ti
00873 yi(ivk,4,p) = wr3*ti - wi3*tr
00874 tr = ur2 - vr2
00875 ti = ui2 - vi2
00876 yr(ivk,2,p) = wr1*tr + wi1*ti
00877 yi(ivk,2,p) = wr1*ti - wi1*tr
00878 tr = ur3 - vr3
00879 ti = ui3 - vi3
00880 yr(ivk,6,p) = wr5*tr + wi5*ti
00881 yi(ivk,6,p) = wr5*ti - wi5*tr
00882 end do
00883 end do
00884 return
00885 end subroutine bftrd6
00886
00887 subroutine fftrd8(nvk,np,w,xr,xi,yr,yi)
00888 implicit none
00889 integer,intent(in):: nvk,np
00890 real(kind=8),intent(in):: xr(nvk,np,8),xi(nvk,np,8)
00891 real(kind=8),intent(out):: yr(nvk,8,np),yi(nvk,8,np)
00892 real(kind=8),intent(in):: w(2,7*np)
00893
00894 real(kind=8),parameter:: sp2 = 0.70710678118654752440d0
00895
00896 real(kind=8):: wr1,wi1,wr2,wi2,wr3,wi3,wr4,wi4,wr5,wi5,wr6,wi6,wr7,wi7
00897 real(kind=8):: tr1,ti1,tr2,ti2,tr3,ti3,tr4,ti4,tr,ti,sr,si
00898 real(kind=8):: ur1,ui1,ur2,ui2,ur3,ui3,ur4,ui4
00899 real(kind=8):: vr1,vi1,vr2,vi2,vr3,vi3,vr4,vi4
00900 integer:: p,ivk,iw
00901
00902 do ivk = 1,nvk
00903 tr1 = xr(ivk,1,1) + xr(ivk,1,5)
00904 ti1 = xi(ivk,1,1) + xi(ivk,1,5)
00905 tr2 = xr(ivk,1,1) - xr(ivk,1,5)
00906 ti2 = xi(ivk,1,1) - xi(ivk,1,5)
00907 tr3 = xr(ivk,1,3) + xr(ivk,1,7)
00908 ti3 = xi(ivk,1,3) + xi(ivk,1,7)
00909 tr4 = xi(ivk,1,3) - xi(ivk,1,7)
00910 ti4 = xr(ivk,1,7) - xr(ivk,1,3)
00911 ur1 = tr1 + tr3
00912 ui1 = ti1 + ti3
00913 ur2 = tr2 + tr4
00914 ui2 = ti2 + ti4
00915 ur3 = tr1 - tr3
00916 ui3 = ti1 - ti3
00917 ur4 = tr2 - tr4
00918 ui4 = ti2 - ti4
00919 tr1 = xr(ivk,1,2) + xr(ivk,1,6)
00920 ti1 = xi(ivk,1,2) + xi(ivk,1,6)
00921 tr2 = xr(ivk,1,2) - xr(ivk,1,6)
00922 ti2 = xi(ivk,1,2) - xi(ivk,1,6)
00923 tr3 = xr(ivk,1,4) + xr(ivk,1,8)
00924 ti3 = xi(ivk,1,4) + xi(ivk,1,8)
00925 tr4 = xi(ivk,1,4) - xi(ivk,1,8)
00926 ti4 = xr(ivk,1,8) - xr(ivk,1,4)
00927 vr1 = tr1 + tr3
00928 vi1 = ti1 + ti3
00929 vr2 = tr2 + tr4
00930 vi2 = ti2 + ti4
00931 vr3 = tr1 - tr3
00932 vi3 = ti1 - ti3
00933 vr4 = tr2 - tr4
00934 vi4 = ti2 - ti4
00935 yr(ivk,1,1) = ur1 + vr1
00936 yi(ivk,1,1) = ui1 + vi1
00937 yr(ivk,5,1) = ur1 - vr1
00938 yi(ivk,5,1) = ui1 - vi1
00939
00940 sr = sp2*(vr2 + vi2)
00941 si = sp2*(vi2 - vr2)
00942 yr(ivk,2,1) = ur2 + sr
00943 yi(ivk,2,1) = ui2 + si
00944 yr(ivk,6,1) = ur2 - sr
00945 yi(ivk,6,1) = ui2 - si
00946
00947 yr(ivk,3,1) = ur3 + vi3
00948 yi(ivk,3,1) = ui3 - vr3
00949 yr(ivk,7,1) = ur3 - vi3
00950 yi(ivk,7,1) = ui3 + vr3
00951
00952 sr = sp2*(vi4 - vr4)
00953 si = -sp2*(vr4 + vi4)
00954 yr(ivk,4,1) = ur4 + sr
00955 yi(ivk,4,1) = ui4 + si
00956 yr(ivk,8,1) = ur4 - sr
00957 yi(ivk,8,1) = ui4 - si
00958 end do
00959
00960 iw = 8
00961 do p = 2,np
00962 wr1 = w(1,iw)
00963 wi1 = w(2,iw)
00964 iw = iw+1
00965 wr2 = w(1,iw)
00966 wi2 = w(2,iw)
00967 iw = iw+1
00968 wr3 = w(1,iw)
00969 wi3 = w(2,iw)
00970 iw = iw+1
00971 wr4 = w(1,iw)
00972 wi4 = w(2,iw)
00973 iw = iw+1
00974 wr5 = w(1,iw)
00975 wi5 = w(2,iw)
00976 iw = iw+1
00977 wr6 = w(1,iw)
00978 wi6 = w(2,iw)
00979 iw = iw+1
00980 wr7 = w(1,iw)
00981 wi7 = w(2,iw)
00982 iw = iw+1
00983 do ivk = 1,nvk
00984 tr1 = xr(ivk,p,1) + xr(ivk,p,5)
00985 ti1 = xi(ivk,p,1) + xi(ivk,p,5)
00986 tr2 = xr(ivk,p,1) - xr(ivk,p,5)
00987 ti2 = xi(ivk,p,1) - xi(ivk,p,5)
00988 tr3 = xr(ivk,p,3) + xr(ivk,p,7)
00989 ti3 = xi(ivk,p,3) + xi(ivk,p,7)
00990 tr4 = xi(ivk,p,3) - xi(ivk,p,7)
00991 ti4 = xr(ivk,p,7) - xr(ivk,p,3)
00992 ur1 = tr1 + tr3
00993 ui1 = ti1 + ti3
00994 ur2 = tr2 + tr4
00995 ui2 = ti2 + ti4
00996 ur3 = tr1 - tr3
00997 ui3 = ti1 - ti3
00998 ur4 = tr2 - tr4
00999 ui4 = ti2 - ti4
01000 tr1 = xr(ivk,p,2) + xr(ivk,p,6)
01001 ti1 = xi(ivk,p,2) + xi(ivk,p,6)
01002 tr2 = xr(ivk,p,2) - xr(ivk,p,6)
01003 ti2 = xi(ivk,p,2) - xi(ivk,p,6)
01004 tr3 = xr(ivk,p,4) + xr(ivk,p,8)
01005 ti3 = xi(ivk,p,4) + xi(ivk,p,8)
01006 tr4 = xi(ivk,p,4) - xi(ivk,p,8)
01007 ti4 = xr(ivk,p,8) - xr(ivk,p,4)
01008 vr1 = tr1 + tr3
01009 vi1 = ti1 + ti3
01010 vr2 = tr2 + tr4
01011 vi2 = ti2 + ti4
01012 vr3 = tr1 - tr3
01013 vi3 = ti1 - ti3
01014 vr4 = tr2 - tr4
01015 vi4 = ti2 - ti4
01016 yr(ivk,1,p) = ur1 + vr1
01017 yi(ivk,1,p) = ui1 + vi1
01018 tr = ur1 - vr1
01019 ti = ui1 - vi1
01020 yr(ivk,5,p) = wr4*tr - wi4*ti
01021 yi(ivk,5,p) = wr4*ti + wi4*tr
01022
01023 sr = sp2*(vr2 + vi2)
01024 si = sp2*(vi2 - vr2)
01025 tr = ur2 + sr
01026 ti = ui2 + si
01027 yr(ivk,2,p) = wr1*tr - wi1*ti
01028 yi(ivk,2,p) = wr1*ti + wi1*tr
01029 tr = ur2 - sr
01030 ti = ui2 - si
01031 yr(ivk,6,p) = wr5*tr - wi5*ti
01032 yi(ivk,6,p) = wr5*ti + wi5*tr
01033
01034 tr = ur3 + vi3
01035 ti = ui3 - vr3
01036 yr(ivk,3,p) = wr2*tr - wi2*ti
01037 yi(ivk,3,p) = wr2*ti + wi2*tr
01038 tr = ur3 - vi3
01039 ti = ui3 + vr3
01040 yr(ivk,7,p) = wr6*tr - wi6*ti
01041 yi(ivk,7,p) = wr6*ti + wi6*tr
01042
01043 sr = sp2*(vi4 - vr4)
01044 si = -sp2*(vr4 + vi4)
01045 tr = ur4 + sr
01046 ti = ui4 + si
01047 yr(ivk,4,p) = wr3*tr - wi3*ti
01048 yi(ivk,4,p) = wr3*ti + wi3*tr
01049 tr = ur4 - sr
01050 ti = ui4 - si
01051 yr(ivk,8,p) = wr7*tr - wi7*ti
01052 yi(ivk,8,p) = wr7*ti + wi7*tr
01053 end do
01054 end do
01055 return
01056 end subroutine fftrd8
01057
01058 subroutine bftrd8(nvk,np,w,xr,xi,yr,yi)
01059 implicit none
01060 integer,intent(in):: nvk,np
01061 real(kind=8),intent(in):: xr(nvk,np,8),xi(nvk,np,8)
01062 real(kind=8),intent(out):: yr(nvk,8,np),yi(nvk,8,np)
01063 real(kind=8),intent(in):: w(2,7*np)
01064
01065 real(kind=8),parameter:: sp2 = 0.70710678118654752440d0
01066
01067 real(kind=8):: wr1,wi1,wr2,wi2,wr3,wi3,wr4,wi4,wr5,wi5,wr6,wi6,wr7,wi7
01068 real(kind=8):: tr1,ti1,tr2,ti2,tr3,ti3,tr4,ti4,tr,ti,sr,si
01069 real(kind=8):: ur1,ui1,ur2,ui2,ur3,ui3,ur4,ui4
01070 real(kind=8):: vr1,vi1,vr2,vi2,vr3,vi3,vr4,vi4
01071 integer:: p,ivk,iw
01072
01073 do ivk = 1,nvk
01074 tr1 = xr(ivk,1,1) + xr(ivk,1,5)
01075 ti1 = xi(ivk,1,1) + xi(ivk,1,5)
01076 tr2 = xr(ivk,1,1) - xr(ivk,1,5)
01077 ti2 = xi(ivk,1,1) - xi(ivk,1,5)
01078 tr3 = xr(ivk,1,3) + xr(ivk,1,7)
01079 ti3 = xi(ivk,1,3) + xi(ivk,1,7)
01080 tr4 = xi(ivk,1,7) - xi(ivk,1,3)
01081 ti4 = xr(ivk,1,3) - xr(ivk,1,7)
01082 ur1 = tr1 + tr3
01083 ui1 = ti1 + ti3
01084 ur2 = tr2 + tr4
01085 ui2 = ti2 + ti4
01086 ur3 = tr1 - tr3
01087 ui3 = ti1 - ti3
01088 ur4 = tr2 - tr4
01089 ui4 = ti2 - ti4
01090 tr1 = xr(ivk,1,2) + xr(ivk,1,6)
01091 ti1 = xi(ivk,1,2) + xi(ivk,1,6)
01092 tr2 = xr(ivk,1,2) - xr(ivk,1,6)
01093 ti2 = xi(ivk,1,2) - xi(ivk,1,6)
01094 tr3 = xr(ivk,1,4) + xr(ivk,1,8)
01095 ti3 = xi(ivk,1,4) + xi(ivk,1,8)
01096 tr4 = xi(ivk,1,8) - xi(ivk,1,4)
01097 ti4 = xr(ivk,1,4) - xr(ivk,1,8)
01098 vr1 = tr1 + tr3
01099 vi1 = ti1 + ti3
01100 vr2 = tr2 + tr4
01101 vi2 = ti2 + ti4
01102 vr3 = tr1 - tr3
01103 vi3 = ti1 - ti3
01104 vr4 = tr2 - tr4
01105 vi4 = ti2 - ti4
01106 yr(ivk,1,1) = ur1 + vr1
01107 yi(ivk,1,1) = ui1 + vi1
01108 yr(ivk,5,1) = ur1 - vr1
01109 yi(ivk,5,1) = ui1 - vi1
01110
01111 sr = sp2*(vr2 - vi2)
01112 si = sp2*(vi2 + vr2)
01113 yr(ivk,2,1) = ur2 + sr
01114 yi(ivk,2,1) = ui2 + si
01115 yr(ivk,6,1) = ur2 - sr
01116 yi(ivk,6,1) = ui2 - si
01117
01118 yr(ivk,3,1) = ur3 - vi3
01119 yi(ivk,3,1) = ui3 + vr3
01120 yr(ivk,7,1) = ur3 + vi3
01121 yi(ivk,7,1) = ui3 - vr3
01122
01123 sr = -sp2*(vr4 + vi4)
01124 si = sp2*(vr4 - vi4)
01125 yr(ivk,4,1) = ur4 + sr
01126 yi(ivk,4,1) = ui4 + si
01127 yr(ivk,8,1) = ur4 - sr
01128 yi(ivk,8,1) = ui4 - si
01129 end do
01130
01131 iw = 8
01132 do p = 2,np
01133 wr1 = w(1,iw)
01134 wi1 = w(2,iw)
01135 iw = iw+1
01136 wr2 = w(1,iw)
01137 wi2 = w(2,iw)
01138 iw = iw+1
01139 wr3 = w(1,iw)
01140 wi3 = w(2,iw)
01141 iw = iw+1
01142 wr4 = w(1,iw)
01143 wi4 = w(2,iw)
01144 iw = iw+1
01145 wr5 = w(1,iw)
01146 wi5 = w(2,iw)
01147 iw = iw+1
01148 wr6 = w(1,iw)
01149 wi6 = w(2,iw)
01150 iw = iw+1
01151 wr7 = w(1,iw)
01152 wi7 = w(2,iw)
01153 iw = iw+1
01154 do ivk = 1,nvk
01155 tr1 = xr(ivk,p,1) + xr(ivk,p,5)
01156 ti1 = xi(ivk,p,1) + xi(ivk,p,5)
01157 tr2 = xr(ivk,p,1) - xr(ivk,p,5)
01158 ti2 = xi(ivk,p,1) - xi(ivk,p,5)
01159 tr3 = xr(ivk,p,3) + xr(ivk,p,7)
01160 ti3 = xi(ivk,p,3) + xi(ivk,p,7)
01161 tr4 = xi(ivk,p,7) - xi(ivk,p,3)
01162 ti4 = xr(ivk,p,3) - xr(ivk,p,7)
01163 ur1 = tr1 + tr3
01164 ui1 = ti1 + ti3
01165 ur2 = tr2 + tr4
01166 ui2 = ti2 + ti4
01167 ur3 = tr1 - tr3
01168 ui3 = ti1 - ti3
01169 ur4 = tr2 - tr4
01170 ui4 = ti2 - ti4
01171 tr1 = xr(ivk,p,2) + xr(ivk,p,6)
01172 ti1 = xi(ivk,p,2) + xi(ivk,p,6)
01173 tr2 = xr(ivk,p,2) - xr(ivk,p,6)
01174 ti2 = xi(ivk,p,2) - xi(ivk,p,6)
01175 tr3 = xr(ivk,p,4) + xr(ivk,p,8)
01176 ti3 = xi(ivk,p,4) + xi(ivk,p,8)
01177 tr4 = xi(ivk,p,8) - xi(ivk,p,4)
01178 ti4 = xr(ivk,p,4) - xr(ivk,p,8)
01179 vr1 = tr1 + tr3
01180 vi1 = ti1 + ti3
01181 vr2 = tr2 + tr4
01182 vi2 = ti2 + ti4
01183 vr3 = tr1 - tr3
01184 vi3 = ti1 - ti3
01185 vr4 = tr2 - tr4
01186 vi4 = ti2 - ti4
01187 yr(ivk,1,p) = ur1 + vr1
01188 yi(ivk,1,p) = ui1 + vi1
01189 tr = ur1 - vr1
01190 ti = ui1 - vi1
01191 yr(ivk,5,p) = wr4*tr + wi4*ti
01192 yi(ivk,5,p) = wr4*ti - wi4*tr
01193
01194 sr = sp2*(vr2 - vi2)
01195 si = sp2*(vi2 + vr2)
01196 tr = ur2 + sr
01197 ti = ui2 + si
01198 yr(ivk,2,p) = wr1*tr + wi1*ti
01199 yi(ivk,2,p) = wr1*ti - wi1*tr
01200 tr = ur2 - sr
01201 ti = ui2 - si
01202 yr(ivk,6,p) = wr5*tr + wi5*ti
01203 yi(ivk,6,p) = wr5*ti - wi5*tr
01204
01205 tr = ur3 - vi3
01206 ti = ui3 + vr3
01207 yr(ivk,3,p) = wr2*tr + wi2*ti
01208 yi(ivk,3,p) = wr2*ti - wi2*tr
01209 tr = ur3 + vi3
01210 ti = ui3 - vr3
01211 yr(ivk,7,p) = wr6*tr + wi6*ti
01212 yi(ivk,7,p) = wr6*ti - wi6*tr
01213
01214 sr = -sp2*(vr4 + vi4)
01215 si = sp2*(vr4 - vi4)
01216 tr = ur4 + sr
01217 ti = ui4 + si
01218 yr(ivk,4,p) = wr3*tr + wi3*ti
01219 yi(ivk,4,p) = wr3*ti - wi3*tr
01220 tr = ur4 - sr
01221 ti = ui4 - si
01222 yr(ivk,8,p) = wr7*tr + wi7*ti
01223 yi(ivk,8,p) = wr7*ti - wi7*tr
01224 end do
01225 end do
01226 return
01227 end subroutine bftrd8
01228
01229 subroutine fftv(nv,n,nf,fac,w,lxyz,dr,di,er,ei,ix)
01230 implicit none
01231 integer,intent(in):: nv,n,lxyz
01232 integer,intent(in):: nf,fac(nfmax)
01233 real(kind=8),intent(in):: w(2,n)
01234 real(kind=8),dimension(lxyz):: dr,di,er,ei
01235 integer:: ix
01236
01237 integer:: iw,nvk,np,uif,if
01238
01239 iw = 1
01240 nvk = nv
01241 np = n/fac(1)
01242 uif = mod(nf,2)
01243 do if = 1,nf-uif,2
01244 select case (fac(if))
01245 case (2)
01246 call fftrd2(nvk,np,w(1,iw),dr,di,er,ei)
01247 case (3)
01248 call fftrd3(nvk,np,w(1,iw),dr,di,er,ei)
01249 case (4)
01250 call fftrd4(nvk,np,w(1,iw),dr,di,er,ei)
01251 case (5)
01252 call fftrd5(nvk,np,w(1,iw),dr,di,er,ei)
01253 case (6)
01254 call fftrd6(nvk,np,w(1,iw),dr,di,er,ei)
01255 case (8)
01256 call fftrd8(nvk,np,w(1,iw),dr,di,er,ei)
01257 case default
01258 write(6,*)'unsupported factor'
01259 stop
01260 end select
01261 iw=iw+np*(fac(if)-1)
01262 nvk=nvk*fac(if)
01263 np=np/fac(if+1)
01264 select case (fac(if+1))
01265 case (2)
01266 call fftrd2(nvk,np,w(1,iw),er,ei,dr,di)
01267 case (3)
01268 call fftrd3(nvk,np,w(1,iw),er,ei,dr,di)
01269 case (4)
01270 call fftrd4(nvk,np,w(1,iw),er,ei,dr,di)
01271 case (5)
01272 call fftrd5(nvk,np,w(1,iw),er,ei,dr,di)
01273 case (6)
01274 call fftrd6(nvk,np,w(1,iw),er,ei,dr,di)
01275 case (8)
01276 call fftrd8(nvk,np,w(1,iw),er,ei,dr,di)
01277 case default
01278 write(6,*)'unsupported factor'
01279 stop
01280 end select
01281 iw=iw+np*(fac(if+1)-1)
01282 nvk=nvk*fac(if+1)
01283 if (if+2 .le. nf) then
01284 np=np/fac(if+2)
01285 end if
01286 end do
01287 if (uif .eq. 1) then
01288 select case (fac(nf))
01289 case (2)
01290 call fftrd2(nvk,np,w(1,iw),dr,di,er,ei)
01291 case (3)
01292 call fftrd3(nvk,np,w(1,iw),dr,di,er,ei)
01293 case (4)
01294 call fftrd4(nvk,np,w(1,iw),dr,di,er,ei)
01295 case (5)
01296 call fftrd5(nvk,np,w(1,iw),dr,di,er,ei)
01297 case (6)
01298 call fftrd6(nvk,np,w(1,iw),dr,di,er,ei)
01299 case (8)
01300 call fftrd8(nvk,np,w(1,iw),dr,di,er,ei)
01301 case default
01302 write(6,*)'unsupported factor'
01303 stop
01304 end select
01305 ix = - ix
01306 end if
01307
01308 return
01309 end subroutine fftv
01310
01311 subroutine bftv(nv,n,nf,fac,w,lxyz,dr,di,er,ei,ix)
01312 implicit none
01313 integer,intent(in):: nv,n,lxyz
01314 integer,intent(in):: nf,fac(nfmax)
01315 real(kind=8),intent(in):: w(2,n)
01316 real(kind=8),dimension(lxyz):: dr,di,er,ei
01317 integer:: ix
01318
01319 integer:: iw,nvk,np,uif,if
01320
01321 iw = 1
01322 nvk = nv
01323 np = n/fac(1)
01324 uif = mod(nf,2)
01325 do if = 1,nf-uif,2
01326 select case (fac(if))
01327 case (2)
01328 call bftrd2(nvk,np,w(1,iw),dr,di,er,ei)
01329 case (3)
01330 call bftrd3(nvk,np,w(1,iw),dr,di,er,ei)
01331 case (4)
01332 call bftrd4(nvk,np,w(1,iw),dr,di,er,ei)
01333 case (5)
01334 call bftrd5(nvk,np,w(1,iw),dr,di,er,ei)
01335 case (6)
01336 call bftrd6(nvk,np,w(1,iw),dr,di,er,ei)
01337 case (8)
01338 call bftrd8(nvk,np,w(1,iw),dr,di,er,ei)
01339 case default
01340 write(6,*)'unsupported factor'
01341 stop
01342 end select
01343 iw=iw+np*(fac(if)-1)
01344 nvk=nvk*fac(if)
01345 np=np/fac(if+1)
01346 select case(fac(if+1))
01347 case (2)
01348 call bftrd2(nvk,np,w(1,iw),er,ei,dr,di)
01349 case (3)
01350 call bftrd3(nvk,np,w(1,iw),er,ei,dr,di)
01351 case (4)
01352 call bftrd4(nvk,np,w(1,iw),er,ei,dr,di)
01353 case (5)
01354 call bftrd5(nvk,np,w(1,iw),er,ei,dr,di)
01355 case (6)
01356 call bftrd6(nvk,np,w(1,iw),er,ei,dr,di)
01357 case (8)
01358 call bftrd8(nvk,np,w(1,iw),er,ei,dr,di)
01359 case default
01360 write(6,*)'unsupported factor'
01361 stop
01362 end select
01363 iw=iw+np*(fac(if+1)-1)
01364 nvk=nvk*fac(if+1)
01365 if (if+2 .le. nf) then
01366 np=np/fac(if+2)
01367 end if
01368 end do
01369 if (uif .eq. 1) then
01370 select case (fac(nf))
01371 case (2)
01372 call bftrd2(nvk,np,w(1,iw),dr,di,er,ei)
01373 case (3)
01374 call bftrd3(nvk,np,w(1,iw),dr,di,er,ei)
01375 case (4)
01376 call bftrd4(nvk,np,w(1,iw),dr,di,er,ei)
01377 case (5)
01378 call bftrd5(nvk,np,w(1,iw),dr,di,er,ei)
01379 case (6)
01380 call bftrd6(nvk,np,w(1,iw),dr,di,er,ei)
01381 case (8)
01382 call bftrd8(nvk,np,w(1,iw),dr,di,er,ei)
01383 case default
01384 write(6,*)'unsupported factor'
01385 stop
01386 end select
01387 ix = - ix
01388 end if
01389
01390 return
01391 end subroutine bftv
01392
01393 subroutine fft_pad_zero(nx,ny,nz,lx,ly,lz,dat)
01394 implicit none
01395 integer,intent(in):: nx,ny,nz,lx,ly,lz
01396 real(kind=8):: dat(lx,ly,lz,2)
01397
01398 integer:: j,k
01399 do k=1,2
01400 do j=nz+1,lz
01401 dat(:,:,j,k) = 0.0d0
01402 end do
01403 do j=nx+1,lx
01404 dat(j,:,:,k) = 0.0d0
01405 end do
01406 do j=ny+1,ly
01407 dat(:,j,:,k) = 0.0d0
01408 end do
01409 end do
01410 return
01411 end subroutine fft_pad_zero
01412
01413 subroutine fft_tp1(lxy,lz,src,dst)
01414 implicit none
01415 integer,intent(in):: lxy,lz
01416 real(kind=8),intent(in):: src(lxy,lz)
01417 real(kind=8),intent(out):: dst(lz,lxy)
01418 integer:: i,j,ui,uj
01419 real(kind=8):: r11,r12,r13,r14
01420 real(kind=8):: r21,r22,r23,r24
01421 real(kind=8):: r31,r32,r33,r34
01422 real(kind=8):: r41,r42,r43,r44
01423
01424 ui = mod(lz,4)
01425 uj = mod(lxy,4)
01426
01427 do i=1,ui
01428 do j=1,lxy
01429 dst(i,j) = src(j,i)
01430 end do
01431 end do
01432
01433 do j=1,uj
01434 do i=ui+1,lz
01435 dst(i,j) = src(j,i)
01436 end do
01437 end do
01438 do i=ui+1,lz,4
01439 do j=uj+1,lxy,4
01440 r11 = src(j,i)
01441 r21 = src(j+1,i)
01442 r31 = src(j+2,i)
01443 r41 = src(j+3,i)
01444 r12 = src(j,i+1)
01445 r22 = src(j+1,i+1)
01446 r32 = src(j+2,i+1)
01447 r42 = src(j+3,i+1)
01448 r13 = src(j,i+2)
01449 r23 = src(j+1,i+2)
01450 r33 = src(j+2,i+2)
01451 r43 = src(j+3,i+2)
01452 r14 = src(j,i+3)
01453 r24 = src(j+1,i+3)
01454 r34 = src(j+2,i+3)
01455 r44 = src(j+3,i+3)
01456
01457 dst(i,j) = r11
01458 dst(i+1,j) = r12
01459 dst(i+2,j) = r13
01460 dst(i+3,j) = r14
01461 dst(i,j+1) = r21
01462 dst(i+1,j+1) = r22
01463 dst(i+2,j+1) = r23
01464 dst(i+3,j+1) = r24
01465 dst(i,j+2) = r31
01466 dst(i+1,j+2) = r32
01467 dst(i+2,j+2) = r33
01468 dst(i+3,j+2) = r34
01469 dst(i,j+3) = r41
01470 dst(i+1,j+3) = r42
01471 dst(i+2,j+3) = r43
01472 dst(i+3,j+3) = r44
01473 end do
01474 end do
01475 return
01476 end subroutine fft_tp1
01477
01478 subroutine fft_tp(lxy,lz,src,dst)
01479 implicit none
01480 integer,intent(in):: lxy,lz
01481 real(kind=8),intent(in):: src(lxy,lz,2)
01482 real(kind=8),intent(out):: dst(lz,lxy,2)
01483
01484 call fft_tp1(lxy,lz,src(1,1,1),dst(1,1,1))
01485 call fft_tp1(lxy,lz,src(1,1,2),dst(1,1,2))
01486 return
01487 end subroutine fft_tp
01488
01489 subroutine fft_tps1(lxy,lz,src,dst,scl)
01490 implicit none
01491 integer,intent(in):: lxy,lz
01492 real(kind=8),intent(in):: src(lxy,lz)
01493 real(kind=8),intent(out):: dst(lz,lxy)
01494 real(kind=8),intent(in):: scl
01495 integer:: i,j,ui,uj
01496 real(kind=8):: r11,r12,r13,r14
01497 real(kind=8):: r21,r22,r23,r24
01498 real(kind=8):: r31,r32,r33,r34
01499 real(kind=8):: r41,r42,r43,r44
01500
01501 ui = mod(lz,4)
01502 uj = mod(lxy,4)
01503
01504 do i=1,ui
01505 do j=1,lxy
01506 dst(i,j) = scl*src(j,i)
01507 end do
01508 end do
01509
01510 do j=1,uj
01511 do i=ui+1,lz
01512 dst(i,j) = scl*src(j,i)
01513 end do
01514 end do
01515 do i=ui+1,lz,4
01516 do j=uj+1,lxy,4
01517 r11 = src(j,i)
01518 r21 = src(j+1,i)
01519 r31 = src(j+2,i)
01520 r41 = src(j+3,i)
01521 r12 = src(j,i+1)
01522 r22 = src(j+1,i+1)
01523 r32 = src(j+2,i+1)
01524 r42 = src(j+3,i+1)
01525 r13 = src(j,i+2)
01526 r23 = src(j+1,i+2)
01527 r33 = src(j+2,i+2)
01528 r43 = src(j+3,i+2)
01529 r14 = src(j,i+3)
01530 r24 = src(j+1,i+3)
01531 r34 = src(j+2,i+3)
01532 r44 = src(j+3,i+3)
01533
01534 dst(i,j) = scl*r11
01535 dst(i+1,j) = scl*r12
01536 dst(i+2,j) = scl*r13
01537 dst(i+3,j) = scl*r14
01538 dst(i,j+1) = scl*r21
01539 dst(i+1,j+1) = scl*r22
01540 dst(i+2,j+1) = scl*r23
01541 dst(i+3,j+1) = scl*r24
01542 dst(i,j+2) = scl*r31
01543 dst(i+1,j+2) = scl*r32
01544 dst(i+2,j+2) = scl*r33
01545 dst(i+3,j+2) = scl*r34
01546 dst(i,j+3) = scl*r41
01547 dst(i+1,j+3) = scl*r42
01548 dst(i+2,j+3) = scl*r43
01549 dst(i+3,j+3) = scl*r44
01550 end do
01551 end do
01552 return
01553 end subroutine fft_tps1
01554
01555 subroutine fft_tps(lxy,lz,src,dst,scl)
01556 implicit none
01557 integer,intent(in):: lxy,lz
01558 real(kind=8),intent(in):: src(lxy,lz,2)
01559 real(kind=8),intent(out):: dst(lz,lxy,2)
01560 real(kind=8),intent(in):: scl
01561
01562 call fft_tps1(lxy,lz,src(1,1,1),dst(1,1,1),scl)
01563 call fft_tps1(lxy,lz,src(1,1,2),dst(1,1,2),scl)
01564 return
01565 end subroutine fft_tps
01566
01567 subroutine fft3_init(mx,my,mz,lmx,lmy,lmz,fs)
01568 implicit none
01569 integer,intent(in):: mx,my,mz,lmx,lmy,lmz
01570 type (fft3_struct):: fs
01571
01572 fs%nx = mx
01573 fs%ny = my
01574 fs%nz = mz
01575 fs%nxyz = mx*my*mz
01576 fs%lx = lmx
01577 fs%ly = lmy
01578 fs%lz = lmz
01579 fs%lxyz = lmx*lmy*lmz
01580 call fact235(mx,fs%nfx,fs%facx)
01581 call fact235(my,fs%nfy,fs%facy)
01582 call fact235(mz,fs%nfz,fs%facz)
01583 allocate(fs%wx(2,mx),fs%wy(2,my),fs%wz(2,mz))
01584 call genw(mx,fs%nfx,fs%facx,fs%wx)
01585 call genw(my,fs%nfy,fs%facy,fs%wy)
01586 call genw(mz,fs%nfz,fs%facz,fs%wz)
01587 return
01588 end subroutine fft3_init
01589
01590 subroutine fft3_done(fs)
01591 implicit none
01592 type (fft3_struct):: fs
01593
01594 deallocate(fs%wx,fs%wy,fs%wz)
01595 return
01596 end subroutine fft3_done
01597
01598 subroutine fft3_fw(fs,dat,tmp)
01599 implicit none
01600 type (fft3_struct):: fs
01601 real(kind=8):: dat(fs%lxyz,2)
01602 real(kind=8):: tmp(fs%lxyz,2)
01603 integer:: ix
01604 real(kind=8):: scl
01605
01606 call fft_pad_zero(fs%nx,fs%ny,fs%nz,fs%lx,fs%ly,fs%lz,tmp)
01607 ix = 1
01608
01609 call fftv(fs%lx*fs%ly, fs%nz, fs%nfz, fs%facz, fs%wz, &
01610 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01611
01612 if (ix .eq. 1) then
01613 call fft_tp(fs%lx*fs%ly, fs%lz, dat, tmp)
01614 ix = - ix
01615
01616 call fftv(fs%lz*fs%lx, fs%ny, fs%nfy, fs%facy, fs%wy, &
01617 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01618 else if (ix .eq. -1) then
01619 call fft_tp(fs%lx*fs%ly, fs%lz, tmp, dat)
01620 ix = - ix
01621
01622 call fftv(fs%lz*fs%lx, fs%ny, fs%nfy, fs%facy, fs%wy, &
01623 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01624 else
01625 write(6,*) 'fft3_fw: y, abs(ix) != 1'
01626 stop
01627 end if
01628
01629 if (ix .eq. 1) then
01630 call fft_tp(fs%lz*fs%lx, fs%ly, dat, tmp)
01631 ix = - ix
01632
01633 call fftv(fs%ly*fs%lz, fs%nx, fs%nfx, fs%facx, fs%wx, &
01634 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01635 else if (ix .eq. -1) then
01636 call fft_tp(fs%lz*fs%lx, fs%ly, tmp, dat)
01637 ix = - ix
01638
01639 call fftv(fs%ly*fs%lz, fs%nx, fs%nfx, fs%facx, fs%wx, &
01640 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01641 else
01642 write(6,*) 'fft3_fw: x, abs(ix) != 1'
01643 stop
01644 end if
01645
01646 scl = 1.0d0/fs%nxyz
01647 if (ix .eq. 1) then
01648 call fft_tp(fs%ly*fs%lz, fs%lx, dat, tmp)
01649
01650 dat = scl*tmp
01651 else if (ix .eq. -1) then
01652 call fft_tp(fs%ly*fs%lz, fs%lx, tmp, dat)
01653
01654 dat = scl*dat
01655 else
01656 write(6,*) 'fft3_fw: f, abs(ix) != 1'
01657 stop
01658 end if
01659 return
01660 end subroutine fft3_fw
01661
01662 subroutine fft3_bw(fs,dat,tmp)
01663 implicit none
01664 type (fft3_struct):: fs
01665 real(kind=8):: dat(fs%lxyz,2)
01666 real(kind=8):: tmp(fs%lxyz,2)
01667 integer:: ix
01668
01669 call fft_pad_zero(fs%nx,fs%ny,fs%nz,fs%lx,fs%ly,fs%lz,tmp)
01670 ix = 1
01671
01672 call bftv(fs%lx*fs%ly, fs%nz, fs%nfz, fs%facz, fs%wz, &
01673 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01674
01675 if (ix .eq. 1) then
01676 call fft_tp(fs%lx*fs%ly, fs%lz, dat, tmp)
01677 ix = - ix
01678
01679 call bftv(fs%lz*fs%lx, fs%ny, fs%nfy, fs%facy, fs%wy, &
01680 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01681 else if (ix .eq. -1) then
01682 call fft_tp(fs%lx*fs%ly, fs%lz, tmp, dat)
01683 ix = - ix
01684
01685 call bftv(fs%lz*fs%lx, fs%ny, fs%nfy, fs%facy, fs%wy, &
01686 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01687 else
01688 write(6,*) 'fft3_bw: y, abs(ix) != 1'
01689 stop
01690 end if
01691
01692 if (ix .eq. 1) then
01693 call fft_tp(fs%lz*fs%lx, fs%ly, dat, tmp)
01694 ix = - ix
01695
01696 call bftv(fs%ly*fs%lz, fs%nx, fs%nfx, fs%facx, fs%wx, &
01697 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01698 else if (ix .eq. -1) then
01699 call fft_tp(fs%lz*fs%lx, fs%ly, tmp, dat)
01700 ix = - ix
01701
01702 call bftv(fs%ly*fs%lz, fs%nx, fs%nfx, fs%facx, fs%wx, &
01703 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01704 else
01705 write(6,*) 'fft3_bw: x, abs(ix) != 1'
01706 stop
01707 end if
01708
01709 if (ix .eq. 1) then
01710 call fft_tp(fs%ly*fs%lz, fs%lx, dat, tmp)
01711 dat = tmp
01712 else if (ix .eq. -1) then
01713 call fft_tp(fs%ly*fs%lz, fs%lx, tmp, dat)
01714 else
01715 write(6,*) 'fft3_bw: f, abs(ix) != 1'
01716 stop
01717 end if
01718 return
01719 end subroutine fft3_bw
01720
01721 subroutine fft3_fwsx(fs,dat,tmp,ix)
01722 implicit none
01723 type (fft3_struct):: fs
01724 real(kind=8):: dat(fs%lxyz,2)
01725 real(kind=8):: tmp(fs%lxyz,2)
01726 integer:: ix
01727 real(kind=8):: scl
01728
01729 call fft_pad_zero(fs%nx,fs%ny,fs%nz,fs%lx,fs%ly,fs%lz,tmp)
01730
01731 if (ix .eq. 1) then
01732 call fftv(fs%lx*fs%ly, fs%nz, fs%nfz, fs%facz, fs%wz, &
01733 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01734 else if (ix .eq. -1) then
01735 call fftv(fs%lx*fs%ly, fs%nz, fs%nfz, fs%facz, fs%wz, &
01736 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01737 else
01738 write(6,*) 'fft3_fwsx: z, abs(ix) != 1'
01739 stop
01740 end if
01741
01742 if (ix .eq. 1) then
01743 call fft_tp(fs%lx*fs%ly, fs%lz, dat, tmp)
01744 ix = - ix
01745
01746 call fftv(fs%lz*fs%lx, fs%ny, fs%nfy, fs%facy, fs%wy, &
01747 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01748 else if (ix .eq. -1) then
01749 call fft_tp(fs%lx*fs%ly, fs%lz, tmp, dat)
01750 ix = - ix
01751
01752 call fftv(fs%lz*fs%lx, fs%ny, fs%nfy, fs%facy, fs%wy, &
01753 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01754 else
01755 write(6,*) 'fft3_fwsx: y, abs(ix) != 1'
01756 stop
01757 end if
01758
01759 scl = 1.0d0/fs%nxyz
01760
01761 if (ix .eq. 1) then
01762 call fft_tps(fs%lz*fs%lx, fs%ly, dat, tmp, scl)
01763 ix = - ix
01764
01765 call fftv(fs%ly*fs%lz, fs%nx, fs%nfx, fs%facx, fs%wx, &
01766 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01767 else if (ix .eq. -1) then
01768 call fft_tps(fs%lz*fs%lx, fs%ly, tmp, dat, scl)
01769 ix = - ix
01770
01771 call fftv(fs%ly*fs%lz, fs%nx, fs%nfx, fs%facx, fs%wx, &
01772 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01773 else
01774 write(6,*) 'fft3_fwsx: x, abs(ix) != 1'
01775 stop
01776 end if
01777
01778 return
01779 end subroutine fft3_fwsx
01780
01781 subroutine fft3_bwsx(fs,dat,tmp,ix)
01782 implicit none
01783 type (fft3_struct):: fs
01784 real(kind=8):: dat(fs%lxyz,2)
01785 real(kind=8):: tmp(fs%lxyz,2)
01786 integer:: ix
01787
01788 call fft_pad_zero(fs%ny,fs%nz,fs%nx,fs%ly,fs%lz,fs%lx,tmp)
01789
01790 if (ix .eq. 1) then
01791 call bftv(fs%lz*fs%lx, fs%ny, fs%nfy, fs%facy, fs%wy, &
01792 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01793 else if (ix .eq. -1) then
01794 call bftv(fs%lz*fs%lx, fs%ny, fs%nfy, fs%facy, fs%wy, &
01795 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01796 else
01797 write(6,*) 'fft3_bwsx: x, abs(ix) != 1'
01798 stop
01799 end if
01800
01801 if (ix .eq. 1) then
01802 call fft_tp(fs%lz*fs%lx, fs%ly, dat, tmp)
01803 ix = - ix
01804
01805 call bftv(fs%ly*fs%lz, fs%nx, fs%nfx, fs%facx, fs%wx, &
01806 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01807 else if (ix .eq. -1) then
01808 call fft_tp(fs%lz*fs%lx, fs%ly, tmp, dat)
01809 ix = - ix
01810
01811 call bftv(fs%ly*fs%lz, fs%nx, fs%nfx, fs%facx, fs%wx, &
01812 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01813 else
01814 write(6,*) 'fft3_bwsx: y, abs(ix) != 1'
01815 stop
01816 end if
01817
01818 if (ix .eq. 1) then
01819 call fft_tp(fs%ly*fs%lz, fs%lx, dat, tmp)
01820 ix = - ix
01821
01822 call bftv(fs%lx*fs%ly, fs%nz, fs%nfz, fs%facz, fs%wz, &
01823 fs%lxyz, tmp(1,1), tmp(1,2), dat(1,1), dat(1,2), ix)
01824 else if (ix .eq. -1) then
01825 call fft_tp(fs%ly*fs%lz, fs%lx, tmp, dat)
01826 ix = - ix
01827
01828 call bftv(fs%lx*fs%ly, fs%nz, fs%nfz, fs%facz, fs%wz, &
01829 fs%lxyz, dat(1,1), dat(1,2), tmp(1,1), tmp(1,2), ix)
01830 else
01831 write(6,*) 'fft3_bwsx: z, abs(ix) != 1'
01832 stop
01833 end if
01834
01835 return
01836 end subroutine fft3_bwsx
01837
01838 end module fft_3d
