Line data    Source code

       1              : block
       2              : !Local variables-------------------------------
       3              : !scalars
       4              :  integer,parameter :: cplex0=0,option0=0,ndat1=1,iflag1=1,me_g0=1
       5              :  integer :: nx,ny,nz,ldx,ldy,ldz,fftalg,fftalga,fftalgc,fftcache
       6              :  integer :: idat,ig,rspad,gspad,ix,iy,iz,ifft
       7              :  real(dp),parameter :: weight1_r=one,weight1_i=one,xnorm1=one
       8              : !arrays
       9              :  integer :: shiftg(3),symm(3,3),dum_kg_kout(0,0)
      10              :  real(dp) :: dum_denpot(0,0,0),dum_fofgout(0,0)
      11       215047 :  real(dp),allocatable :: fofgin(:,:),fofr(:,:,:,:),ftarr(:,:,:,:)
      12              : 
      13              : ! *************************************************************************
      14              : 
      15       215047 :  fftalg = ngfft(7); fftcache = ngfft(8)
      16       215047 :  fftalga = fftalg/100; fftalgc = mod(fftalg, 10)
      17       215047 :  nx = ngfft(1); ny = ngfft(2); nz = ngfft(3)
      18              : 
      19            0 :  select case (fftalga)
      20              :  case (FFT_FFTW3)
      21            0 :    ldx = nx; ldy = ny; ldz = nz ! No augmentation, the caller does not support it.
      22              :    call fftw3_fftug(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, nspinor*ndat, istwf_k, mgfft, &
      23            0 :                     kg_k, gbound_k, ug, ur)
      24              : 
      25              :  case (FFT_DFTI)
      26       215006 :    ldx = nx; ldy = ny; ldz = nz ! No augmentation, the caller does not support it.
      27              :    call dfti_fftug(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, nspinor*ndat, istwf_k, mgfft, &
      28       215006 :                    kg_k, gbound_k, ug, ur)
      29              : 
      30              :  case (FFT_SG, FFT_SG2002)
      31           41 :    ldx=ngfft(4); ldy=ngfft(5); ldz=ngfft(6) ! Here augmentation is supported.
      32          123 :    ABI_MALLOC(fofgin,(2,npw_k))
      33          205 :    ABI_MALLOC(fofr,(2,ldx,ldy,ldz))
      34              : 
      35           41 :    if (istwf_k<=2) then
      36              : 
      37           38 :      do idat=1,nspinor*ndat
      38           25 :        gspad = (idat-1)*npw_k
      39           25 :        rspad = (idat-1)*nfft
      40              :        ! Have to convert from CPLX to REAL.
      41              :        ! Cannot use ZCOPY or C tricks here because this is a template used
      42              :        ! both for single and double precision.
      43      1111805 :        do ig=1,npw_k
      44      1111780 :          fofgin(1,ig) = DBLE (ug(ig+gspad))
      45      1111805 :          fofgin(2,ig) = AIMAG(ug(ig+gspad))
      46              :        end do
      47              : 
      48              :        !this is not thread-safe! use ndat1
      49              :        call sg_fftrisc(cplex0,dum_denpot,fofgin,dum_fofgout,fofr,gbound_k,gbound_k,istwf_k,&
      50           25 :                        kg_k,dum_kg_kout,mgfft,ndat1,ngfft,npw_k,0,ldx,ldy,ldz,option0,weight1_r,weight1_i)
      51              : 
      52              :        ! Fill the output array on the ngfft(1:3) mesh.
      53         2488 :        do iz=1,nz
      54       242975 :          do iy=1,ny
      55     23887950 :            do ix=1,nx
      56     23645000 :              ifft = ix + (iy-1)*nx + (iz-1)*nx*ny + rspad
      57     23885500 :              ur(ifft) = DCMPLX(fofr(1,ix,iy,iz), fofr(2,ix,iy,iz))
      58              :            end do
      59              :          end do
      60              :        end do
      61              :        !
      62              :      end do ! idat
      63              : 
      64              :    else
      65              :      !
      66              :      ! sg_fftrisc does not accept istwf_k>2 when option=0.
      67              :      ! in this case nspinor is always 1.
      68              :      !
      69          112 :      ABI_MALLOC(ftarr,(2,ldx,ldy,ldz))
      70           98 :      do idat=1,nspinor*ndat
      71           70 :        gspad = (idat-1)*npw_k
      72           70 :        rspad = (idat-1)*nfft
      73              : 
      74      2199120 :        do ig=1,npw_k ! Have to convert from CPLX to REAL.
      75      2199050 :          fofgin(1,ig) = DBLE (ug(ig+gspad))
      76      2199120 :          fofgin(2,ig) = AIMAG(ug(ig+gspad))
      77              :        end do
      78              :        !
      79              :        ! Reconstruct the full G-sphere.
      80           70 :        call sphere(fofgin,ndat1,npw_k,fofr,nx,ny,nz,ldx,ldy,ldz,kg_k,istwf_k,iflag1,me_g0,shiftg,symm,xnorm1)
      81              : 
      82              :        ! ndat=1 here.
      83           70 :        call sg_fftpad(fftcache,mgfft,nx,ny,nz,ldx,ldy,ldz,1,gbound_k,+1,fofr,ftarr)
      84              : 
      85         7098 :        do iz=1,nz  ! Fill the output array on the ngfft(1:3) mesh.
      86       707070 :          do iy=1,ny
      87     70707000 :            do ix=1,nx
      88     70000000 :              ifft = ix + (iy-1)*nx + (iz-1)*nx*ny + rspad
      89     70700000 :              ur(ifft) = DCMPLX(ftarr(1,ix,iy,iz), ftarr(2,ix,iy,iz))
      90              :            end do
      91              :          end do
      92              :        end do
      93              :      end do
      94           28 :      ABI_FREE(ftarr)
      95              :    end if
      96              : 
      97           41 :    ABI_FREE(fofgin)
      98           41 :    ABI_FREE(fofr)
      99              : 
     100              :  case default
     101       215047 :    ABI_ERROR(sjoin("Wrong fftalga:", itoa(fftalga)))
     102              :  end select
     103              : 
     104              : end block