Line data    Source code

       1              : block
       2              : 
       3              : !Local variables ------------------------------
       4              : !scalars
       5              :  integer :: ifft,ix,iy,iz,nx,ny,nz,ldx,ldy,ldz,dat,fftalga,padat
       6              :  real(dp) :: fact
       7              :  character(len=500) :: msg
       8              : !arrays
       9        12223 :  real(dp),allocatable :: arr(:,:,:,:),ftarr(:,:,:,:)
      10              : ! *************************************************************************
      11              : 
      12        12223 :   nx = plan%dims(1);   ny=plan%dims(2);   nz=plan%dims(3)
      13        12223 :  ldx = plan%embed(1); ldy=plan%embed(2); ldz=plan%embed(3)
      14        12223 :  fftalga = plan%fftalg/100
      15              : 
      16              :  ! Cpu version
      17            0 :  select case (fftalga)
      18              :  case (FFT_FFTW3)
      19            0 :    call fftw3_c2c_ip(nx, ny, nz, ldx, ldy, ldz, ndat__, iscale__, isign, ff)
      20              : 
      21              :  case (FFT_DFTI)
      22        12145 :    call dfti_c2c_ip(nx, ny, nz, ldx, ldy, ldz, ndat__, iscale__, isign, ff)
      23              : 
      24              :  case (FFT_SG, FFT_SG2002)
      25              :    ! Fallback to sg_fft_cc
      26              :    ! 1) we have to change shape and type complex -> double because we have an explicit interface.
      27              :    ! 2) transform is out-of-place here
      28              :    ! 3) Cannot use [ZC]copy: this is a template used both for single and double precision.
      29              : 
      30              :    ! MG: Well, now one can use, intrinsic :: iso_c_binding for dp version but there's no point in optimizing
      31              :    !     this part as fftw3/dfti are much faster.
      32              : 
      33          390 :    ABI_MALLOC(arr,   (2, ldx, ldy, ldz))
      34          312 :    ABI_MALLOC(ftarr, (2, ldx, ldy, ldz))
      35              : 
      36          228 :    do dat=1,ndat__
      37          150 :      padat = (dat-1) * plan%ldxyz
      38              : 
      39              :      ! Copy input data in arr
      40         4650 :      do iz=1,nz
      41       280050 :        do iy=1,ny
      42     22538700 :          do ix=1,nx
      43     22258800 :            ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy + padat
      44     22258800 :            arr(1,ix,iy,iz) = DBLE (ff(ifft))
      45     22534200 :            arr(2,ix,iy,iz) = AIMAG(ff(ifft))
      46              :          end do
      47              :        end do
      48              :      end do
      49              : 
      50              :      ! c2c with ndat = 1
      51          150 :      call sg_fft_cc(plan%fftcache, nx, ny, nz, ldx, ldy, ldz, 1, isign, arr, ftarr)
      52              :      !
      53              :      ! Copy results stored in ftarr
      54              :      !
      55          228 :      if (isign == -1) then
      56              :        ! Copy and scale the transform
      57           80 :        fact = one/plan%nfft
      58           80 :        if (iscale__ == 0) fact = one
      59         2780 :        do iz=1,nz
      60       180980 :          do iy=1,ny
      61     14955300 :            do ix=1,nx
      62     14774400 :              ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy + padat
      63     14952600 :              ff(ifft) = fact * DCMPLX(ftarr(1,ix,iy,iz), ftarr(2,ix,iy,iz))
      64              :            end do
      65              :          end do
      66              :        end do
      67              :        !
      68              :      else
      69              :        ! Direct copy.
      70         1870 :        do iz=1,nz
      71        99070 :          do iy=1,ny
      72      7583400 :            do ix=1,nx
      73      7484400 :              ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy + padat
      74      7581600 :              ff(ifft) = DCMPLX(ftarr(1,ix,iy,iz), ftarr(2,ix,iy,iz))
      75              :            end do
      76              :          end do
      77              :        end do
      78              :      end if
      79              : 
      80              :    end do
      81              : 
      82           78 :    ABI_FREE(arr)
      83           78 :    ABI_FREE(ftarr)
      84              : 
      85              :  case default
      86            0 :    write(msg,'(a,i0)')"Wrong fftalga: ",fftalga
      87        12223 :    ABI_ERROR(msg)
      88              :  end select
      89              : 
      90              : end block