TABLE OF CONTENTS


ABINIT/m_fftcore [ Modules ]

[ Top ] [ Modules ]

NAME

  m_fftcore

FUNCTION

  Low-level tools for FFT (sequential and MPI parallel version)
  It also provides helper functions to set up the list of G vectors
  inside a sphere or to count them.

COPYRIGHT

  Copyright (C) 2014-2018 ABINIT group (SG, XG, AR, MG, MT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

TODO

  1) Pass distribfft instead of MPI_enreg to simplify the API and facilitate code-reuse.

  2) Merge this module with m_distribfft

  3) Get rid of paral_kgb and MPI_type! This is a low-level module that may be called by other
     code in which paral_kgb is meaningless! FFT tables and a MPI communicator are sufficient.

SOURCE

 29 #if defined HAVE_CONFIG_H
 30 #include "config.h"
 31 #endif
 32 
 33 #include "abi_common.h"
 34 
 35 module m_fftcore
 36 
 37  use defs_basis
 38  use m_abicore
 39  use m_errors
 40  use m_xmpi
 41  use m_sort
 42 
 43  use m_time,         only : timab
 44  use defs_abitypes,  only : MPI_type
 45  use m_mpinfo,       only : destroy_mpi_enreg, initmpi_seq
 46 
 47  implicit none
 48 
 49  private
 50 
 51  public :: fftalg_isavailable  ! True if the FFT library specified by fftalg is available.
 52  public :: fftalg_has_mpi      ! True if fftalg provides MPI-FFTs.
 53  public :: fftalg_for_npfft    ! Returns the default value for fftalg given the number of processors for the FFT.
 54  public :: fftalg_info         ! Returns strings with info on the FFT library specified by fftalg.
 55  public :: get_cache_kb        ! Returns the cache size in Kbs (based on CPP variables).
 56  public :: ngfft_seq           ! initialize ngfft(18) from the FFT divisions (assume sequential FFT)
 57  public :: print_ngfft         ! Print the content of ngfft(18) in explicative format.
 58  public :: bound               ! Find distance**2 to boundary point of fft box nearest to kpt
 59  public :: getng               ! From ecut and metric tensor in reciprocal space, computes recommended ngfft(1:3)
 60  public :: sphereboundary      ! Finds the boundary of the basis sphere of G vectors
 61  public :: sphere
 62  public :: sphere_fft      ! Insert cg inside box.
 63  public :: sphere_fft1     ! TODO: This one should be replaced by sphere_fft.
 64  public :: change_istwfk   ! Change the istwfk mode of a set of wavefunctions (sequential version, same k-point)
 65  public :: kpgsph          ! Set up the G vector list
 66  public :: kpgcount        ! Give the number of G vector in each direction
 67  public :: get_kg          ! Helper function to calculate the set of G-vectors at a given kpoint (no MPI FFT)
 68  public :: kgindex         ! Compute the index of each plane wave on a FFT grid.
 69 
 70  ! Low-level tools for MPI FFT
 71  public :: switch
 72  public :: switch_cent
 73  public :: switchreal
 74  public :: switchreal_cent
 75  public :: scramble
 76  public :: fill
 77  public :: fill_cent
 78  public :: unfill
 79  public :: unfill_cent
 80  public :: unmpiswitch
 81  public :: unswitch
 82  public :: unswitchreal_cent
 83  public :: unmpiswitch_cent
 84  public :: unscramble
 85  public :: unswitch_cent
 86  public :: unswitchreal
 87  public :: mpiswitch
 88  public :: mpiswitch_cent
 89 
 90  public :: mpifft_fg2dbox
 91  public :: mpifft_fg2dbox_dpc
 92  public :: mpifft_dbox2fg
 93  public :: mpifft_dbox2fg_dpc
 94  public :: mpifft_dbox2fr
 95  public :: mpifft_dbox2fr_dpc
 96  public :: mpifft_fr2dbox
 97  public :: mpifft_fr2dbox_dpc
 98  public :: mpifft_collect_datar           ! Collect a real-space MPI-FFT distributed array on each proc.
 99 
100  public :: indfftrisc
101 
102  public :: addrho
103  public :: multpot
104 
105 ! *************************************************************************
106 
107 !----------------------------------------------------------------------
108 ! Private variables
109 
110 #define FFTALGA_SIZE 5
111  character(len=*),private,parameter :: fftalga2name(1:FFTALGA_SIZE)= &
112 & (/"Goedecker     ", &
113 &   "Vendor FFT    ", &
114 &   "FFTW3         ", &
115 &   "Goedecker2002 ", &
116 &   "DFTI          " /)
117 
118 #define FFTALGB_SIZE 1
119  character(len=*),private,parameter :: fftalgb2name(0:FFTALGB_SIZE)= &
120 & (/"C2C",&
121 &   "R2C"/)
122 
123 #define FFTALGC_SIZE 2
124  character(len=*),private,parameter :: fftalgc2name(0:FFTALGC_SIZE)= &
125 & (/"No pad         ",&
126 &   "zero-pad       ",&
127 &   "zero-pad+cache "/)
128 
129 contains

m_fftcore/addrho [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  addrho

FUNCTION

   Add the contribution to the density generated by n1dfft x-y planes

INPUTS

  icplexwf=1 if u(r) is real, 2 otherwise.
  includelast
  nd1,nd2=Leading dimensions of rhopart.
  n2=FFT dimension along y.
  lot=2nd Leading dimension of zw (cache blocking factor).
  n1dfft=Number of 1D FFTs along y performed.
  zw(2,lot,n2)=Array with the x-y planes (wavefunction in real space).
  weigth=Weight factor for the density.

SIDE EFFECTS

   rhopart(nd1,nd2)=density in the x-y plane, accumulated in output.

PARENTS

SOURCE

4997 pure subroutine addrho(icplexwf,includelast,nd1,nd2,n2,lot,n1dfft,zw,rhopart,weight_r,weight_i)
4998 
4999 
5000 !This section has been created automatically by the script Abilint (TD).
5001 !Do not modify the following lines by hand.
5002 #undef ABI_FUNC
5003 #define ABI_FUNC 'addrho'
5004 !End of the abilint section
5005 
5006  implicit none
5007 
5008 !Arguments ------------------------------------
5009  integer,intent(in) :: icplexwf,includelast,nd1,nd2,n2,lot,n1dfft
5010  real(dp),intent(in) :: zw(2,lot,n2)
5011  real(dp),intent(inout) :: rhopart(nd1,nd2)
5012  real(dp),intent(in) :: weight_i,weight_r
5013 
5014 !Local variables-------------------------------
5015  integer :: i2,j
5016 
5017 ! *************************************************************************
5018 
5019  if (icplexwf==2) then
5020    ! Complex wavefunction in real space.
5021    do i2=1,n2-1,2
5022      do j=1,n1dfft
5023        rhopart(j,i2) =   rhopart(j,i2) +   weight_r*zw(1,j,i2)**2+weight_i*zw(2,j,i2)**2
5024        rhopart(j,i2+1) = rhopart(j,i2+1) + weight_r*zw(1,j,i2+1)**2+weight_i*zw(2,j,i2+1)**2
5025      end do
5026    end do
5027 
5028    if (2*(n2/2)/=n2) then
5029      do j=1,n1dfft
5030        rhopart(j,n2  )=rhopart(j,n2  )+weight_r*zw(1,j,n2  )**2+weight_i*zw(2,j,n2  )**2
5031      end do
5032    end if
5033  else
5034    ! The wavefunction is real, in real space
5035    if (includelast==1) then
5036      do i2=1,n2
5037        do j=1,n1dfft
5038          rhopart(2*j-1,i2)=rhopart(2*j-1,i2)+zw(1,j,i2)**2*weight_r
5039          rhopart(2*j  ,i2)=rhopart(2*j  ,i2)+zw(2,j,i2)**2*weight_i
5040        end do
5041      end do
5042    else
5043      do i2=1,n2
5044        do j=1,n1dfft-1
5045          rhopart(2*j-1,i2)=rhopart(2*j-1,i2)+zw(1,j,i2)**2*weight_r
5046          rhopart(2*j  ,i2)=rhopart(2*j  ,i2)+zw(2,j,i2)**2*weight_i
5047        end do
5048        rhopart(2*n1dfft-1,i2)=rhopart(2*n1dfft-1,i2)+zw(1,n1dfft,i2)**2*weight_r
5049      end do
5050    end if
5051 
5052  end if
5053 
5054 end subroutine addrho

m_fftcore/bound [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 bound

FUNCTION

 For given kpt, ngfft, and gmet,
  Find distance**2 to boundary point of fft box nearest to kpt
  Find distance**2 to boundary point of fft box farthest to kpt

INPUTS

  kpt(3)=real input k vector (reduced coordinates)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  gmet(3,3)=reciprocal space metric (currently in Bohr**-2)

OUTPUT

  dsqmax=maximum distance**2 from k to boundary in Bohr**-2.
  dsqmin=minimum distance**2 from k to boundary in Bohr**-2.
  gbound(3)=coords of G on boundary (correspnding to gsqmin)
  plane=which plane min occurs in (1,2, or 3 for G1,etc).

NOTES

 Potential trouble: this routine was written assuming kpt lies inside
 first Brillouin zone.  No measure is taken to fold input kpt back
 into first zone.  Given arbitrary kpt, this will cause trouble.

PARENTS

      getcut,getng

CHILDREN

SOURCE

582 subroutine bound(dsqmax,dsqmin,gbound,gmet,kpt,ngfft,plane)
583 
584 
585 !This section has been created automatically by the script Abilint (TD).
586 !Do not modify the following lines by hand.
587 #undef ABI_FUNC
588 #define ABI_FUNC 'bound'
589 !End of the abilint section
590 
591  implicit none
592 
593 !Arguments ------------------------------------
594 !scalars
595  integer,intent(out) :: plane
596  real(dp),intent(out) :: dsqmax,dsqmin
597 !arrays
598  integer,intent(in) :: ngfft(18)
599  integer,intent(out) :: gbound(3)
600  real(dp),intent(in) :: gmet(3,3),kpt(3)
601 
602 !Local variables-------------------------------
603 !scalars
604  integer :: i1,i1min,i2,i2min,i3,i3min
605  real(dp) :: dsm,dsp
606  character(len=500) :: message
607 
608 ! *************************************************************************
609 
610 ! dsq(i1,i2,i3)=gmet(1,1)*(kpt(1)+dble(i1))**2&
611 !& +gmet(2,2)*(kpt(2)+dble(i2))**2&
612 !& +gmet(3,3)*(kpt(3)+dble(i3))**2&
613 !& +2._dp*(gmet(1,2)*(kpt(1)+dble(i1))*(kpt(2)+dble(i2))&
614 !& +gmet(2,3)*(kpt(2)+dble(i2))*(kpt(3)+dble(i3))&
615 !& +gmet(3,1)*(kpt(3)+dble(i3))*(kpt(1)+dble(i1)))
616 
617 !Set plane to impossible value
618  plane=0
619 
620 !look at +/- g1 planes:
621  dsqmax=zero
622  dsqmin=dsq(ngfft(1)/2,-ngfft(2)/2,-ngfft(3)/2,gmet,kpt)+0.01_dp
623  do i2=-ngfft(2)/2,ngfft(2)/2
624    do i3=-ngfft(3)/2,ngfft(3)/2
625      dsp = dsq(ngfft(1)/2, i2, i3,gmet,kpt)
626      dsm = dsq( - ngfft(1)/2, i2, i3,gmet,kpt)
627      if (dsp>dsqmax) dsqmax = dsp
628      if (dsm>dsqmax) dsqmax = dsm
629      if (dsp<dsqmin) then
630        dsqmin = dsp
631        i1min = ngfft(1)/2
632        i2min = i2
633        i3min = i3
634        plane=1
635      end if
636      if (dsm<dsqmin) then
637        dsqmin = dsm
638        i1min =  - ngfft(1)/2
639        i2min = i2
640        i3min = i3
641        plane=1
642      end if
643    end do
644  end do
645 !
646 !+/- g2 planes:
647  do i1=-ngfft(1)/2,ngfft(1)/2
648    do i3=-ngfft(3)/2,ngfft(3)/2
649      dsp = dsq(i1,ngfft(2)/2,i3,gmet,kpt)
650      dsm = dsq(i1,-ngfft(2)/2,i3,gmet,kpt)
651      if (dsp>dsqmax) dsqmax = dsp
652      if (dsm>dsqmax) dsqmax = dsm
653      if (dsp<dsqmin) then
654        dsqmin = dsp
655        i1min = i1
656        i2min = ngfft(2)/2
657        i3min = i3
658        plane=2
659      end if
660      if (dsm<dsqmin) then
661        dsqmin = dsm
662        i1min = i1
663        i2min =  - ngfft(2)/2
664        i3min = i3
665        plane=2
666      end if
667    end do
668  end do
669 !
670 !+/- g3 planes:
671  do i1=-ngfft(1)/2,ngfft(1)/2
672    do i2=-ngfft(2)/2,ngfft(2)/2
673      dsp = dsq(i1,i2,ngfft(3)/2,gmet,kpt)
674      dsm = dsq(i1,i2,-ngfft(3)/2,gmet,kpt)
675      if (dsp>dsqmax) dsqmax = dsp
676      if (dsm>dsqmax) dsqmax = dsm
677      if (dsp<dsqmin) then
678        dsqmin = dsp
679        i1min = i1
680        i2min = i2
681        i3min = ngfft(3)/2
682        plane=3
683      end if
684      if (dsm<dsqmin) then
685        dsqmin = dsm
686        i1min = i1
687        i2min = i2
688        i3min =  - ngfft(3)/2
689        plane=3
690      end if
691    end do
692  end do
693 
694  if (plane==0) then
695 !  Trouble: missed boundary somehow
696    write(message, '(a,a,a,3f9.4,a,3i5,a,a,a,a,a)' )&
697 &   'Trouble finding boundary of G sphere for',ch10,&
698 &   'kpt=',kpt(:),' and ng=',ngfft(1:3),ch10,&
699 &   'Action : check that kpt lies',&
700 &   'reasonably within first Brillouin zone; ',ch10,&
701 &   'else code bug, contact ABINIT group.'
702    MSG_BUG(message)
703  end if
704 
705  gbound(1)=i1min
706  gbound(2)=i2min
707  gbound(3)=i3min
708 
709  contains
710 
711    function dsq(i1,i2,i3,gmet,kpt)
712 
713 
714 !This section has been created automatically by the script Abilint (TD).
715 !Do not modify the following lines by hand.
716 #undef ABI_FUNC
717 #define ABI_FUNC 'dsq'
718 !End of the abilint section
719 
720      integer :: i1,i2,i3
721      real(dp) :: dsq
722      real(dp) :: kpt(3),gmet(3,3)
723 
724      dsq=gmet(1,1)*(kpt(1)+dble(i1))**2&
725 &      +gmet(2,2)*(kpt(2)+dble(i2))**2&
726 &      +gmet(3,3)*(kpt(3)+dble(i3))**2&
727 &      +2._dp*(gmet(1,2)*(kpt(1)+dble(i1))*(kpt(2)+dble(i2))&
728 &      +gmet(2,3)*(kpt(2)+dble(i2))*(kpt(3)+dble(i3))&
729 &      +gmet(3,1)*(kpt(3)+dble(i3))*(kpt(1)+dble(i1)))
730    end function dsq
731 
732 end subroutine bound

m_fftcore/change_istwfk [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 change_istwfk

FUNCTION

 This function allows one to change the time-reversal storage mode (istwfk)
 of a *full* set of u(G). It does not support MPI-FFT!

INPUTS

 from_npw=number of G vectors in input from_cg
 from_kg_k(3,npw)=integer coordinates of the G vectors of from_cg
 from_istwfk=option parameter that describes the storage in from_cg
 to_npw=number of G vectors in output to_cg
 to_kg_k(3,npw)=integer coordinates of the G vectors in to_cg
 to_istwfk=option parameter that describes the storage in to_cg
 n1,n2,n3=physical dimension of the box (must be large enough to contain the sphere, no check is done)
 ndat=number of wavefunctions
 from_cg(2,from_npw*ndat)= Input u(g) values

 OUTPUTS
 to_cg(2,to_npw*ndat)= Output u(g) defined on the list of vectors to_kg_k with time-reversal mode to_istwfk

PARENTS

      m_fft

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

2127 subroutine change_istwfk(from_npw,from_kg,from_istwfk,to_npw,to_kg,to_istwfk,n1,n2,n3,ndat,from_cg,to_cg)
2128 
2129 
2130 !This section has been created automatically by the script Abilint (TD).
2131 !Do not modify the following lines by hand.
2132 #undef ABI_FUNC
2133 #define ABI_FUNC 'change_istwfk'
2134 !End of the abilint section
2135 
2136  implicit none
2137 
2138 !Arguments ------------------------------------
2139 !scalars
2140  integer,intent(in) :: from_npw,from_istwfk,to_npw,to_istwfk,n1,n2,n3,ndat
2141 !arrays
2142  integer,intent(in) :: from_kg(3,from_npw),to_kg(3,to_npw)
2143  real(dp),intent(inout) :: from_cg(2,from_npw*ndat) ! out due to sphere!
2144  real(dp),intent(inout) :: to_cg(2,to_npw*ndat)
2145 
2146 !Local variables-------------------------------
2147 !scalars
2148  integer :: n4,n5,n6
2149  real(dp),parameter :: xnorm1=one
2150 !arrays
2151  integer,parameter :: shiftg0(3)=0,me_g0=1
2152  integer,parameter :: symmE(3,3)=reshape([1,0,0,0,1,0,0,0,1],[3,3])
2153  real(dp),allocatable :: cfft(:,:,:,:)
2154 
2155 ! *************************************************************************
2156 
2157  n4=2*(n1/2)+1
2158  n5=2*(n2/2)+1
2159  n6=2*(n3/2)+1
2160 
2161  ABI_MALLOC(cfft, (2,n4,n5,n6*ndat))
2162 
2163  ! iflag=1 ==> insert from_cg into cfft.
2164  call sphere(from_cg,ndat,from_npw,cfft,n1,n2,n3,n4,n5,n6,from_kg,from_istwfk,+1,me_g0,shiftg0,symmE,xnorm1)
2165 
2166  ! iflag=-1 ==> extract to_cg from cfft.
2167  call sphere(to_cg,ndat,to_npw,cfft,n1,n2,n3,n4,n5,n6,to_kg,to_istwfk,-1,me_g0,shiftg0,symmE,xnorm1)
2168 
2169  ABI_FREE(cfft)
2170 
2171 end subroutine change_istwfk

m_fftcore/fftalg_for_npfft [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 fftalg_for_npfft

FUNCTION

  Returns the default value of fftalg given the number of MPI nodes
  to be used in the FFTs.

INPUTS

  nproc_fft=Number of processors used for MPI FFT

OUTPUT

  fftalg=Integer used to select the FFT library.

PARENTS

SOURCE

255 pure function fftalg_for_npfft(nproc_fft) result(fftalg)
256 
257 
258 !This section has been created automatically by the script Abilint (TD).
259 !Do not modify the following lines by hand.
260 #undef ABI_FUNC
261 #define ABI_FUNC 'fftalg_for_npfft'
262 !End of the abilint section
263 
264  implicit none
265 
266 !Arguments ------------------------------------
267 !scalars
268  integer,intent(in) :: nproc_fft
269  integer :: fftalg
270 
271 ! *************************************************************************
272 
273  ! Default  for the sequential case.
274  fftalg = 112
275 
276  ! Use Goedecker2002 if fftalg does not support MPI (e.g 112)
277  if (nproc_fft > 1) fftalg = 401
278  !if (nproc_fft > 1) fftalg = 402
279 
280 #ifdef HAVE_FFT_FFTW3
281  fftalg = 312
282 #elif defined HAVE_FFT_DFTI
283  fftalg = 512
284  if (nproc_fft > 1) fftalg = 401  ! MPI-FFT with DFTI is not implemented yet
285 #endif
286 
287  !if (nproc_fft > 1) fftalg = 401 ! This is to revert to the old behavior.
288 
289 end function fftalg_for_npfft

m_fftcore/fftalg_has_mpi [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 fftalg_has_mpi

FUNCTION

  True if the FFT library specified by fftalg is available.

INPUTS

  fftalg=Input variable.

PARENTS

SOURCE

203 pure function fftalg_has_mpi(fftalg) result(ans)
204 
205 
206 !This section has been created automatically by the script Abilint (TD).
207 !Do not modify the following lines by hand.
208 #undef ABI_FUNC
209 #define ABI_FUNC 'fftalg_has_mpi'
210 !End of the abilint section
211 
212  implicit none
213 
214 !Arguments ------------------------------------
215 !scalars
216  integer,intent(in) :: fftalg
217  logical :: ans
218 
219 !Local variables-------------------------------
220 !scalars
221  integer :: fftalga,fftalgb,fftalgc
222 
223 ! *************************************************************************
224 
225  ans = .False.
226  fftalga = fftalg/100; fftalgb = mod(fftalg,100)/10; fftalgc = mod(fftalg,10)
227 
228  if (fftalga == FFT_FFTW3) ans = .True.
229  !if (fftalga == FFT_DFTI)  ans = .True.
230  if (fftalga == FFT_SG2002) ans = .True.
231 
232 end function fftalg_has_mpi

m_fftcore/fftalg_info [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 fftalg_info

FUNCTION

  Returns info on the FFT library specified by fftalg (ngfft(7))

INPUTS

  fftalg=Input variable.

OUTPUT

  library=String with the name of FFT library
  cplex_mode= String defining whether the FFT library supports real<-->complex transforms.
  padding_mode=Padding mode.

PARENTS

      m_fft,m_fft_prof

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

317 subroutine fftalg_info(fftalg,library,cplex_mode,padding_mode)
318 
319 
320 !This section has been created automatically by the script Abilint (TD).
321 !Do not modify the following lines by hand.
322 #undef ABI_FUNC
323 #define ABI_FUNC 'fftalg_info'
324 !End of the abilint section
325 
326  implicit none
327 
328 !Arguments ------------------------------------
329 !scalars
330  integer,intent(in) :: fftalg
331  character(len=*),intent(out) :: library,cplex_mode,padding_mode
332 
333 !Local variables-------------------------------
334 !scalars
335  integer :: fftalga,fftalgb,fftalgc
336 
337 ! *************************************************************************
338 
339  library = "Unknown"; cplex_mode = "Unknown"; padding_mode  = "Unknown"
340 
341  fftalga=fftalg/100
342  if (fftalga>0 .and. fftalga<=FFTALGA_SIZE) library = fftalga2name(fftalga)
343 
344  fftalgb=mod(fftalg,100)/10
345  if (fftalgb>=0 .and. fftalgb<=FFTALGB_SIZE) cplex_mode = fftalgb2name(fftalgb)
346 
347  fftalgc=mod(fftalg,10)
348  if (fftalgc>=0 .and. fftalgc<=FFTALGC_SIZE) padding_mode = fftalgc2name(fftalgc)
349 
350 end subroutine fftalg_info

m_fftcore/fftalg_isavailable [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 fftalg_isavailable

FUNCTION

  Returns TRUE if the FFT library specified by fftalg (ngfft(7)) is available

INPUTS

  fftalg=Input variable.

PARENTS

SOURCE

148 pure function fftalg_isavailable(fftalg) result(ans)
149 
150 
151 !This section has been created automatically by the script Abilint (TD).
152 !Do not modify the following lines by hand.
153 #undef ABI_FUNC
154 #define ABI_FUNC 'fftalg_isavailable'
155 !End of the abilint section
156 
157  implicit none
158 
159 !Arguments ------------------------------------
160 !scalars
161  integer,intent(in) :: fftalg
162  logical :: ans
163 
164 !Local variables-------------------------------
165 !scalars
166  integer :: fftalga,fftalgb,fftalgc
167 
168 ! *************************************************************************
169 
170  ans = .TRUE.
171  fftalga = fftalg/100
172  fftalgb = mod(fftalg,100)/10
173  fftalgc = mod(fftalg,10)
174 
175  ! Optional FFT libraries.
176 #ifndef HAVE_FFT_FFTW3
177  if (fftalga == FFT_FFTW3) ans = .FALSE.
178 #endif
179 
180 #ifndef HAVE_FFT_DFTI
181  if (fftalga == FFT_DFTI) ans = .FALSE.
182 #endif
183 
184 end function fftalg_isavailable

m_fftcore/fill [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  fill

FUNCTION

   Receives a set of z-lines in reciprocal space,
   insert the values in the cache work array zw (no padding)

INPUTS

  nd1,nd3=Dimensions of the input array zf.
  lot=Cache blocking factor.
  n1dfft=Number of 1D FFTs to perform
  n3=Dimension of the transform along z
  zf(2,nd1,nd3)=Input array

OUTPUT

  zw(2,lot,n3)=Cache work array with the z-lines.

PARENTS

SOURCE

2593 pure subroutine fill(nd1,nd3,lot,n1dfft,n3,zf,zw)
2594 
2595 
2596 !This section has been created automatically by the script Abilint (TD).
2597 !Do not modify the following lines by hand.
2598 #undef ABI_FUNC
2599 #define ABI_FUNC 'fill'
2600 !End of the abilint section
2601 
2602  implicit none
2603 
2604 !Arguments ------------------------------------
2605  integer,intent(in) :: nd1,nd3,lot,n1dfft,n3
2606  real(dp),intent(in) :: zf(2,nd1,nd3)
2607  real(dp),intent(inout) :: zw(2,lot,n3)
2608 
2609 ! local variables
2610  integer :: i1,i3
2611 
2612 ! *************************************************************************
2613 
2614  do i3=1,n3
2615    do i1=1,n1dfft
2616      zw(1,i1,i3)=zf(1,i1,i3)
2617      zw(2,i1,i3)=zf(2,i1,i3)
2618    end do
2619  end do
2620 
2621 end subroutine fill

m_fftcore/fill_cent [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 fill_cent

FUNCTION

   Receives a set of z-lines in reciprocal space,
   insert the values in cache work array defined on the FFT box
   and pads the central of the frequency region with zeros.

INPUTS

  md1,md3=Leading dimension of zf along x and z
  lot=second dimension of zw (cache blocking factor)
  n1dfft=Number of 1d transforms to be performed along z.
  max3=Max G_z in the small box enclosing the G-sphere.
  m3=Number of points in the *small* box enclosing the G-sphere
  n3=Dimension of the FFT transform along z
  zf(2,md1,md3)=x-z planes in reciprocal space

OUTPUT

   zw(2,lot,n3)= Filled cache work array.
     zw(:,1:n1dfft,n3) contains the lines to be transformed along.

PARENTS

SOURCE

2652 pure subroutine fill_cent(md1,md3,lot,n1dfft,max3,m3,n3,zf,zw)
2653 
2654 
2655 !This section has been created automatically by the script Abilint (TD).
2656 !Do not modify the following lines by hand.
2657 #undef ABI_FUNC
2658 #define ABI_FUNC 'fill_cent'
2659 !End of the abilint section
2660 
2661  implicit none
2662 
2663 !Arguments ------------------------------------
2664  integer,intent(in) :: md1,md3,lot,n1dfft,max3,m3,n3
2665  real(dp),intent(in) :: zf(2,md1,md3)
2666  real(dp),intent(inout) :: zw(2,lot,n3)
2667 
2668 !Local variables-------------------------------
2669 !scalars
2670  integer :: i1,i3
2671 
2672 ! *************************************************************************
2673 
2674  ! Here, zero and positive frequencies
2675  do i3=1,max3+1
2676    do i1=1,n1dfft
2677      zw(1,i1,i3)=zf(1,i1,i3)
2678      zw(2,i1,i3)=zf(2,i1,i3)
2679    end do
2680  end do
2681 
2682  ! Fill the center region with zeros
2683  do i3=max3+2,n3-m3+max3+1
2684    do i1=1,n1dfft
2685      zw(1,i1,i3)=zero
2686      zw(2,i1,i3)=zero
2687    end do
2688  end do
2689 
2690  ! Here, negative frequencies
2691  do i3=max3+2,m3
2692    do i1=1,n1dfft
2693      zw(1,i1,i3+n3-m3)=zf(1,i1,i3)
2694      zw(2,i1,i3+n3-m3)=zf(2,i1,i3)
2695    end do
2696  end do
2697 
2698 end subroutine fill_cent

m_fftcore/get_cache_kb [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 get_cache_kb

FUNCTION

  Returns the cache size in KB to be used for cache blocking algorithms in the FFT routines.
  The value is derived from the values of the CPP options defined in config.h

TODO

   Use C to get the real cache size.
   See http://stackoverflow.com/questions/12594208/c-program-to-determine-levels-size-of-cache

PARENTS

SOURCE

371 pure function get_cache_kb()
372 
373 
374 !This section has been created automatically by the script Abilint (TD).
375 !Do not modify the following lines by hand.
376 #undef ABI_FUNC
377 #define ABI_FUNC 'get_cache_kb'
378 !End of the abilint section
379 
380  implicit none
381 
382 !Local variables-------------------------------
383 !scalars
384  integer :: get_cache_kb
385 
386 ! *************************************************************************
387 
388  ! Default value
389  get_cache_kb = 16
390  !get_cache_kb = 32
391  !get_cache_kb = 256
392 
393 end function get_cache_kb

m_fftcore/get_kg [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  get_kg

FUNCTION

  Helper function to calculate the set of G-vectors at a given kpoint.
  without taking advantage of FFT parallelism and G-vector distributions.

INPUTS

  kpoint(3)=The k-point in reduced coordinates.
  ecut=Cutoff energy for planewave basis set.
  gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
  istwfk=Options defining if time-reversal is used to decrease the number of G"s.

OUTPUT

  npw_k=Total number of G-vectors in the full G-sphere.

SIDE EFFECTS

  kg_k(:,:):
   allocated inside the routine.
   output: kg_k(3,npw_k) contains the list of G-vectors.

PARENTS

      fftprof,m_ebands,m_fft,m_fft_prof,m_gkk,m_io_kss,m_phgamma,m_phpi
      m_shirley,m_sigmaph,m_wfd,m_wfk,outkss

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

4828 subroutine get_kg(kpoint,istwf_k,ecut,gmet,npw_k,kg_k)
4829 
4830 
4831 !This section has been created automatically by the script Abilint (TD).
4832 !Do not modify the following lines by hand.
4833 #undef ABI_FUNC
4834 #define ABI_FUNC 'get_kg'
4835 !End of the abilint section
4836 
4837  implicit none
4838 
4839 !Arguments ------------------------------------
4840 !scalars
4841  integer,intent(in) :: istwf_k
4842  integer,intent(out) :: npw_k
4843  real(dp),intent(in) :: ecut
4844 !arrays
4845  integer,allocatable,intent(out) :: kg_k(:,:)
4846  real(dp),intent(in) :: gmet(3,3),kpoint(3)
4847 
4848 !Local variables-------------------------------
4849 !scalars
4850  integer,parameter :: mkmem_=1,exchn2n3d0=0,ikg0=0
4851  integer :: npw_k_test
4852  type(MPI_type) :: MPI_enreg_seq
4853 !arrays
4854  integer :: kg_dum(3,0)
4855 
4856 ! *********************************************************************
4857 
4858  call initmpi_seq(MPI_enreg_seq)
4859  !
4860  ! * Calculate the number of G-vectors for this k-point.
4861  call kpgsph(ecut,exchn2n3d0,gmet,ikg0,0,istwf_k,kg_dum,kpoint,0,MPI_enreg_seq,0,npw_k)
4862  !
4863  ! * Allocate and calculate the set of G-vectors.
4864  ABI_MALLOC(kg_k,(3,npw_k))
4865  call kpgsph(ecut,exchn2n3d0,gmet,ikg0,0,istwf_k,kg_k,kpoint,mkmem_,MPI_enreg_seq,npw_k,npw_k_test)
4866 
4867  call destroy_mpi_enreg(MPI_enreg_seq)
4868 
4869 end subroutine get_kg

m_fftcore/getng [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 getng

FUNCTION

 From ecut and metric tensor in reciprocal space, computes recommended ngfft(1:3)
 Also computes the recommended value of nfft and mgfft
 Pay attention that the FFT grid must be compatible with the symmetry operations (see irrzg.f).

INPUTS

 boxcutmin=minimum value of boxcut admitted (boxcut is the ratio
  between the radius of the sphere contained in the FFT box, and the
  radius of the planewave sphere) : usually 2.0 .
 ecut=energy cutoff in Hartrees
 gmet(3,3)=reciprocal space metric (bohr**-2).
 kpt(3)=input k vector in terms of reciprocal lattice primitive translations
 me_fft=index of the processor in the FFT set (use 0 if sequential)
 nproc_fft=number of processors in the FFT set (use 1 if sequential)
 nsym=number of symmetry elements in group
 paral_fft=0 if no FFT parallelisation ; 1 if FFT parallelisation
 symrel(3,3,nsym)=symmetry matrices in real space (integers)

OUTPUT

 mgfft= max(ngfft(1),ngfft(2),ngfft(3))
 nfft=number of points in the FFT box=ngfft(1)*ngfft(2)*ngfft(3)/nproc_fft

SIDE EFFECTS

 Input/Output
 ngfft(1:18)=integer array with FFT box dimensions and other
   information on FFTs. On input ngfft(1:3) contains
   optional trial values. If ngfft(1:3)/minbox is greater than value
   calculated to avoid wrap-around error and ngfft obeys constraint
   placed by the FFT routine that is used
   then ngfft(1:3) is left unchanged. Otherwise set to value computed in now.

 Note that there is the possibility of an undetected error if we
 are dealing with a cubic unit cell and ngfft(1), ngfft(2) and ngfft(3)
 are different. In the future we should handle this case.

 ngfft(4),ngfft(5),ngfft(6)= modified values to avoid cache trashing,
        presently : ngfft(4)=ngfft(1)+1 if ngfft(1) is even ;
                    ngfft(5)=ngfft(2)+1 if ngfft(2) is even.
           in the other cases, ngfft(4:6)=ngfft(1:3).
   Other choices may better, but this is left for the future.
 ngfft(7)=choice for FFT algorithm, see the input variable fftalg
 ngfft(8)=size of the cache, in bytes (not used here presently).!!
   other ngfft slots are used for parallelism see ~abinit/doc/variables/vargs.htm#ngfft
 [unit] = Output Unit number (DEFAULT std_out)

PARENTS

      fftprof,m_fft,m_fft_prof,memory_eval,mpi_setup,mrgscr,scfcv

CHILDREN

      bound,print_ngfft,sort_int,wrtout

SOURCE

 792 subroutine getng(boxcutmin,ecut,gmet,kpt,me_fft,mgfft,nfft,ngfft,nproc_fft,nsym,paral_fft,symrel,&
 793 &                ngfftc,use_gpu_cuda,unit) ! optional
 794 
 795  use defs_fftdata,  only : mg
 796 
 797 !This section has been created automatically by the script Abilint (TD).
 798 !Do not modify the following lines by hand.
 799 #undef ABI_FUNC
 800 #define ABI_FUNC 'getng'
 801 !End of the abilint section
 802 
 803  implicit none
 804 
 805 !Arguments ------------------------------------
 806 !scalars
 807  integer,intent(in) :: me_fft,nproc_fft,nsym,paral_fft
 808  integer,intent(out) :: mgfft,nfft
 809  integer,optional,intent(in) :: unit,use_gpu_cuda
 810  real(dp),intent(in) :: boxcutmin,ecut
 811 !arrays
 812  integer,intent(in) :: symrel(3,3,nsym)
 813  integer,intent(in),optional :: ngfftc(3)
 814  integer,intent(inout) :: ngfft(18)
 815  real(dp),intent(in) :: gmet(3,3),kpt(3)
 816 
 817 !Local variables-------------------------------
 818 !scalars
 819  integer,save :: first=1,msrch(3),previous_paral_mode=0
 820  integer :: element,ii,index,isrch,isrch1,isrch2,isrch3,isym,jj,mu,paral_fft_
 821  integer :: plane,testok,tobechecked,ount,fftalga
 822  real(dp),parameter :: minbox=0.75_dp
 823  real(dp) :: dsqmax,dsqmin,ecutmx,prodcurrent,prodtrial,xx,yy
 824  logical :: testdiv
 825  character(len=500) :: message
 826  integer,parameter :: largest_ngfft=mg ! Goedecker FFT: any powers of 2, 3, and 5 - must be coherent with defs_fftdata.F90
 827  integer,parameter :: maxpow2 =16      ! int(log(largest_ngfft+half)/log(two))
 828  integer,parameter :: maxpow3 =6       ! int(log(largest_ngfft+half)/log(three))
 829  integer,parameter :: maxpow5 =6       ! int(log(largest_ngfft+half)/log(five))
 830  integer,parameter :: maxpow7 =0
 831  integer,parameter :: maxpow11=0
 832  integer,parameter :: mmsrch=(maxpow2+1)*(maxpow3+1)*(maxpow5+1)*(maxpow7+1)*(maxpow11+1)
 833  integer,save :: iperm(mmsrch),srch(mmsrch,3)
 834  integer(i8b) :: li_srch(mmsrch)
 835  integer :: divisor(3,3),gbound(3),imax(3),imin(3),ngcurrent(3)
 836  integer :: ngmax(3),ngsav(3),ngtrial(3)
 837 
 838 ! *************************************************************************
 839 
 840  ount = std_out; if (PRESENT(unit)) ount = unit
 841 
 842  fftalga = ngfft(7)/100
 843 
 844 !If not yet done, compute recommended (boxcut>=2) fft grid dimensions
 845 !In case we switch for paral to sequential mode, recompute srch.
 846 !This is the case e.g. when computing ngfftdiel in sequential mode
 847 !after an initial computation of ngfft in parallel
 848 
 849  paral_fft_=paral_fft;if (nproc_fft==0) paral_fft_=0
 850 
 851  if(first==1.or.paral_fft_ /= previous_paral_mode) then
 852    first=0; previous_paral_mode=paral_fft_
 853    srch(:,:)=0
 854 
 855    ! Factors of 2
 856    srch(1,1)=1
 857    do ii=1,maxpow2
 858      srch(ii+1,1)=srch(ii,1)*2
 859    end do
 860 
 861    ! Factors of 3
 862    index=maxpow2+1
 863    if(maxpow3>0)then
 864      do ii=1,max(1,maxpow3)
 865        srch(1+ii*index:(ii+1)*index,1)=3*srch(1+(ii-1)*index:ii*index,1)
 866      end do
 867    end if
 868 
 869    ! Factors of 5
 870    index=(maxpow3+1)*index
 871    if(maxpow5>0)then
 872      do ii=1,max(1,maxpow5)
 873        li_srch = 0
 874        li_srch(1+ii*index:(ii+1)*index)=5_i8b*srch(1+(ii-1)*index:ii*index,1)
 875        where (li_srch > huge(maxpow3)) li_srch = huge(maxpow3)
 876        srch(1+ii*index:(ii+1)*index,1)=li_srch(1+ii*index:(ii+1)*index)
 877      end do
 878    end if
 879 
 880    ! Factors of 7
 881    index=(maxpow5+1)*index
 882    if(maxpow7>0)then
 883      do ii=1,max(1,maxpow7)
 884        srch(1+ii*index:(ii+1)*index,1)=7*srch(1+(ii-1)*index:ii*index,1)
 885      end do
 886    end if
 887 
 888    ! Factors of 11
 889    index=(maxpow7+1)*index
 890    if(maxpow11>0)then
 891      do ii=1,max(1,maxpow11)
 892        srch(1+ii*index:(ii+1)*index,1)=11*srch(1+(ii-1)*index:ii*index,1)
 893      end do
 894    end if
 895 
 896    call sort_int(mmsrch,srch(:,1),iperm)
 897 
 898    do ii=1,mmsrch
 899      if(srch(ii,1)>largest_ngfft)exit
 900    end do
 901    msrch(1)=ii-1
 902 
 903    ! In case of FFT parallelism, one need ngfft(2) and ngfft(3) to be multiple of nproc_fft
 904    if(paral_fft_==1)then
 905      msrch(2)=0
 906      do ii=1,msrch(1)
 907        if(modulo(srch(ii,1),nproc_fft)==0) then
 908          msrch(2)=msrch(2)+1
 909          srch(msrch(2),2)=srch(ii,1)
 910        end if
 911      end do
 912      !write(message,'(a,i0,a,i0,2a,i0)')&
 913      ! 'The second and third dimension of the FFT grid: ',ngfft(2),", ",ngfft(3),ch10,&
 914      ! 'were imposed to be multiple of the number of processors for the FFT: ', nproc_fft
 915      !if (ount /= dev_null) MSG_COMMENT(message)
 916    else
 917      msrch(2)=msrch(1)
 918      srch(:,2)=srch(:,1)
 919    end if
 920 
 921    ! The second and third search list have the same constraint
 922    msrch(3)=msrch(2)
 923    srch(:,3)=srch(:,2)
 924 
 925 !  The set of allowed ngfft values has been found
 926  end if ! first==1
 927 
 928 !Save input values of ngfft
 929  ngsav(1:3) = ngfft(1:3)
 930 
 931 !As an initial guess for ngfft, use the provided coarse mesh grid
 932  if (PRESENT(ngfftc)) then
 933    ngfft(1:3)=ngfftc(1:3)
 934    call wrtout(ount,' Using supplied coarse mesh as initial guess.','COLL')
 935  else
 936    ngfft(1:3)=2
 937  end if
 938 
 939 !Enlarge the initial guess until the set of ngfft entirely comprises the sphere
 940  do
 941 
 942    call bound(dsqmax,dsqmin,gbound,gmet,kpt,ngfft,plane)
 943 
 944    ! Exit the infinite do-loop if the sphere is inside the FFT box
 945    if (dsqmin>=(half*boxcutmin**2*ecut/pi**2)) exit
 946 
 947    ! Fix nearest boundary
 948    do ii=1,msrch(plane)-1
 949      if (srch(ii,plane)>=ngfft(plane)) then
 950 !      redefine ngfft(plane) to next higher choice
 951        ngfft(plane)=srch(ii+1,plane)
 952        exit ! Exit the loop over ii
 953      end if
 954 
 955      if (ii==msrch(plane)-1)then
 956        ! Here, we are in trouble
 957        write(message, '(a,i12,5a)' ) &
 958 &       'ngfft is bigger than allowed value =',ngfft(plane),'.',ch10,&
 959 &       'This indicates that desired ngfft is larger than getng',ch10,&
 960 &       'can handle. The code has to be changed and compiled.'
 961        MSG_BUG(message)
 962      end if
 963    end do
 964 
 965  end do ! End of the infinite do-loop : will either "exit", or abort
 966 
 967 !ecutmx=maximum ecut consistent with chosen ngfft
 968  ecutmx=0.5_dp*pi**2*dsqmin
 969 
 970 !Print results
 971  write(message, '(a,1p,e14.6,a,3i8,a,a,e14.6)' ) &
 972 & ' For input ecut=',ecut,' best grid ngfft=',ngfft(1:3),ch10,&
 973 & '       max ecut=',ecutmx
 974  call wrtout(ount,message,'COLL')
 975 
 976 ! The FFT grid is compatible with the symmetries if for each
 977 ! symmetry isym, each ii and each jj, the quantity
 978 ! (ngfft(jj)*symrel(jj,ii,isym))/ngfft(ii) is an integer.
 979 ! This relation is immediately verified for diagonal elements, since
 980 ! symrel is an integer. It is also verified if symrel(ii,jj,isym) is zero.
 981 
 982 !Compute the biggest (positive) common divisor of each off-diagonal element of the symmetry matrices
 983  divisor(:,:)=0; tobechecked=0
 984 
 985  do ii=1,3
 986    do jj=1,3
 987      if(jj==ii)cycle
 988      do isym=1,nsym
 989        if(symrel(jj,ii,isym)==0 .or. divisor(jj,ii)==1 )cycle
 990        tobechecked=1
 991        element=abs(symrel(jj,ii,isym))
 992        testdiv= ( divisor(jj,ii)==0 .or. divisor(jj,ii)==element .or. element==1)
 993        if(testdiv)then
 994          divisor(jj,ii)=element
 995        else
 996 !        Must evaluate common divisor between non-trivial numbers
 997          do
 998            if(divisor(jj,ii)<element)element=element-divisor(jj,ii)
 999            if(divisor(jj,ii)>element)divisor(jj,ii)=divisor(jj,ii)-element
1000            if(divisor(jj,ii)==element)exit
1001          end do
1002        end if
1003      end do
1004    end do
1005  end do
1006 
1007 !Check whether there is a problem
1008  testok=1
1009  if(tobechecked==1)then
1010    do ii=1,3
1011      do jj=1,3
1012        xx=divisor(jj,ii)*ngfft(jj)
1013        yy=xx/ngfft(ii)
1014        if(abs(yy-nint(yy))>tol8)testok=0
1015      end do
1016    end do
1017  end if
1018 
1019 !There is definitely a problem
1020  if(testok==0)then
1021 !  Use a brute force algorithm
1022 !  1) Because one knows that three identical numbers will satisfy
1023 !  the constraint, use the maximal ngfft value to define current triplet
1024 !  and associate total number of grid points
1025    ngcurrent(1:3)=maxval(ngfft(1:3))
1026 !  Takes into account the fact that ngfft(2) and ngfft(3) must
1027 !  be multiple of nproc_fft
1028    if(mod(ngcurrent(1),nproc_fft)/=0)ngcurrent(1:3)=ngcurrent(1:3)*max(1,nproc_fft)
1029    prodcurrent=ngcurrent(1)**3+1.0d-3
1030 !  2) Define maximal values for each component, limited
1031 !  by the maximal value of the list
1032    ngmax(1)=min(int(prodcurrent/(ngfft(2)*ngfft(3))),srch(msrch(1),1))
1033    ngmax(2)=min(int(prodcurrent/(ngfft(1)*ngfft(3))),srch(msrch(2),2))
1034    ngmax(3)=min(int(prodcurrent/(ngfft(1)*ngfft(2))),srch(msrch(3),3))
1035 !  3) Get minimal and maximal search indices
1036    do ii=1,3
1037      do isrch=1,msrch(ii)
1038        index=srch(isrch,ii)
1039        if(index==ngfft(ii))imin(ii)=isrch
1040 !      One cannot suppose that imax belongs to the allowed list,
1041 !      so must use <= instead of == , to determine largest index
1042        if(index<=ngmax(ii))imax(ii)=isrch
1043      end do
1044    end do
1045 !  4) Compute product of trial ngffts
1046 !  DEBUG
1047 !  write(ount,*)' getng : enter triple loop '
1048 !  write(ount,*)'imin',imin(1:3)
1049 !  write(ount,*)'imax',imax(1:3)
1050 !  write(ount,*)'ngcurrent',ngcurrent(1:3)
1051 !  ENDDEBUG
1052    do isrch1=imin(1),imax(1)
1053      ngtrial(1)=srch(isrch1,1)
1054      do isrch2=imin(2),imax(2)
1055        ngtrial(2)=srch(isrch2,2)
1056        do isrch3=imin(3),imax(3)
1057          ngtrial(3)=srch(isrch3,3)
1058          prodtrial=real(ngtrial(1))*real(ngtrial(2))*real(ngtrial(3))+1.0d-3
1059          if(prodtrial>prodcurrent-1.0d-4)exit
1060 !        The trial product is lower or equal to the current product,
1061 !        so now, checks whether the symmetry constraints are OK
1062          testok=1
1063          do ii=1,3
1064            do jj=1,3
1065              xx=divisor(jj,ii)*ngtrial(jj)
1066              yy=xx/ngtrial(ii)
1067              if(abs(yy-nint(yy))>tol8)testok=0
1068            end do
1069          end do
1070 !        DEBUG
1071 !        write(ount,'(a,3i6,a,i3,a,es16.6)' )' getng : current trial triplet',ngtrial(1:3),&
1072 !        &     ' testok=',testok,' prodtrial=',prodtrial
1073 !        ENDDEBUG
1074          if(testok==0)cycle
1075 !        The symmetry constraints are fulfilled, so update current values
1076          ngcurrent(1:3)=ngtrial(1:3)
1077          prodcurrent=prodtrial
1078        end do
1079      end do
1080    end do
1081 
1082    ngfft(1:3)=ngcurrent(1:3)
1083    call bound(dsqmax,dsqmin,gbound,gmet,kpt,ngfft,plane)
1084 !  ecutmx=maximum ecut consistent with chosen ngfft
1085    ecutmx=0.5_dp*pi**2*dsqmin
1086 !  Give results
1087    write(message, '(a,3i8,a,a,e14.6)' ) &
1088 &   ' However, must be changed due to symmetry =>',ngfft(1:3),ch10,&
1089 &   '       with max ecut=',ecutmx
1090    call wrtout(ount,message,'COLL')
1091 
1092    if (prodcurrent>huge(ii)) then
1093      write(message, '(5a)' )&
1094 &     'The best FFT grid will lead to indices larger',ch10,&
1095 &     'than the largest representable integer on this machine.',ch10,&
1096 &     'Action: try to deal with smaller problems. Also contact ABINIT group.'
1097      MSG_ERROR(message)
1098    end if
1099 
1100  end if ! testok==0
1101 
1102 !Possibly use the input values of ngfft
1103  if ( int( dble(ngsav(1)) / minbox ) >= ngfft(1) .and.&
1104 &     int( dble(ngsav(2)) / minbox ) >= ngfft(2) .and.&
1105 &     int( dble(ngsav(3)) / minbox ) >= ngfft(3) ) then
1106 
1107 !  Must check whether the values are in the allowed list
1108    testok=0
1109    do mu=1,3
1110      do ii=1,msrch(mu)
1111        if(srch(ii,mu)==ngsav(mu))then
1112          testok=testok+1
1113          exit
1114        end if
1115      end do
1116    end do
1117    if(testok==3)then
1118      write(message,'(a,3(a,i1,a,i3),a)') ' input values of',&
1119 &     (' ngfft(',mu,') =',ngsav(mu),mu=1,3),' are alright and will be used'
1120      call wrtout(ount,message,'COLL')
1121      do mu = 1,3
1122        ngfft(mu) = ngsav(mu)
1123      end do
1124    end if
1125 
1126  end if
1127 
1128 !mgfft needs to be set to the maximum of ngfft(1),ngfft(2),ngfft(3)
1129  mgfft = maxval(ngfft(1:3))
1130 
1131  if (paral_fft_==1) then
1132    ! For the time being, one need ngfft(2) and ngfft(3) to be multiple of nproc_fft
1133    if(modulo(ngfft(2),nproc_fft)/=0)then
1134      write(message,'(3a,i5,a,i5)')&
1135 &     'The second dimension of the FFT grid, ngfft(2), should be',&
1136 &     'a multiple of the number of processors for the FFT, nproc_fft.',&
1137 &     'However, ngfft(2)=',ngfft(2),' and nproc_fft=',nproc_fft
1138      MSG_BUG(message)
1139    end if
1140    if(modulo(ngfft(3),nproc_fft)/=0)then
1141      write(message,'(3a,i5,a,i5)')&
1142 &     'The third dimension of the FFT grid, ngfft(3), should be',&
1143 &     'a multiple of the number of processors for the FFT, nproc_fft.',&
1144 &     'However, ngfft(3)=',ngfft(3),' and nproc_fft=',nproc_fft
1145      MSG_BUG(message)
1146    end if
1147 
1148  else if (paral_fft_/=0) then
1149    write(message,'(a,i0)')'paral_fft_ should be 0 or 1, but its value is ',paral_fft_
1150    MSG_BUG(message)
1151  end if
1152 
1153 ! Compute effective number of FFT points (for this MPI node if parall FFT)
1154  nfft=product(ngfft(1:3))/max(1,nproc_fft)
1155 
1156 !Set up fft array dimensions ngfft(4,5,6) to avoid cache conflicts
1157  ngfft(4)=2*(ngfft(1)/2)+1
1158  ngfft(5)=2*(ngfft(2)/2)+1
1159  ngfft(6)=ngfft(3)
1160  if (any(fftalga == [FFT_FFTW3, FFT_DFTI])) then
1161    ! FFTW3 supports leading dimensions but at the price of a larger number of FFTs
1162    ! to be executed along z when the zero-padded version is used.
1163    ! One should check whether the augmentation is beneficial for FFTW3.
1164    ngfft(4)=2*(ngfft(1)/2)+1
1165    ngfft(5)=2*(ngfft(2)/2)+1
1166    !ngfft(4)=ngfft(1)
1167    !ngfft(5)=ngfft(2)
1168    ngfft(6)=ngfft(3)
1169  end if
1170 
1171  if (present(use_gpu_cuda)) then
1172    if (use_gpu_cuda==1) then
1173      ngfft(4)=ngfft(1)
1174      ngfft(5)=ngfft(2)
1175      ngfft(6)=ngfft(3)
1176    end if
1177  end if
1178 
1179  ngfft(14:18)=0 ! ngfft(14) to be filled outside of getng
1180 
1181  if (paral_fft_==0) then
1182    ngfft(9)=0     ! paral_fft_
1183    ngfft(10)=1    ! nproc_fft
1184    ngfft(11)=0    ! me_fft
1185    ngfft(12)=0    ! n2proc
1186    ngfft(13)=0    ! n3proc
1187  else
1188    ngfft(9)=1     ! paral_fft_
1189    ngfft(10)=nproc_fft
1190    ngfft(11)=me_fft
1191    ngfft(12)=ngfft(2)/nproc_fft
1192    ngfft(13)=ngfft(3)/nproc_fft
1193  end if
1194 
1195 
1196  call print_ngfft(ngfft,header="FFT mesh",unit=ount,mode_paral="COLL")
1197 
1198 end subroutine getng

m_fftcore/indfftrisc [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 indfftrisc

FUNCTION

 Take the data for sphere boundary and list of planewave in sphere (kg_k), manipulate them
 for convenient use in fourwf, and output them in indpw_k

INPUTS

 gbound(2*mgfft+4)=sphere boundary data
 kg_k(3,npw_k)=reduced planewave coordinates
 mgfft=maximum size of 1D FFTs
 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 npw_k=number of G vectors in basis at this k point

OUTPUT

 indpw_k(4,npw_k)=array which gives fft box index for given basis sphere
   in a representation that is directly usable by sg_fftrisc.f
 ngb=number of FFTs along z

PARENTS

      m_sgfft

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

4171 subroutine indfftrisc(gbound,indpw_k,kg_k,mgfft,ngb,ngfft,npw_k)
4172 
4173 
4174 !This section has been created automatically by the script Abilint (TD).
4175 !Do not modify the following lines by hand.
4176 #undef ABI_FUNC
4177 #define ABI_FUNC 'indfftrisc'
4178 !End of the abilint section
4179 
4180  implicit none
4181 
4182 !Arguments ------------------------------------
4183 !scalars
4184  integer,intent(in) :: mgfft,npw_k
4185  integer,intent(out) :: ngb
4186 !arrays
4187  integer,intent(in) :: gbound(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
4188  integer,intent(out) :: indpw_k(4,npw_k)
4189 
4190 !Local variables-------------------------------
4191 !scalars
4192  integer :: g1,g2,i1,i2,i3,igb,index,ipw,n1,n2,n3
4193 !arrays
4194  integer,allocatable :: index2d(:,:)
4195 
4196 ! *************************************************************************
4197 
4198  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
4199 
4200 !First, generate a 2d index for each column of data
4201  ABI_ALLOCATE(index2d,(n1,n2))
4202  index2d(:,:)=0
4203  index=1
4204  igb=3
4205  do g2=0,gbound(2) ! g2max
4206    do g1=0,gbound(igb+1)  ! g1max
4207      index2d(g1+1,g2+1)=index
4208      index=index+1
4209    end do
4210    if(gbound(igb)<=-1)then ! g1min
4211      do g1=gbound(igb)+n1,n1-1
4212        index2d(g1+1,g2+1)=index
4213        index=index+1
4214      end do
4215    end if
4216    igb=igb+2
4217  end do
4218 
4219  if(gbound(1)<=-1)then ! g2min
4220    do g2=gbound(1)+n2,n2-1
4221      do g1=0,gbound(igb+1)
4222        index2d(g1+1,g2+1)=index
4223        index=index+1
4224      end do
4225      if(gbound(igb)<=-1)then
4226        do g1=gbound(igb)+n1,n1-1
4227          index2d(g1+1,g2+1)=index
4228          index=index+1
4229        end do
4230      end if
4231      igb=igb+2
4232    end do
4233  end if
4234 
4235  ngb=index-1
4236 
4237 
4238 !The 2d index has been generated
4239 !Now, contract indpw_k(1,ipw) and indpw_k(2,ipw) into indpw_k(4,ipw)
4240 !indpw_k(1,ipw) and indpw_k(2,ipw) are used to hold inverse of index2d,
4241 !and for them, the second index does not fill 1:npw . It is only
4242 !the number of z-transform FFTs.
4243 
4244 !$OMP PARALLEL DO PRIVATE(i1,i2,i3)
4245  do ipw=1,npw_k
4246    i1=kg_k(1,ipw); if(i1<0)i1=i1+n1 ; i1=i1+1
4247    i2=kg_k(2,ipw); if(i2<0)i2=i2+n2 ; i2=i2+1
4248    i3=kg_k(3,ipw); if(i3<0)i3=i3+n3 ; i3=i3+1
4249    indpw_k(4,ipw)=index2d(i1,i2)
4250    indpw_k(3,ipw)=i3
4251  end do
4252 
4253  do i1=1,n1
4254    do i2=1,n2
4255      index=index2d(i1,i2)
4256      if(index/=0)then
4257        indpw_k(1,index)=i1
4258        indpw_k(2,index)=i2
4259      end if
4260    end do
4261  end do
4262 
4263  ABI_DEALLOCATE(index2d)
4264 
4265 end subroutine indfftrisc

m_fftcore/kgindex [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 kgindex

FUNCTION

 Compute the index of each plane wave on a FFT grid.

INPUTS

  kg_k(3,npw_k)=dimensionless coords of G vecs (integer)
  mpi_enreg=information about MPI parallelization
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  npw_k=number of planewaves

OUTPUT

  indpw_k(npw_k)=linear list number (in fft box) of given G vector for the current processor (local adress)
                =0 if kg_k(ipw) is not treated by this procesor
  mask(npw_k)=True if  kg_k(ipw) belongs to this procesor, false otherwise.

NOTES

   mpi_enreg is not necessary in this case (the info is also in ngfft), but much more easy to read...

PARENTS

      m_fft_prof,m_gsphere,m_screening,m_shirley,m_wfd,prcref,prcref_PMA

CHILDREN

SOURCE

4900 subroutine kgindex(indpw_k,kg_k,mask,mpi_enreg,ngfft,npw_k)
4901 
4902 
4903 !This section has been created automatically by the script Abilint (TD).
4904 !Do not modify the following lines by hand.
4905 #undef ABI_FUNC
4906 #define ABI_FUNC 'kgindex'
4907 !End of the abilint section
4908 
4909  implicit none
4910 
4911 !Arguments ------------------------------------
4912 !scalars
4913  integer,intent(in) :: npw_k
4914  type(MPI_type),intent(in) :: mpi_enreg
4915 !arrays
4916  integer,intent(in) :: kg_k(3,npw_k),ngfft(18)
4917  integer,intent(out) :: indpw_k(npw_k)
4918  logical,intent(out) :: mask(npw_k)
4919 !Local variables-------------------------------
4920 !scalars
4921  integer :: ig,ig1,ig2,ig3,me_fft,n1,n2,n3,nd2
4922  character(len=500) :: msg
4923  !arrays
4924  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
4925  !integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
4926 
4927 ! *************************************************************************
4928 
4929  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
4930 
4931 !Use the following indexing (N means ngfft of the adequate direction)
4932 !0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= kg
4933 !1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index
4934 
4935  me_fft=mpi_enreg%me_fft
4936  nd2=(n2-1)/mpi_enreg%nproc_fft+1
4937 
4938  if (n2== mpi_enreg%distribfft%n2_coarse) then
4939    fftn2_distrib => mpi_enreg%distribfft%tab_fftdp2_distrib
4940    ffti2_local => mpi_enreg%distribfft%tab_fftdp2_local
4941  else if (n2 == mpi_enreg%distribfft%n2_fine) then
4942    fftn2_distrib => mpi_enreg%distribfft%tab_fftdp2dg_distrib
4943    ffti2_local => mpi_enreg%distribfft%tab_fftdp2dg_local
4944  else
4945    MSG_BUG("Unable to find an allocated distrib for this fft grid")
4946  end if
4947 
4948  !call ptabs_fourwf(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
4949 
4950  do ig=1,npw_k
4951    ig1=modulo(kg_k(1,ig),n1)
4952    ig2=modulo(kg_k(2,ig),n2)
4953    ig3=modulo(kg_k(3,ig),n3)
4954    if(me_fft==fftn2_distrib(ig2+1)) then
4955      ig2=ffti2_local(ig2+1) - 1
4956      indpw_k(ig)=ig1+1+n1*(ig2+nd2*ig3)
4957      mask(ig)=.true.
4958    else
4959      indpw_k(ig)=0
4960      mask(ig)=.false.
4961    end if
4962    if ( ANY(kg_k(:,ig)>ngfft(1:3)/2) .or. ANY(kg_k(:,ig)<-(ngfft(1:3)-1)/2) ) then
4963      write(msg,'(a,i0,a)')" The G-vector with ig: ",ig," falls outside the FFT box."
4964      MSG_ERROR(msg)
4965    end if
4966  end do
4967 
4968 end subroutine kgindex

m_fftcore/kpgcount [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 kpgcount

FUNCTION

  Give the minimum and maximum number of G vectors in each direction:
  for each k-point, compute the number of G vectors in each direction,
  then store the min and max value over the set of k-points.

INPUTS

  ecut=planewave kinetic energy cutoff (hartrees)
  exchn2n3d=if 1, n2 and n3 are exchanged
  gmet(3,3)=reciprocal space metric (bohr^-2)
  istwfk(nkpt)=option parameter that describes the storage of wfs
  kpt(3,nkpt)=reduced coords of k point (in terms of recip latt vecs)
  nkpt=number of k points

OUTPUT

  ngmax(3)=maximum number of G vectors in each direction (x,y,z)
  ngmin(3)=minimum number of G vectors in each direction (x,y,z)

NOTES

  This routine has been extracted from kpgsph...

PARENTS

      finddistrproc

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

4714 subroutine kpgcount(ecut,exchn2n3d,gmet,istwfk,kpt,ngmax,ngmin,nkpt)
4715 
4716 
4717 !This section has been created automatically by the script Abilint (TD).
4718 !Do not modify the following lines by hand.
4719 #undef ABI_FUNC
4720 #define ABI_FUNC 'kpgcount'
4721 !End of the abilint section
4722 
4723  implicit none
4724 
4725 !Arguments ------------------------------------
4726 !scalars
4727  integer,intent(in) :: exchn2n3d,nkpt
4728  integer,intent(out) :: ngmax(3),ngmin(3)
4729  real(dp),intent(in) :: ecut
4730 !arrays
4731  integer,intent(in) :: istwfk(nkpt)
4732  real(dp),intent(in) :: gmet(3,3),kpt(3,nkpt)
4733 
4734 !Local variables-------------------------------
4735 !scalars
4736  integer :: ii,ikpt,istwf_k,kmax,ng1,ng2,ng3,nmin
4737  real(dp) :: gmet_trace,gscut,xx
4738 !arrays
4739  integer :: ngrid(3),nmax(3)
4740  real(dp) :: minor(3),numer(3)
4741 
4742 ! *************************************************************************
4743 
4744  DBG_ENTER("COLL")
4745 
4746  gscut=0.5_dp*ecut*piinv**2
4747  gmet_trace=gmet(1,1)+gmet(2,2)+gmet(3,3)
4748  minor(1)=gmet(2,2)*gmet(3,3)-gmet(2,3)**2
4749  minor(2)=gmet(1,1)*gmet(3,3)-gmet(1,3)**2
4750  minor(3)=gmet(2,2)*gmet(1,1)-gmet(1,2)**2
4751  numer(1)=gmet(1,2)**2*gmet(3,3)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3)+gmet(1,3)**2*gmet(2,2)
4752  numer(2)=gmet(2,3)**2*gmet(1,1)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3)+gmet(2,1)**2*gmet(3,3)
4753  numer(3)=gmet(3,2)**2*gmet(1,1)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3)+gmet(1,3)**2*gmet(2,2)
4754 
4755  ngmin(:)=1000000;ngmax(:)=0
4756  do ikpt=1,nkpt
4757    istwf_k=istwfk(ikpt)
4758 
4759    do ii=1,3
4760      xx=gmet(ii,ii)*minor(ii)-numer(ii)
4761      if(xx<tol10*gmet_trace**3.or.minor(ii)<tol10*gmet_trace**2)then
4762        MSG_BUG('The metric tensor seem incorrect')
4763      end if
4764      kmax=sqrt(gscut*minor(ii)/xx)
4765      nmax(ii)=floor(kmax-kpt(ii,ikpt)+tol10)
4766      nmin=ceiling(-kmax-kpt(ii,ikpt)-tol10)
4767      ngrid(ii)=nmax(ii)-nmin+1
4768    end do
4769 
4770    ng1=ngrid(1);if(istwf_k==2.or.istwf_k==3) ng1=nmax(1)+1
4771    if(exchn2n3d==0)then
4772      ng3=ngrid(3)
4773      ng2=ngrid(2);if(istwf_k>=2) ng2=nmax(2)+1
4774      if(istwf_k>=2.and.istwf_k<=5) ng2=ng2-1
4775    else
4776      ng2=ngrid(2)
4777      ng3=ngrid(3);if(istwf_k>=2) ng3=nmax(3)+1
4778      if(istwf_k==2.or.istwf_k==3.or.istwf_k==6.or.istwf_k==7) ng3=ng3-1
4779    end if
4780 
4781    if(ng1<ngmin(1)) ngmin(1)=ng1
4782    if(ng2<ngmin(2)) ngmin(2)=ng2
4783    if(ng3<ngmin(3)) ngmin(3)=ng3
4784    if(ng1>ngmax(1)) ngmax(1)=ng1
4785    if(ng2>ngmax(2)) ngmax(2)=ng2
4786    if(ng3>ngmax(3)) ngmax(3)=ng3
4787 
4788  end do
4789 
4790  DBG_EXIT("COLL")
4791 
4792 end subroutine kpgcount

m_fftcore/kpgsph [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 kpgsph

FUNCTION

 Use reciprocal space metric gmet(3,3) to set up the list
 of G vectors inside a sphere out to $ (1/2)*(2*\pi*(k+G))^2=ecut $.
 If mkmem=0 and mpw=0, then only count the number of planewaves

INPUTS

  ecut=planewave kinetic energy cutoff (hartrees)
  exchn2n3d=if 1, n2 and n3 are exchanged
  gmet(3,3)=reciprocal space metric (bohr^-2)
  ikg=shift to be given to the location of the output data in the array kg
  ikpt=number of the k-point
  istwf_k=option parameter that describes the storage of wfs
  kpt(3)=reduced coords of k point (in terms of recip latt vecs)
  mkmem =maximum number of k points which can fit in core memory
  mpi_enreg=information about MPI parallelization
  mpw=maximum number of planewaves as dimensioned in calling routine

OUTPUT

  kg(3,mpw*mkmem)=dimensionless coords of resulting G vecs (integer)
  npw=resulting number of planewaves inside ecut centered at kpt

SIDE EFFECTS

  mpi_enreg
    %me_g0=if 1, the plane wave G(0 0 0) is in the set of plane waves (and is the first)
    TODO: other SIDE EFFECTS on mpi_enreg should be described !!!

NOTES

  Must take into account the time-reversal symmetry when istwf_k is not 1.

PARENTS

      getmpw,initberry,initorbmag,kpgio,ks_ddiago,m_fft,m_fftcore,m_gsphere
      m_hamiltonian,m_wfd,mkpwind_k,wfconv

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

4313 subroutine kpgsph(ecut,exchn2n3d,gmet,ikg,ikpt,istwf_k,kg,kpt,mkmem,mpi_enreg,mpw,npw)
4314 
4315 
4316 !This section has been created automatically by the script Abilint (TD).
4317 !Do not modify the following lines by hand.
4318 #undef ABI_FUNC
4319 #define ABI_FUNC 'kpgsph'
4320 !End of the abilint section
4321 
4322  implicit none
4323 
4324 !Arguments ------------------------------------
4325 !scalars
4326  integer,intent(in) :: exchn2n3d,ikg,ikpt,istwf_k,mkmem,mpw
4327  integer,intent(out) :: npw
4328  real(dp),intent(in) :: ecut
4329  type(MPI_type),intent(inout) :: mpi_enreg
4330 !arrays
4331  integer,intent(inout) :: kg(3,mpw*mkmem)
4332  real(dp),intent(in) :: gmet(3,3),kpt(3)
4333 
4334 !Local variables-------------------------------
4335 !scalars
4336  integer :: i1,ig,ig1p,ig1pmax,ig2,ig2p,ig2pmax,ig2pmin,ig3,ig3p,ig3pmax
4337  integer :: ig3pmin,igtot,ii,ikpt_this_proc,in,ind,np_band,np_fft,npw_before,npw_remain,npw_split
4338  integer, save :: alloc_size=0
4339  real(dp) :: gap_pw,gmet11,gmet_trace,gmin,gs_fact,gs_part,gscut,v1,v2,v3,xx
4340  logical :: ipw_ok
4341  character(len=500) :: message
4342 !arrays
4343  integer :: ngrid(3),nmax(3),nmin(3),n2,ierr
4344  integer,allocatable :: array_ipw(:),ig1arr(:),ig2arr(:)
4345  integer,allocatable :: ig3arr(:),kg_ind(:),kg_small(:,:)
4346  integer, allocatable :: npw_gather(:),npw_disp(:) ,kg_ind_gather(:),kg_small_gather(:,:)
4347  real(dp) :: kmax(3),minor(3),numer(3),tsec(2)
4348  real(dp),allocatable :: kg1arr(:),kg2arr(:),kg3arr(:)
4349 
4350 ! *************************************************************************
4351 
4352  DBG_ENTER("COLL")
4353 
4354  call timab(23,1,tsec)
4355  if(istwf_k<1 .or. istwf_k>9)then
4356    write(message,'(3a,i0,a)' )&
4357 &   'The variable istwf_k must be between 1 and 9, while',ch10,&
4358 &   'the argument of the routine istwf_k =',istwf_k,'.'
4359    MSG_BUG(message)
4360  end if
4361 
4362  if(ikg+mpw>mkmem*mpw)then
4363    write(message,'(5a,i0,a,i0,a,i0,4a)')&
4364 &   'The variables ikg, mkmem, and mpw  must satisfy ikg<=(mkmem-1)*mpw,',ch10,&
4365 &   'while the arguments of the routine are',ch10,&
4366 &   'ikg =',ikg,', mkmem =',mkmem,', mpw =',mpw,ch10,&
4367 &   'Probable cause: Known error in invars1 for parallel spin-polarized case.',ch10,&
4368 &   'Temporary solution: Change the number of parallel processes.'
4369    MSG_BUG(message)
4370  end if
4371 
4372  np_band=0
4373  if (mpw>0) then
4374    np_band=1; if(mpi_enreg%paral_kgb==1) np_band=max(1,mpi_enreg%nproc_band)
4375    alloc_size=max(alloc_size,(mpw+1)*np_band)
4376    ABI_ALLOCATE(kg_small,(3, alloc_size))
4377    ABI_ALLOCATE(kg_ind,(alloc_size))
4378    kg_ind(:)=0
4379  end if
4380 
4381 !A larger array, that will be split on the correct processor
4382 !G**2 cutoff, gscut=Ecut/2 /Pi^2
4383 
4384  gscut=0.5_dp*ecut*piinv**2
4385 
4386 !In reduced coordinates, determine maximal value of k+G and G for each direction
4387 
4388  minor(1)=gmet(2,2)*gmet(3,3)-gmet(2,3)**2
4389  numer(1)=gmet(1,2)**2*gmet(3,3)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3) +gmet(1,3)**2*gmet(2,2)
4390  minor(2)=gmet(1,1)*gmet(3,3)-gmet(1,3)**2
4391  numer(2)=gmet(2,3)**2*gmet(1,1)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3) +gmet(2,1)**2*gmet(3,3)
4392  minor(3)=gmet(2,2)*gmet(1,1)-gmet(1,2)**2
4393  numer(3)=gmet(3,2)**2*gmet(1,1)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3) +gmet(1,3)**2*gmet(2,2)
4394 
4395 !Take the trace of the gmet tensor as dimensional reference
4396  gmet_trace=gmet(1,1)+gmet(2,2)+gmet(3,3)
4397 
4398  do ii=1,3
4399    xx=gmet(ii,ii)*minor(ii)-numer(ii)
4400    if(xx<tol10*gmet_trace**3 .or. minor(ii)<tol10*gmet_trace**2)then
4401      MSG_BUG('The metric tensor seem incorrect')
4402    end if
4403    kmax(ii)=sqrt(gscut*minor(ii)/xx)
4404    nmax(ii)=floor(kmax(ii)-kpt(ii)+tol10)
4405    nmin(ii)=ceiling(-kmax(ii)-kpt(ii)-tol10)
4406    ngrid(ii)=nmax(ii)-nmin(ii)+1
4407  end do
4408 !perform looping over fft box grid of size ngfft(1)*ngfft(2)*ngfft(3):
4409  ig=0;ind=0
4410  in=0
4411  gmet11=gmet(1,1)
4412 
4413 !Set up standard search sequence for grid points, in standard storage mode :
4414 !0 1 2 3 ... nmax nmin ... -1
4415 !If the mode is not standard, then some part of the FFT grid must be selected
4416 !
4417  ABI_ALLOCATE(ig1arr,(ngrid(1)))
4418  ABI_ALLOCATE(ig2arr,(ngrid(2)))
4419  ABI_ALLOCATE(ig3arr,(ngrid(3)))
4420  ABI_ALLOCATE(kg1arr,(ngrid(1)))
4421  ABI_ALLOCATE(kg2arr,(ngrid(2)))
4422  ABI_ALLOCATE(kg3arr,(ngrid(3)))
4423 
4424  do ig1p=1,ngrid(1)
4425    ig1arr(ig1p)=ig1p-1
4426    if (ig1p-1>nmax(1)) ig1arr(ig1p)=ig1p-ngrid(1)-1
4427    kg1arr(ig1p)=kpt(1)+dble(ig1arr(ig1p))
4428  end do
4429 
4430 !For the second direction, the number of points might depend on istwf_k
4431 !---------------------------------------------------------------------
4432  ig2pmax=ngrid(2)
4433  if(istwf_k>=2 .and. exchn2n3d==0) ig2pmax=nmax(2)+1
4434  ABI_ALLOCATE(array_ipw,(-ig2pmax:ig2pmax))
4435  array_ipw(:)=0
4436  do ig2p=1,ig2pmax
4437    ig2arr(ig2p)=ig2p-1
4438    if (ig2p-1>nmax(2)) ig2arr(ig2p)=ig2p-ngrid(2)-1
4439    kg2arr(ig2p)=kpt(2)+dble(ig2arr(ig2p))
4440  end do
4441 
4442 !For the third direction, the number of points might depend on istwf_k
4443 !---------------------------------------------------------------------
4444  ig3pmax=ngrid(3)
4445  if (istwf_k>=2 .and. exchn2n3d==1) ig3pmax=nmax(3)+1
4446 
4447  do ig3p=1,ig3pmax
4448    ig3arr(ig3p)=ig3p-1
4449    if(ig3p-1>nmax(3)) ig3arr(ig3p)=ig3p-ngrid(3)-1
4450    kg3arr(ig3p)=kpt(3)+dble(ig3arr(ig3p))
4451  end do
4452 
4453 !Performs loop on all grid points.
4454 !---------------------------------------------------------------------
4455  igtot = 0
4456  if(exchn2n3d==0)then
4457    mpi_enreg%me_g0=0
4458    do ig3p=1,ngrid(3)
4459      ig3=ig3arr(ig3p)
4460      v3=kg3arr(ig3p)
4461      ig2pmin=1
4462      if( istwf_k>=2 .and. istwf_k<=5 .and. ig3<0)then
4463        ig2pmin=2
4464      end if
4465 !    ig2pmax was initialized previously
4466      do ig2p=ig2pmin,ig2pmax
4467        ig2=ig2arr(ig2p)
4468 !      PAY ATTENTION : the proc 0 must have me_g0=1
4469        ipw_ok = .true.
4470        if(mpi_enreg%paral_kgb==1 ) then
4471           n2 =mpi_enreg%distribfft%n2_coarse
4472           ipw_ok = ipw_ok.and.(mpi_enreg%me_fft == mpi_enreg%distribfft%tab_fftwf2_distrib( modulo(ig2,n2) + 1))
4473        end if
4474        if (ig2==0 .and. ipw_ok) mpi_enreg%me_g0=1
4475        v2=kg2arr(ig2p)
4476        gs_part=gmet(2,2)*v2*v2+gmet(3,3)*v3*v3+2.0_dp*gmet(2,3)*v2*v3
4477        gs_fact=2.0_dp*(gmet(1,2)*v2+gmet(3,1)*v3)
4478        ig1pmax=ngrid(1)
4479        if( (istwf_k==2.or.istwf_k==3) .and. ig3p==1 .and. ig2p==1)ig1pmax=nmax(1)+1
4480        do ig1p=1,ig1pmax
4481          v1=kg1arr(ig1p)
4482          gmin=gs_part+v1*(gs_fact+v1*gmet11)
4483 !        If inside sphere:
4484          if (gmin<=gscut) then
4485            if (ipw_ok) then
4486              ig=ig+1  ! inside sphere
4487              igtot=igtot+1
4488              if (mpw>0.and.ig<=alloc_size) then
4489 !              Keep coords of pw:
4490                kg_small(1,ig)=ig1arr(ig1p)
4491                kg_small(2,ig)=ig2
4492                kg_small(3,ig)=ig3
4493                kg_ind(ig)=igtot
4494              end if
4495              array_ipw(ig2)=array_ipw(ig2)+1
4496            else
4497              igtot=igtot+1
4498            end if
4499          end if
4500        end do !ig1p
4501      end do !ig2p
4502    end do !ig3p
4503 
4504  else ! if (exchn2n3d/=0)
4505 
4506 !  ig2pmax was initialized previously
4507    mpi_enreg%me_g0=0
4508    do ig2p=1,ngrid(2)
4509      ig2=ig2arr(ig2p)
4510 !    PAY ATTENTION : the proc 0 must have me_g0=1
4511      ipw_ok = .true.
4512      if(mpi_enreg%paral_kgb==1 ) then
4513         n2 =mpi_enreg%distribfft%n2_coarse
4514         ipw_ok = ipw_ok.and.(mpi_enreg%me_fft == mpi_enreg%distribfft%tab_fftwf2_distrib( modulo(ig2,n2) + 1))
4515      end if
4516      if(ig2==0 .and. istwf_k>=2 .and. ipw_ok) mpi_enreg%me_g0=1
4517      v2     =kg2arr(ig2p)
4518      ig3pmin=1
4519      if( (istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7) .and. ig2<0)then
4520        ig3pmin=2
4521      end if
4522      do ig3p=ig3pmin,ig3pmax
4523        ig3=ig3arr(ig3p)
4524        v3=kg3arr(ig3p)
4525        gs_part=gmet(2,2)*v2*v2+gmet(3,3)*v3*v3+2.0_dp*gmet(2,3)*v2*v3
4526        gs_fact=2.0_dp*(gmet(1,2)*v2+gmet(3,1)*v3)
4527        ig1pmax=ngrid(1)
4528        if( (istwf_k==2.or.istwf_k==3) .and. ig3p==1 .and. ig2p==1)ig1pmax=nmax(1)+1
4529        do ig1p=1,ig1pmax
4530          v1=kg1arr(ig1p)
4531          gmin=gs_part+v1*(gs_fact+v1*gmet11)
4532 !        If inside sphere:
4533          if (gmin<=gscut) then
4534            if (ipw_ok) then
4535              ig=ig+1  ! inside sphere
4536              igtot=igtot+1
4537 !            Make sure not to overrun array, or simply do not store if mpw=0
4538              if (mpw>0.and.ig<=alloc_size) then
4539 !              Keep coords of pw:
4540                kg_small(1,ig)=ig1arr(ig1p)
4541                kg_small(2,ig)=ig2
4542                kg_small(3,ig)=ig3
4543                kg_ind(ig)=igtot
4544              end if
4545            else
4546              igtot=igtot+1
4547            end if
4548          end if
4549 
4550        end do ! ig1p
4551      end do ! ig3p
4552 !    end if ! if the ig2 plane is to be treated by this processor
4553    end do ! ig2p
4554 
4555  end if ! exchn2n3d==0 or ==1
4556 
4557 !Total number of G vectors at this k point is assigned: npw
4558 !when getcell = 1 it can be that ig exceeds mpw, the bound on kp_small
4559 !here is a workaround:
4560  if (mpw*np_band > 0 .and. ig > alloc_size) then
4561    npw = mpw*np_band
4562  else
4563    npw=ig
4564  end if
4565  alloc_size = max(alloc_size,npw)
4566 
4567  ABI_DEALLOCATE(ig1arr)
4568  ABI_DEALLOCATE(ig2arr)
4569  ABI_DEALLOCATE(ig3arr)
4570  ABI_DEALLOCATE(kg1arr)
4571  ABI_DEALLOCATE(kg2arr)
4572  ABI_DEALLOCATE(kg3arr)
4573 
4574 !BandFFT: plane-wave load balancing
4575  if (mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_fft>1.and. &
4576 &    mpi_enreg%pw_unbal_thresh>zero.and. &
4577 &    istwf_k==1) then
4578 !  Check for reequilibration
4579    np_fft=max(1,mpi_enreg%nproc_fft)
4580    ABI_ALLOCATE(npw_gather,(np_fft)) ! Count pw before balancing
4581    call xmpi_allgather(npw,npw_gather,mpi_enreg%comm_fft,ierr)
4582    gap_pw = 100._dp*(maxval(npw_gather(:))-minval(npw_gather))/(1.*sum(npw_gather(:))/np_fft)
4583    write(message,'(a,f5.2)' ) ' Relative gap for number of plane waves between process (%): ',gap_pw
4584    call wrtout(std_out,message,'COLL')
4585    if(gap_pw > mpi_enreg%pw_unbal_thresh) then ! Effective reequilibration
4586      write(message,'(a,f5.2,a,i4,a,f5.2,a)') &
4587 &        'Plane-wave unbalancing (',gap_pw,'%) for kpt ',ikpt,' is higher than threshold (',&
4588 &        mpi_enreg%pw_unbal_thresh,'%); a plane-wave balancing procedure is activated!'
4589      call wrtout(std_out,message,'COLL')
4590      !Get optimal number
4591      npw_split=sum(npw_gather(:))
4592      npw=npw_split/np_fft
4593      npw_remain=modulo(npw_split,np_fft)
4594      if(mpi_enreg%me_fft < npw_remain) npw=npw+1
4595      ig=npw
4596 #ifdef DEBUG_MODE
4597      write(message,*) 'New npw_fft = ', npw
4598      call wrtout(std_out,message,'COLL')
4599 #endif
4600      alloc_size = max(alloc_size,npw)
4601      if(mpw>0 ) then  !Step for requilibration between fft process
4602        ABI_ALLOCATE(npw_disp,(np_fft))
4603        npw_disp=0
4604        do i1=2,np_fft
4605          npw_disp(i1) = npw_disp(i1-1) + npw_gather(i1-1)
4606        end do
4607        ABI_ALLOCATE(kg_ind_gather,(npw_split))
4608        ABI_ALLOCATE(kg_small_gather,(3,npw_split))
4609        call xmpi_allgatherv(kg_ind, npw_gather(mpi_enreg%me_fft+1) , &
4610 &            kg_ind_gather,npw_gather,npw_disp,mpi_enreg%comm_fft,ierr)
4611        call xmpi_allgatherv(kg_small,3*npw_gather(mpi_enreg%me_fft+1), &
4612 &            kg_small_gather,3*npw_gather, 3*npw_disp,mpi_enreg%comm_fft,ierr)
4613        npw_before=mpi_enreg%me_fft*(npw_split/np_fft)+min(npw_remain,mpi_enreg%me_fft)
4614        kg_small(:,1:npw)=kg_small_gather(:,npw_before+1:npw_before+npw)
4615        kg_ind(  1:npw)=kg_ind_gather(npw_before+1:npw_before+npw)
4616 #ifdef DEBUG_MODE
4617        call wrtout(std_out,"keeping values done",'COLL')
4618 #endif
4619        ABI_DEALLOCATE(npw_disp)
4620        ABI_DEALLOCATE(kg_ind_gather)
4621        ABI_DEALLOCATE(kg_small_gather)
4622      end if
4623    end if!End of reequilibration step for paral KGB
4624    ABI_DEALLOCATE(npw_gather)
4625  end if
4626 
4627 !BandFFT: band load balancing
4628  if(mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_band>0) then
4629    np_band=max(1,mpi_enreg%nproc_band)
4630    npw_split=ig;npw=npw_split/np_band
4631    npw_remain=modulo(npw_split,np_band)
4632    if(mpi_enreg%me_band < npw_remain) npw=npw+1
4633    if(mpw > 0) then ! This is for the case when we only compute npw and put mpw=0
4634      npw_before=mpi_enreg%me_band*(npw_split/np_band)+min(npw_remain,mpi_enreg%me_band)
4635      kg_small(:,1:npw)=kg_small(:,npw_before+1:npw_before+npw)
4636      kg_ind  (  1:npw)=kg_ind  (  npw_before+1:npw_before+npw)
4637    end if
4638  end if
4639  if(mpw > 0) then
4640    do i1=1,npw
4641      kg(:,i1+ikg)=kg_small(:,i1)
4642      if (allocated(mpi_enreg%my_kgtab)) then
4643        ikpt_this_proc=mpi_enreg%my_kpttab(ikpt)
4644        mpi_enreg%my_kgtab(i1,ikpt_this_proc) = kg_ind(i1)
4645      end if
4646    end do
4647    ABI_DEALLOCATE(kg_small)
4648    ABI_DEALLOCATE(kg_ind)
4649  end if
4650 
4651  ABI_DEALLOCATE(array_ipw)
4652 
4653 !Take care of the me_g0 flag
4654  if(mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_band>0) then
4655    if(mpi_enreg%me_band==0.and.mpi_enreg%me_g0==1) then
4656 !    In this case, the processors had the 0 G vector before the new distribution, and still keeps it
4657      mpi_enreg%me_g0=1
4658    else
4659 !    All other cases
4660      mpi_enreg%me_g0=0
4661    end if
4662  end if
4663 
4664 !Check that npw is not zero
4665  if(mpi_enreg%paral_kgb==1.and.npw==0) then
4666    write(message,'(5a)' )&
4667 &   'Please decrease the number of npband*npfft MPI processes!',ch10,&
4668 &   'One of the MPI process has no plane-wave to handle.',ch10,&
4669 &   'Action: decrease npband and/or npfft.'
4670    MSG_ERROR(message)
4671  endif
4672 
4673  call timab(23,2,tsec)
4674 
4675  DBG_EXIT("COLL")
4676 
4677 end subroutine kpgsph

m_fftcore/mpifft_collect_datar [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_collect_datar

FUNCTION

 Collect a real-space MPI-FFT distributed array on each proc.

INPUTS

  ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
  cplex=1 if real array, 2 for complex
  nfft=Number of FFT points treated by this MPI proc
  nspden=Second dimension of rhor
  rhor(cplex*nfft,nspden)=Array in real space (MPI-FFT distributed)
  fftn3_distrib(n3)=rank of the processors which own fft planes in 3rd dimension.
  fftn3_local(n3)=local i3 indices
  comm_fft=MPI-FFT communicator
  [master]=MPI rank, Optional. If present, the global array is available only on master node.

OUTPUT

   rhor_glob(cplex*nfft_tot,nspden)=Global array

PARENTS

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

5210 subroutine mpifft_collect_datar(ngfft,cplex,nfft,nspden,rhor,comm_fft,fftn3_distrib,ffti3_local,rhor_glob,master)
5211 
5212 
5213 !This section has been created automatically by the script Abilint (TD).
5214 !Do not modify the following lines by hand.
5215 #undef ABI_FUNC
5216 #define ABI_FUNC 'mpifft_collect_datar'
5217 !End of the abilint section
5218 
5219  implicit none
5220 
5221 !Arguments ------------------------------------
5222 !scalars
5223  integer,intent(in) :: cplex,nfft,nspden,comm_fft
5224  integer,optional,intent(in) :: master
5225 !arrays
5226  integer,intent(in) :: ngfft(18),fftn3_distrib(ngfft(3)),ffti3_local(ngfft(3))
5227  real(dp),intent(in) :: rhor(cplex*nfft,nspden)
5228  real(dp),intent(out) :: rhor_glob(cplex*product(ngfft(1:3)),nspden)
5229 
5230 !Local variables-------------------------------
5231  integer :: ispden,i1,i2,i3,me_fft,i3_local,my_fftbase,glob_fftbase
5232  integer :: n1,n2,n3,ierr,nfft_tot
5233 
5234 ! *************************************************************************
5235 
5236  nfft_tot = product(ngfft(1:3)); me_fft = xmpi_comm_rank(comm_fft)
5237  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
5238 
5239  if (nfft_tot == nfft) then
5240    ! full rhor on each node, just do a copy
5241    rhor_glob = rhor
5242  else
5243    ! if MPI-FFT we have to gather the full rhor on each node.
5244    rhor_glob = zero
5245    do ispden=1,nspden
5246      do i3=1,n3
5247        if (me_fft == fftn3_distrib(i3)) then
5248          i3_local = ffti3_local(i3)
5249          do i2=1,n2
5250            my_fftbase =   cplex * ( (i2-1)*n1 + (i3_local-1)*n1*n2 )
5251            glob_fftbase = cplex * ( (i2-1)*n1 + (i3-1)*n1*n2 )
5252            do i1=1,cplex * n1
5253              rhor_glob(i1+glob_fftbase,ispden) = rhor(i1+my_fftbase,ispden)
5254            end do
5255          end do
5256        end if
5257      end do
5258    end do
5259    if (present(master)) then
5260      call xmpi_sum_master(rhor_glob,master,comm_fft,ierr)
5261    else
5262      call xmpi_sum(rhor_glob,comm_fft,ierr)
5263    end if
5264  end if
5265 
5266 end subroutine mpifft_collect_datar

m_fftcore/mpifft_dbox2fg [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_dbox2fg

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

3690 pure subroutine mpifft_dbox2fg(n1,n2,n3,n4,nd2proc,n6,ndat,fftn2_distrib,ffti2_local,me_fft,workf,nfft,fofg)
3691 
3692 
3693 !This section has been created automatically by the script Abilint (TD).
3694 !Do not modify the following lines by hand.
3695 #undef ABI_FUNC
3696 #define ABI_FUNC 'mpifft_dbox2fg'
3697 !End of the abilint section
3698 
3699  implicit none
3700 
3701 !Arguments ------------------------------------
3702 !scalars
3703  integer,intent(in) :: n1,n2,n3,n4,nd2proc,n6,ndat,me_fft,nfft
3704 !arrays
3705  integer,intent(in) :: fftn2_distrib(n2),ffti2_local(n2)
3706  real(dp),intent(in) :: workf(2,n4,n6,nd2proc*ndat)
3707  real(dp),intent(out) :: fofg(2,nfft*ndat)
3708 
3709 !Local variables-------------------------------
3710  integer :: idat,i1,i2,i3,i2_local,i2_ldat,fgbase
3711  real(dp) :: xnorm
3712 
3713 ! *************************************************************************
3714 
3715  xnorm=one/dble(n1*n2*n3)
3716 
3717  ! Transfer fft output to the original fft box
3718  do idat=1,ndat
3719    do i2=1,n2
3720      if( fftn2_distrib(i2) == me_fft) then
3721        i2_local = ffti2_local(i2)
3722        i2_ldat = i2_local + (idat-1) * nd2proc
3723        do i3=1,n3
3724          fgbase = n1*(i2_local - 1 + nd2proc*(i3-1)) + (idat - 1) * nfft
3725          do i1=1,n1
3726            fofg(1,i1+fgbase)=workf(1,i1,i3,i2_ldat)*xnorm
3727            fofg(2,i1+fgbase)=workf(2,i1,i3,i2_ldat)*xnorm
3728          end do
3729        end do
3730      end if
3731    end do
3732  end do
3733 
3734 end subroutine mpifft_dbox2fg

m_fftcore/mpifft_dbox2fg_dpc [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_dbox2fg_dpc

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

3753 pure subroutine mpifft_dbox2fg_dpc(n1,n2,n3,n4,nd2proc,n6,ndat,fftn2_distrib,ffti2_local,me_fft,workf,nfft,fofg)
3754 
3755 
3756 !This section has been created automatically by the script Abilint (TD).
3757 !Do not modify the following lines by hand.
3758 #undef ABI_FUNC
3759 #define ABI_FUNC 'mpifft_dbox2fg_dpc'
3760 !End of the abilint section
3761 
3762  implicit none
3763 
3764 !Arguments ------------------------------------
3765 !scalars
3766  integer,intent(in) :: n1,n2,n3,n4,nd2proc,n6,ndat,me_fft,nfft
3767 !arrays
3768  integer,intent(in) :: fftn2_distrib(n2),ffti2_local(n2)
3769  complex(dpc),intent(in) :: workf(n4,n6,nd2proc*ndat)
3770  real(dp),intent(out) :: fofg(2,nfft*ndat)
3771 
3772 !Local variables-------------------------------
3773  integer :: idat,i1,i2,i3,i2_local,i2_ldat,fgbase
3774  real(dp) :: xnorm
3775 
3776 ! *************************************************************************
3777 
3778  xnorm=one/dble(n1*n2*n3)
3779 
3780  ! Transfer fft output to the original fft box
3781  do idat=1,ndat
3782    do i2=1,n2
3783      if( fftn2_distrib(i2) == me_fft) then
3784        i2_local = ffti2_local(i2)
3785        i2_ldat = i2_local + (idat-1) * nd2proc
3786        do i3=1,n3
3787          fgbase = n1*(i2_local - 1 + nd2proc*(i3-1)) + (idat - 1) * nfft
3788          do i1=1,n1
3789            fofg(1,i1+fgbase)=REAL (workf(i1,i3,i2_ldat))*xnorm
3790            fofg(2,i1+fgbase)=AIMAG(workf(i1,i3,i2_ldat))*xnorm
3791          end do
3792        end do
3793      end if
3794    end do
3795  end do
3796 
3797 end subroutine mpifft_dbox2fg_dpc

m_fftcore/mpifft_dbox2fr [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_dbox2fr

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

3816 pure subroutine mpifft_dbox2fr(n1,n2,n3,n4,n5,nd3proc,ndat,fftn3_distrib,ffti3_local,me_fft,workr,cplex,nfft,fofr)
3817 
3818 
3819 !This section has been created automatically by the script Abilint (TD).
3820 !Do not modify the following lines by hand.
3821 #undef ABI_FUNC
3822 #define ABI_FUNC 'mpifft_dbox2fr'
3823 !End of the abilint section
3824 
3825  implicit none
3826 
3827 !Arguments ------------------------------------
3828 !scalars
3829  integer,intent(in) :: n1,n2,n3,n4,n5,nd3proc,ndat,me_fft,nfft,cplex
3830 !!arrays
3831  integer,intent(in) :: fftn3_distrib(n3),ffti3_local(n3)
3832  real(dp),intent(in) :: workr(2,n4,n5,nd3proc*ndat)
3833  real(dp),intent(out) :: fofr(cplex*nfft*ndat)
3834 
3835 !Local variables-------------------------------
3836  integer :: idat,i1,i2,i3,i3_local,i3_ldat,frbase
3837 
3838 ! *************************************************************************
3839 
3840  select case (cplex)
3841  case (1)
3842 
3843    do idat=1,ndat
3844      do i3=1,n3
3845        if( fftn3_distrib(i3) == me_fft) then
3846          i3_local = ffti3_local(i3)
3847          i3_ldat = i3_local + (idat - 1) * nd3proc
3848          do i2=1,n2
3849            frbase=n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * nfft
3850            do i1=1,n1
3851              fofr(i1+frbase)=workr(1,i1,i2,i3_ldat)
3852            end do
3853          end do
3854        end if
3855      end do
3856    end do
3857 
3858  case (2)
3859 
3860    do idat=1,ndat
3861      do i3=1,n3
3862        if (fftn3_distrib(i3) == me_fft) then
3863          i3_local = ffti3_local(i3)
3864          i3_ldat = i3_local + (idat - 1) * nd3proc
3865          do i2=1,n2
3866            frbase=2*n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * cplex * nfft
3867            !if (frbase > cplex*nfft*ndat - 2*n1) then
3868            !   write(std_out,*)i2,i3_local,frbase,cplex*nfft*ndat
3869            !   MSG_ERROR("frbase")
3870            !end if
3871            do i1=1,n1
3872              fofr(2*i1-1+frbase)=workr(1,i1,i2,i3_ldat)
3873              fofr(2*i1  +frbase)=workr(2,i1,i2,i3_ldat)
3874            end do
3875          end do
3876        end if
3877      end do
3878    end do
3879 
3880  case default
3881    !MSG_BUG("Wrong cplex")
3882    fofr = huge(one)
3883  end select
3884 
3885 end subroutine mpifft_dbox2fr

m_fftcore/mpifft_dbox2fr_dpc [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_dbox2fr_dpc

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

3904 pure subroutine mpifft_dbox2fr_dpc(n1,n2,n3,n4,n5,nd3proc,ndat,fftn3_distrib,ffti3_local,me_fft,workr,cplex,nfft,fofr)
3905 
3906 
3907 !This section has been created automatically by the script Abilint (TD).
3908 !Do not modify the following lines by hand.
3909 #undef ABI_FUNC
3910 #define ABI_FUNC 'mpifft_dbox2fr_dpc'
3911 !End of the abilint section
3912 
3913  implicit none
3914 
3915 !Arguments ------------------------------------
3916 !scalars
3917  integer,intent(in) :: n1,n2,n3,n4,n5,nd3proc,ndat,me_fft,nfft,cplex
3918 !!arrays
3919  integer,intent(in) :: fftn3_distrib(n3),ffti3_local(n3)
3920  complex(dpc),intent(in) :: workr(n4,n5,nd3proc*ndat)
3921  real(dp),intent(out) :: fofr(cplex*nfft*ndat)
3922 
3923 !Local variables-------------------------------
3924  integer :: idat,i1,i2,i3,i3_local,i3_ldat,frbase
3925 
3926 ! *************************************************************************
3927 
3928  select case (cplex)
3929  case (1)
3930 
3931    do idat=1,ndat
3932      do i3=1,n3
3933        if( fftn3_distrib(i3) == me_fft) then
3934          i3_local = ffti3_local(i3)
3935          i3_ldat = i3_local + (idat - 1) * nd3proc
3936          do i2=1,n2
3937            frbase=n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * nfft
3938            do i1=1,n1
3939              fofr(i1+frbase)=REAL(workr(i1,i2,i3_ldat))
3940            end do
3941          end do
3942        end if
3943      end do
3944    end do
3945 
3946  case (2)
3947 
3948    do idat=1,ndat
3949      do i3=1,n3
3950        if (fftn3_distrib(i3) == me_fft) then
3951          i3_local = ffti3_local(i3)
3952          i3_ldat = i3_local + (idat - 1) * nd3proc
3953          do i2=1,n2
3954            frbase=2*n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * cplex * nfft
3955            do i1=1,n1
3956              fofr(2*i1-1+frbase)=REAL (workr(i1,i2,i3_ldat))
3957              fofr(2*i1  +frbase)=AIMAG(workr(i1,i2,i3_ldat))
3958            end do
3959          end do
3960        end if
3961      end do
3962    end do
3963 
3964  case default
3965    !MSG_BUG("Wrong cplex")
3966    fofr = huge(one)
3967  end select
3968 
3969 end subroutine mpifft_dbox2fr_dpc

m_fftcore/mpifft_fg2dbox [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_fg2dbox

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

3573 pure subroutine mpifft_fg2dbox(nfft,ndat,fofg,n1,n2,n3,n4,nd2proc,n6,fftn2_distrib,ffti2_local,me_fft,workf)
3574 
3575 
3576 !This section has been created automatically by the script Abilint (TD).
3577 !Do not modify the following lines by hand.
3578 #undef ABI_FUNC
3579 #define ABI_FUNC 'mpifft_fg2dbox'
3580 !End of the abilint section
3581 
3582  implicit none
3583 
3584 !Arguments ------------------------------------
3585 !scalars
3586  integer,intent(in) :: nfft,ndat,n1,n2,n3,n4,nd2proc,n6,me_fft
3587 !arrays
3588  integer,intent(in) :: fftn2_distrib(n2),ffti2_local(n2)
3589  real(dp),intent(in) :: fofg(2,nfft*ndat)
3590  real(dp),intent(inout) :: workf(2,n4,n6,nd2proc*ndat)
3591 
3592 !Local variables-------------------------------
3593  integer :: idat,i1,i2,i3,i2_local,i2_ldat,fgbase
3594 
3595 ! *************************************************************************
3596 
3597  do idat=1,ndat
3598    do i3=1,n3
3599      do i2=1,n2
3600        if (fftn2_distrib(i2) == me_fft) then
3601          i2_local = ffti2_local(i2)
3602          i2_ldat = i2_local + (idat-1) * nd2proc
3603          fgbase= n1*(i2_local-1 + nd2proc*(i3-1)) + (idat-1) * nfft
3604          do i1=1,n1
3605            workf(1,i1,i3,i2_ldat)=fofg(1,i1+fgbase)
3606            workf(2,i1,i3,i2_ldat)=fofg(2,i1+fgbase)
3607          end do
3608        end if
3609      end do
3610    end do
3611  end do
3612 
3613 end subroutine mpifft_fg2dbox

m_fftcore/mpifft_fg2dbox_dpc [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_fg2dbox_dpc

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

3632 pure subroutine mpifft_fg2dbox_dpc(nfft,ndat,fofg,n1,n2,n3,n4,nd2proc,n6,fftn2_distrib,ffti2_local,me_fft,workf)
3633 
3634 
3635 !This section has been created automatically by the script Abilint (TD).
3636 !Do not modify the following lines by hand.
3637 #undef ABI_FUNC
3638 #define ABI_FUNC 'mpifft_fg2dbox_dpc'
3639 !End of the abilint section
3640 
3641  implicit none
3642 
3643 !Arguments ------------------------------------
3644 !scalars
3645  integer,intent(in) :: nfft,ndat,n1,n2,n3,n4,nd2proc,n6,me_fft
3646 !arrays
3647  integer,intent(in) :: fftn2_distrib(n2),ffti2_local(n2)
3648  real(dp),intent(in) :: fofg(2,nfft*ndat)
3649  complex(dpc),intent(inout) :: workf(n4,n6,nd2proc*ndat)
3650 
3651 !Local variables-------------------------------
3652  integer :: idat,i1,i2,i3,i2_local,i2_ldat,fgbase
3653 
3654 ! *************************************************************************
3655 
3656  do idat=1,ndat
3657    do i3=1,n3
3658      do i2=1,n2
3659        if (fftn2_distrib(i2) == me_fft) then
3660          i2_local = ffti2_local(i2)
3661          i2_ldat = i2_local + (idat-1) * nd2proc
3662          fgbase= n1*(i2_local-1 + nd2proc*(i3-1)) + (idat-1) * nfft
3663          do i1=1,n1
3664            workf(i1,i3,i2_ldat)=CMPLX(fofg(1,i1+fgbase), fofg(2,i1+fgbase), kind=dpc)
3665          end do
3666        end if
3667      end do
3668    end do
3669  end do
3670 
3671 end subroutine mpifft_fg2dbox_dpc

m_fftcore/mpifft_fr2dbox [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_fr2dbox

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

3988 pure subroutine mpifft_fr2dbox(cplex,nfft,ndat,fofr,n1,n2,n3,n4,n5,nd3proc,fftn3_distrib,ffti3_local,me_fft,workr)
3989 
3990 
3991 !This section has been created automatically by the script Abilint (TD).
3992 !Do not modify the following lines by hand.
3993 #undef ABI_FUNC
3994 #define ABI_FUNC 'mpifft_fr2dbox'
3995 !End of the abilint section
3996 
3997  implicit none
3998 
3999 !Arguments ------------------------------------
4000 !scalars
4001  integer,intent(in) :: cplex,nfft,ndat,n1,n2,n3,n4,n5,nd3proc,me_fft
4002 !!arrays
4003  integer,intent(in) :: fftn3_distrib(n3),ffti3_local(n3)
4004  real(dp),intent(in) :: fofr(cplex*nfft*ndat)
4005  real(dp),intent(inout) :: workr(2,n4,n5,nd3proc*ndat)
4006 
4007 !Local variables-------------------------------
4008  integer :: idat,i1,i2,i3,i3_local,i3_ldat,frbase
4009 
4010 ! *************************************************************************
4011 
4012  select case (cplex)
4013  case (1)
4014 
4015    do idat=1,ndat
4016      do i3=1,n3
4017        if( me_fft == fftn3_distrib(i3) ) then
4018          i3_local = ffti3_local(i3)
4019          i3_ldat = i3_local + (idat-1) * nd3proc
4020          do i2=1,n2
4021            frbase=n1*(i2-1+n2*(i3_local-1)) + (idat-1) * nfft
4022            do i1=1,n1
4023              workr(1,i1,i2,i3_ldat)=fofr(i1+frbase)
4024              workr(2,i1,i2,i3_ldat)=zero
4025            end do
4026          end do
4027        end if
4028      end do
4029    end do
4030 
4031  case (2)
4032 
4033    do idat=1,ndat
4034      do i3=1,n3
4035        if( me_fft == fftn3_distrib(i3) ) then
4036          i3_local = ffti3_local(i3)
4037          i3_ldat = i3_local + (idat-1) * nd3proc
4038          do i2=1,n2
4039            frbase=2*n1*(i2-1+n2*(i3_local-1)) + (idat-1) * cplex * nfft
4040            do i1=1,n1
4041              workr(1,i1,i2,i3_ldat)=fofr(2*i1-1+frbase)
4042              workr(2,i1,i2,i3_ldat)=fofr(2*i1  +frbase)
4043            end do
4044          end do
4045        end if
4046      end do
4047    end do
4048 
4049  case default
4050    !MSG_BUG("Wrong cplex")
4051    workr = huge(one)
4052  end select
4053 
4054 end subroutine mpifft_fr2dbox

m_fftcore/mpifft_fr2dbox_dpc [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  mpifft_fr2dbox_dpc

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

4073 pure subroutine mpifft_fr2dbox_dpc(cplex,nfft,ndat,fofr,n1,n2,n3,n4,n5,nd3proc,fftn3_distrib,ffti3_local,me_fft,workr)
4074 
4075 
4076 !This section has been created automatically by the script Abilint (TD).
4077 !Do not modify the following lines by hand.
4078 #undef ABI_FUNC
4079 #define ABI_FUNC 'mpifft_fr2dbox_dpc'
4080 !End of the abilint section
4081 
4082  implicit none
4083 
4084 !Arguments ------------------------------------
4085 !scalars
4086  integer,intent(in) :: cplex,nfft,ndat,n1,n2,n3,n4,n5,nd3proc,me_fft
4087 !!arrays
4088  integer,intent(in) :: fftn3_distrib(n3),ffti3_local(n3)
4089  real(dp),intent(in) :: fofr(cplex*nfft*ndat)
4090  complex(dpc),intent(inout) :: workr(n4,n5,nd3proc*ndat)
4091 
4092 !Local variables-------------------------------
4093  integer :: idat,i1,i2,i3,i3_local,i3_ldat,frbase
4094 
4095 ! *************************************************************************
4096 
4097  select case (cplex)
4098  case (1)
4099 
4100    do idat=1,ndat
4101      do i3=1,n3
4102        if( me_fft == fftn3_distrib(i3) ) then
4103          i3_local = ffti3_local(i3)
4104          i3_ldat = i3_local + (idat-1) * nd3proc
4105          do i2=1,n2
4106            frbase=n1*(i2-1+n2*(i3_local-1)) + (idat-1) * nfft
4107            do i1=1,n1
4108              workr(i1,i2,i3_ldat)=CMPLX(fofr(i1+frbase), zero, kind=dpc)
4109            end do
4110          end do
4111        end if
4112      end do
4113    end do
4114 
4115  case (2)
4116 
4117    do idat=1,ndat
4118      do i3=1,n3
4119        if( me_fft == fftn3_distrib(i3) ) then
4120          i3_local = ffti3_local(i3)
4121          i3_ldat = i3_local + (idat-1) * nd3proc
4122          do i2=1,n2
4123            frbase=2*n1*(i2-1+n2*(i3_local-1)) + (idat-1) * cplex * nfft
4124            do i1=1,n1
4125              workr(i1,i2,i3_ldat)=CMPLX(fofr(2*i1-1+frbase), fofr(2*i1  +frbase), kind=dpc)
4126            end do
4127          end do
4128        end if
4129      end do
4130    end do
4131 
4132  case default
4133    !MSG_BUG("Wrong cplex")
4134    workr = huge(one)
4135  end select
4136 
4137 end subroutine mpifft_fr2dbox_dpc

m_fftcore/mpiswitch [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 mpiswitch

FUNCTION

INPUTS

OUTPUT

PARENTS

      m_fftw3,m_sg2002

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

3348 pure subroutine mpiswitch(j3,n1dfft,Jp2st,J2st,lot,n1,nd2proc,nd3proc,nproc,ioption,zmpi1,zw)
3349 
3350 
3351 !This section has been created automatically by the script Abilint (TD).
3352 !Do not modify the following lines by hand.
3353 #undef ABI_FUNC
3354 #define ABI_FUNC 'mpiswitch'
3355 !End of the abilint section
3356 
3357  implicit none
3358 
3359 !Arguments ------------------------------------
3360  integer,intent(in) :: j3,n1dfft,lot,n1,nd2proc,nd3proc,nproc,ioption
3361  integer,intent(inout) :: Jp2st,J2st
3362  real(dp),intent(in) :: zmpi1(2,n1,nd2proc,nd3proc,nproc)
3363  real(dp),intent(inout) :: zw(2,lot,n1)
3364 
3365 !Local variables-------------------------------
3366  integer :: Jp2,J2,I1,ind,jj2,mfft,jjp2
3367 
3368 ! *************************************************************************
3369  mfft=0
3370 
3371  if (ioption /= 1) then
3372    do Jp2=Jp2st,nproc
3373      do J2=J2st,nd2proc
3374        mfft=mfft+1
3375        if (mfft.gt.n1dfft) then
3376          Jp2st=Jp2
3377          J2st=J2
3378          return
3379        end if
3380        do I1=1,n1
3381          zw(1,mfft,I1)=zmpi1(1,I1,J2,j3,Jp2)
3382          zw(2,mfft,I1)=zmpi1(2,I1,J2,j3,Jp2)
3383        end do
3384      end do
3385      J2st=1
3386    end do
3387 
3388  else
3389    do Jp2=Jp2st,nproc
3390      do J2=J2st,nd2proc
3391        mfft=mfft+1
3392        if (mfft.gt.n1dfft) then
3393          Jp2st=Jp2
3394          J2st=J2
3395          return
3396        end if
3397        ind=(Jp2-1) * nd2proc + J2
3398        jj2=(ind-1)/nproc +1
3399 
3400        !jjp2=modulo(ind,nproc) +1
3401        jjp2=modulo(ind-1,nproc)+1
3402 
3403        !in other words: mfft=(jj2-1)*nproc+jjp2 (modulo case)
3404        !istead of mfft=(Jjp2-1) * nd2proc + Jj2 (slice case)
3405        !with 1<=jjp2<=nproc, jj2=1,nd2proc
3406        do I1=1,n1
3407          ! zw(1,mfft,I1)=zmpi1(1,I1,J2,j3,Jp2)
3408          ! zw(2,mfft,I1)=zmpi1(2,I1,J2,j3,Jp2)
3409          zw(1,mfft,I1)=zmpi1(1,I1,jj2,j3,jjp2)
3410          zw(2,mfft,I1)=zmpi1(2,I1,jj2,j3,jjp2)
3411        end do
3412      end do
3413      J2st=1
3414    end do
3415  end if
3416 
3417 end subroutine mpiswitch

m_fftcore/mpiswitch_cent [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 mpiswitch_cent

FUNCTION

   Perform the local rotation

     input: I1,J2,j3,Jp2,(jp3)
     output: J2,Jp2,I1,j3,(jp3)

   and fill the central region of the frequency spectrum with zeros

INPUTS

OUTPUT

PARENTS

      m_fftw3,m_sg2002

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

3447 pure subroutine mpiswitch_cent(j3,n1dfft,Jp2stb,J2stb,lot,max1,md1,m1,n1,md2proc,&
3448 &  nd3proc,nproc,ioption,zmpi1,zw,max2,m2,n2)
3449 
3450 
3451 !This section has been created automatically by the script Abilint (TD).
3452 !Do not modify the following lines by hand.
3453 #undef ABI_FUNC
3454 #define ABI_FUNC 'mpiswitch_cent'
3455 !End of the abilint section
3456 
3457  implicit none
3458 
3459 !Arguments ------------------------------------
3460  integer,intent(in) :: j3,n1dfft,lot,max1,md1,m1,n1,md2proc,nd3proc,nproc,ioption
3461  integer,intent(in) :: m2,max2,n2
3462  integer,intent(inout) :: Jp2stb,J2stb
3463  real(dp),intent(in) :: zmpi1(2,md1,md2proc,nd3proc,nproc)
3464  real(dp),intent(inout) :: zw(2,lot,n1)
3465 
3466 !Local variables-------------------------------
3467  integer :: mfft,jp2,j2,jjp2,jj2,i1,ind
3468 
3469 ! *************************************************************************
3470 
3471  ABI_UNUSED((/m2,max2,n2/))
3472 
3473  mfft=0
3474 
3475  if (ioption /= 1) then
3476    do Jp2=Jp2stb,nproc
3477      do J2=J2stb,md2proc
3478 
3479        mfft=mfft+1
3480        if (mfft.gt.n1dfft) then
3481          Jp2stb=Jp2
3482          J2stb=J2
3483          !MSG_WARNING("Returning from mpiswithc_cent")
3484          return
3485        end if
3486 
3487        ! Here, zero and positive frequencies
3488        ! In zmpi1, they are stored from 1 to max1+1
3489        do I1=1,max1+1
3490          zw(1,mfft,I1)=zmpi1(1,I1,J2,j3,Jp2)
3491          zw(2,mfft,I1)=zmpi1(2,I1,J2,j3,Jp2)
3492        end do
3493 
3494        ! Fill the center region with zeros
3495        do I1=max1+2,n1-m1+max1+1
3496          zw(1,mfft,I1)=zero
3497          zw(2,mfft,I1)=zero
3498        end do
3499 
3500        ! Here, negative frequencies
3501        ! In zmpi1, they are stored from 1 to m1half
3502        do I1=max1+2,m1
3503          zw(1,mfft,I1+n1-m1)=zmpi1(1,I1,J2,j3,Jp2)
3504          zw(2,mfft,I1+n1-m1)=zmpi1(2,I1,J2,j3,Jp2)
3505        end do
3506      end do
3507      J2stb=1
3508    end do
3509 
3510  else
3511    do Jp2=Jp2stb,nproc
3512      do J2=J2stb,md2proc
3513 
3514        mfft=mfft+1
3515        if (mfft.gt.n1dfft) then
3516          Jp2stb=Jp2
3517          J2stb=J2
3518          !MSG_WARNING("Returning from mpiswithc_cent")
3519          return
3520        end if
3521 
3522        ind=(Jp2-1) * md2proc + J2
3523        jj2=(ind-1)/nproc +1
3524 
3525        !jjp2=modulo(ind,nproc) +1
3526        jjp2=modulo(ind-1,nproc)+1
3527 
3528        ! I gather consecutive I2 indexes in mfft in the modulo case
3529        ! Here, zero and positive frequencies
3530        ! In zmpi1, they are stored from 1 to max1+1
3531        do I1=1,max1+1
3532          zw(1,mfft,I1)=zmpi1(1,I1,Jj2,j3,Jjp2)
3533          zw(2,mfft,I1)=zmpi1(2,I1,Jj2,j3,Jjp2)
3534        end do
3535 
3536        ! Fill the center region with zeros
3537        do I1=max1+2,n1-m1+max1+1
3538          zw(1,mfft,I1)=zero
3539          zw(2,mfft,I1)=zero
3540        end do
3541 
3542        ! Here, negative frequencies
3543        ! In zmpi1, they are stored from 1 to m1half
3544        do I1=max1+2,m1
3545          zw(1,mfft,I1+n1-m1)=zmpi1(1,I1,Jj2,j3,Jjp2)
3546          zw(2,mfft,I1+n1-m1)=zmpi1(2,I1,Jj2,j3,Jjp2)
3547        end do
3548 
3549      end do
3550      J2stb=1
3551    end do
3552  end if
3553 
3554 end subroutine mpiswitch_cent

m_fftcore/multpot [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  multpot

FUNCTION

INPUTS

  icplexwf=1 if u(r) is real, 2 otherwise.
  icplex=1 if v(r) is real, 2 otherwise.
  includelast
  nd1,nd2=Leading dimensions of pot(icplex*nd1,nd2)
  n2
  lot
  n1dfft

OUTPUT

PARENTS

      m_fftw3,m_sg2002

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

5084 subroutine multpot(icplexwf,icplex,includelast,nd1,nd2,n2,lot,n1dfft,pot,zw)
5085 
5086 
5087 !This section has been created automatically by the script Abilint (TD).
5088 !Do not modify the following lines by hand.
5089 #undef ABI_FUNC
5090 #define ABI_FUNC 'multpot'
5091 !End of the abilint section
5092 
5093  implicit none
5094 
5095  !Arguments ------------------------------------
5096  integer,intent(in) :: icplexwf,icplex,includelast,nd1,nd2,n2,lot,n1dfft
5097  real(dp),intent(in) :: pot(icplex*nd1,nd2)
5098  real(dp),intent(inout) :: zw(2,lot,n2)
5099 
5100 !Local variables-------------------------------
5101  integer :: i2,j
5102  real(dp) :: rew,imw
5103 
5104 ! *************************************************************************
5105 
5106  if (icplexwf==1) then
5107    ! Real u(r)
5108 
5109    if (icplex==2) then
5110      MSG_BUG('multpot: icplexwf=1 and icplex=2')
5111    else
5112      ! TO BE SPEEDED UP : should use the same trick as Stefan
5113      if(includelast==1)then
5114        do i2=1,n2
5115          do j=1,n1dfft
5116            zw(1,j,i2)=zw(1,j,i2)*pot(2*j-1,i2)
5117            zw(2,j,i2)=zw(2,j,i2)*pot(2*j  ,i2)
5118          end do
5119        end do
5120      else
5121        do i2=1,n2
5122          do j=1,n1dfft-1
5123            zw(1,j,i2)=zw(1,j,i2)*pot(2*j-1,i2)
5124            zw(2,j,i2)=zw(2,j,i2)*pot(2*j  ,i2)
5125          end do
5126          zw(1,n1dfft,i2)=zw(1,n1dfft,i2)*pot(2*n1dfft-1,i2)
5127        end do
5128      end if
5129    end if
5130 
5131  else if (icplexwf==2) then
5132    ! Complex u(r)
5133 
5134    if (icplex==1) then
5135 
5136      do i2=1,n2-1,2
5137        do j=1,n1dfft
5138          zw(1,j,i2)=zw(1,j,i2)*pot(j,i2)
5139          zw(2,j,i2)=zw(2,j,i2)*pot(j,i2)
5140          zw(1,j,i2+1)=zw(1,j,i2+1)*pot(j,i2+1)
5141          zw(2,j,i2+1)=zw(2,j,i2+1)*pot(j,i2+1)
5142        end do
5143      end do
5144 
5145      if (2*(n2/2)/=n2) then
5146        do j=1,n1dfft
5147          zw(1,j,n2)=zw(1,j,n2)*pot(j,n2)
5148          zw(2,j,n2)=zw(2,j,n2)*pot(j,n2)
5149        end do
5150      end if
5151 
5152    else
5153 
5154      do i2=1,n2-1,2
5155        do j=1,n1dfft
5156          rew = zw(1,j,i2); imw = zw(2,j,i2)
5157          zw(1,j,i2) = rew*pot(2*j-1,i2) - imw*pot(2*j,i2)
5158          zw(2,j,i2) = imw*pot(2*j-1,i2) + rew*pot(2*j,i2)
5159 
5160          rew = zw(1,j,i2+1); imw = zw(2,j,i2+1)
5161          zw(1,j,i2+1) = rew*pot(2*j-1,i2+1) - imw*pot(2*j,i2+1)
5162          zw(2,j,i2+1) = imw*pot(2*j-1,i2+1) + rew*pot(2*j,i2+1)
5163        end do
5164      end do
5165 
5166      if (2*(n2/2)/=n2) then
5167        do j=1,n1dfft
5168          rew = zw(1,j,n2); imw = zw(2,j,n2)
5169          zw(1,j,n2) = rew*pot(2*j-1,n2) - imw*pot(2*j,n2)
5170          zw(2,j,n2) = imw*pot(2*j-1,n2) + rew*pot(2*j,n2)
5171        end do
5172      end if
5173 
5174    end if
5175  end if
5176 
5177 end subroutine multpot

m_fftcore/ngfft_seq [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 ngfft_seq

FUNCTION

 Helper function used to initialize ngfft(18) from the FFT divisions
 in the case of sequential execution.

INPUTS

  n123(3)=FFT divisions.

OUTPUT

  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft.

PARENTS

      m_dvdb,m_phgamma,m_wfk,mrgdv

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

420 pure subroutine ngfft_seq(ngfft, n123)
421 
422 
423 !This section has been created automatically by the script Abilint (TD).
424 !Do not modify the following lines by hand.
425 #undef ABI_FUNC
426 #define ABI_FUNC 'ngfft_seq'
427 !End of the abilint section
428 
429  implicit none
430 
431 !Arguments ------------------------------------
432 !scalars
433 !arrays
434  integer,intent(in) :: n123(3)
435  integer,intent(out) :: ngfft(18)
436 
437 !Local variables-------------------------------
438 !scalars
439  integer :: fftalg
440 
441 ! *************************************************************************
442 
443  ! Default  for the sequential case.
444  fftalg = 112
445 #ifdef HAVE_FFT_FFTW3
446  fftalg = 312
447 #elif defined HAVE_FFT_DFTI
448  fftalg = 512
449 #endif
450 
451  ngfft(1:3) = n123
452  ngfft(4) = 2*(ngfft(1)/2)+1
453  ngfft(5) = 2*(ngfft(2)/2)+1
454  ngfft(6) = ngfft(3)
455  ngfft(7)= fftalg           ! fftalg
456  ngfft(8)= get_cache_kb()   ! cache_kb
457  ngfft(9)= 0                ! paral_fft_
458  ngfft(10)=1                ! nproc_fft
459  ngfft(11)=0                ! me_fft
460  ngfft(12)=0                ! n2proc
461  ngfft(13)=0                ! n3proc
462  ngfft(14:18)=0             ! not used
463 
464 end subroutine ngfft_seq

m_fftcore/print_ngfft [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 print_ngfft

FUNCTION

  Print the content of ngfft(18) in explicative format.

INPUTS

  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft.
  [unit]=unit number for output (defaults to std_out).
  [prtvol]=verbosity level (defaults to 0).
  [mode_paral]=either "COLL" or "PERS" ("COLL" is default).

OUTPUT

  Only writing

PARENTS

      bethe_salpeter,eph,getng,m_fft,m_fft_prof,m_wfd,screening,setup_bse
      sigma,wfk_analyze

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

494 subroutine print_ngfft(ngfft,header,unit,mode_paral,prtvol)
495 
496 
497 !This section has been created automatically by the script Abilint (TD).
498 !Do not modify the following lines by hand.
499 #undef ABI_FUNC
500 #define ABI_FUNC 'print_ngfft'
501 !End of the abilint section
502 
503  implicit none
504 
505 !Arguments ------------------------------------
506 !scalars
507  integer,intent(in),optional :: prtvol,unit
508  character(len=*),intent(in),optional :: header
509  character(len=4),intent(in),optional :: mode_paral
510 !arrays
511  integer,intent(in) :: ngfft(18)
512 
513 !Local variables-------------------------------
514 !scalars
515  integer :: my_unt,my_prtvol
516  character(len=4) :: my_mode
517  character(len=500) :: msg
518 
519 ! *************************************************************************
520 
521  my_prtvol=0;       if (PRESENT(prtvol    )) my_prtvol=prtvol
522  my_unt   =std_out; if (PRESENT(unit      )) my_unt   =unit
523  my_mode  ='COLL';  if (PRESENT(mode_paral)) my_mode  =mode_paral
524 
525  msg=ch10//' ==== FFT mesh description (ngfft) ==== '
526  if (PRESENT(header)) msg=ch10//' ==== '//TRIM(ADJUSTL(header))//' ==== '
527  call wrtout(my_unt,msg,my_mode)
528  write(msg,'(2(a,3i5,a),a,i5,2a,i5)')&
529 &  '  FFT mesh divisions ........................ ',ngfft(1),ngfft(2),ngfft(3),ch10,&
530 &  '  Augmented FFT divisions ................... ',ngfft(4),ngfft(5),ngfft(6),ch10,&
531 &  '  FFT algorithm ............................. ',ngfft(7),ch10,&
532 &  '  FFT cache size ............................ ',ngfft(8)
533  call wrtout(my_unt,msg,my_mode)
534 
535  if (my_prtvol>0) then
536    write(msg,'(6(a,i5,a),a,4i5)')&
537 &    '  FFT parallelization level ................. ',ngfft(9),ch10,&
538 &    '  Number of processors in my FFT group ...... ',ngfft(10),ch10,&
539 &    '  Index of me in my FFT group ............... ',ngfft(11),ch10,&
540 &    '  No of xy planes in R space treated by me .. ',ngfft(12),ch10,&
541 &    '  No of xy planes in G space treated by me .. ',ngfft(13),ch10,&
542 &    '  MPI communicator for FFT .................. ',ngfft(14),ch10,&
543 &    '  Value of ngfft(15:18) ..................... ',ngfft(15:18)
544    call wrtout(my_unt,msg,my_mode)
545  end if
546 
547 end subroutine print_ngfft

m_fftcore/scramble [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  scramble

FUNCTION

  This routine performs the local rotation

     input:  G1,R3,G2,(Gp2)
     output: G1,G2,R3,(Gp2)

INPUTS

  i1=Index of x in the small box enclosing the G-sphere.
  j2
  lot=Cache blocking factor
  n1dfft=Number of 1D FFTs performed.
  md1,md2proc,nnd3=Used to dimension zmpi2
  n3=Dimension of the transform along z.
  zw(2,lot,n3): zw(:,1:n1dfft,n3) contains the lines transformed along z

 OUTPTU
 zmpi2(2,md1,md2proc,nnd3)

PARENTS

SOURCE

2537 pure subroutine scramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zw,zmpi2)
2538 
2539 
2540 !This section has been created automatically by the script Abilint (TD).
2541 !Do not modify the following lines by hand.
2542 #undef ABI_FUNC
2543 #define ABI_FUNC 'scramble'
2544 !End of the abilint section
2545 
2546  implicit none
2547 
2548 !Arguments ------------------------------------
2549  integer,intent(in) :: i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3
2550  real(dp),intent(in) :: zw(2,lot,n3)
2551  real(dp),intent(inout) :: zmpi2(2,md1,md2proc,nnd3)
2552 
2553 !Local variables-------------------------------
2554 !scalars
2555  integer :: i3,i
2556 
2557 ! *************************************************************************
2558 
2559  do i3=1,n3
2560    do i=0,n1dfft-1
2561      zmpi2(1,i1+i,j2,i3)=zw(1,i+1,i3)
2562      zmpi2(2,i1+i,j2,i3)=zw(2,i+1,i3)
2563    end do
2564  end do
2565 
2566 end subroutine scramble

m_fftcore/sphere [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 sphere

FUNCTION

 Array cg is defined in sphere with npw points. Insert cg inside box
 of n1*n2*n3 points to define array cfft for fft box.
 corresponds to given element in cg.  rest of cfft is filled with 0 s.

 iflag=1==>insert cg into cfft.
 iflag=2==>insert cg into cfft, where the second and third dimension
 have been switched (needed for new 2002 SGoedecker FFT)
 iflag=-1==> extract cg from cfft.
 iflag=-2==> extract cg from cfft, where the second and third dimension
 have been switched (needed for new 2002 SGoedecker FFT)
  (WARNING : iflag=-2 cannot use symmetry operations)

 There is also the possibility to apply a symmetry operation,
 as well as to make a shift in reciprocal space, or to multiply
 by a constant factor, in the case iflag=-1.
 Multiplication by a constant factor is also possible in the case iflag=-2.

INPUTS

 cg(2,npw*ndat)= contains values for npw G vectors in basis sphere
 ndat=number of FFT to do in //
 npw=number of G vectors in basis at this k point
 cfft(2,n4,n5,n6) = fft box
 n1,n2,n3=physical dimension of the box (cfft)
 n4,n5,n6=memory dimension of cfft
 kg_k(3,npw)=integer coordinates of G vectors in basis sphere
 istwf_k=option parameter that describes the storage of wfs
 iflag=option parameter. Possible values: -1, -2, 1, 2
 me_g0=1 if this node has G=0.
 shiftg(3)=The shift in reciprocal space.
 symm(3,3)=symmetry operation in reciprocal space to be applied.
 xnorm=Normalization factor.

SIDE EFFECTS

 Input/Output
 iflag=1 and 2, insert cg(input) into cfft(output)
 iflag=-1 and -2, extract cg(output) from cfft(input)

NOTES

 cg and cfft are assumed to be of type COMPLEX, although this routine treats
 them as real of twice the length to avoid nonstandard complex*16.
 If istwf_k differs from 1, then special storage modes must be taken
 into account, for symmetric wavefunctions coming from k=(0 0 0) or other
 special k points.

TODO

 1) Order arguments
 2) Split the two cases: from and to sphere (merge with cg_box2gpsh and cg_gsph2box?)
 3) If symmetries are used with or without shiftg, it might happen that the FFT mesh
    is not large enough to accomodate the rotated G, in this case one should return ierr /= 0

PARENTS

      cg_rotate,fourwf,m_fft,m_fftcore,m_fftw3,m_io_kss,wfconv

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

1540 subroutine sphere(cg,ndat,npw,cfft,n1,n2,n3,n4,n5,n6,kg_k,istwf_k,iflag,me_g0,shiftg,symm,xnorm)
1541 
1542 
1543 !This section has been created automatically by the script Abilint (TD).
1544 !Do not modify the following lines by hand.
1545 #undef ABI_FUNC
1546 #define ABI_FUNC 'sphere'
1547 !End of the abilint section
1548 
1549  implicit none
1550 
1551 !Arguments ------------------------------------
1552 !scalars
1553  integer,intent(in) :: iflag,istwf_k,n1,n2,n3,n4,n5,n6,ndat,npw,me_g0
1554  real(dp),intent(in) :: xnorm
1555 !arrays
1556  integer,intent(in) :: kg_k(3,npw),shiftg(3),symm(3,3)
1557  real(dp),intent(inout) :: cfft(2,n4,n5,n6*ndat),cg(2,npw*ndat)
1558 
1559 !Local variables-------------------------------
1560 !scalars
1561  integer :: i1,i1inv,i2,i2inv,i3,i3inv,id1,id2,id3,idat,ipw
1562  integer :: j1,j2,j3,l1,l2,l3,npwmin,use_symmetry,i3dat,i3invdat,i2invdat,ipwdat,i2dat
1563  character(len=500) :: msg
1564 !arrays
1565  integer :: identity(3,3)
1566  integer :: i1inver(n1),i2inver(n2),i3inver(n3)
1567 
1568 ! *************************************************************************
1569 
1570  DBG_ENTER("COLL")
1571 
1572  !In the case of special k-points, invariant under time-reversal,
1573  !but not Gamma, initialize the inverse coordinates. !Remember indeed that
1574  !
1575  !  u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^* and therefore:
1576  !  u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
1577 
1578  if (istwf_k>=2) then
1579    if(istwf_k==2 .or. istwf_k==4 .or. istwf_k==6 .or. istwf_k==8)then
1580      i1inver(1)=1
1581      do i1=2,n1
1582        i1inver(i1)=n1+2-i1
1583      end do
1584    else
1585      do i1=1,n1
1586        i1inver(i1)=n1+1-i1
1587      end do
1588    end if
1589    if(istwf_k>=2 .and. istwf_k<=5)then
1590      i2inver(1)=1
1591      do i2=2,n2
1592        i2inver(i2)=n2+2-i2
1593      end do
1594    else
1595      do i2=1,n2
1596        i2inver(i2)=n2+1-i2
1597      end do
1598    end if
1599    if(istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7)then
1600      i3inver(1)=1
1601      do i3=2,n3
1602        i3inver(i3)=n3+2-i3
1603      end do
1604    else
1605      do i3=1,n3
1606        i3inver(i3)=n3+1-i3
1607      end do
1608    end if
1609  end if
1610 
1611  if (iflag==1 .or. iflag==2) then
1612    ! Insert cg into cfft with extra 0 s around outside:
1613    cfft = zero
1614 
1615    ! Take care of each plane wave, and complete cfft if needed
1616    if (istwf_k==1) then
1617 
1618      if (iflag==1) then
1619 !$OMP PARALLEL DO PRIVATE(i1,i2,i3) IF (ndat>1)
1620        do idat=1,ndat
1621          do ipw=1,npw
1622            i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1623            i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1624            i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1625 
1626            cfft(1,i1,i2,i3+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1627            cfft(2,i1,i2,i3+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1628          end do
1629        end do
1630      end if
1631 
1632      if (iflag==2) then
1633 !$OMP PARALLEL DO PRIVATE(i1,i2,i3) IF (ndat>1)
1634        do idat=1,ndat
1635          do ipw=1,npw
1636            i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1637            i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1638            i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1639 
1640            cfft(1,i1,i3,i2+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1641            cfft(2,i1,i3,i2+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1642          end do
1643        end do
1644      end if
1645 
1646    else if (istwf_k>=2) then
1647 
1648      npwmin=1
1649      if (istwf_k==2 .and. me_g0==1) then
1650        ! If gamma point, then cfft must be completed
1651        do idat=1,ndat
1652          cfft(1,1,1,1+n6*(idat-1))=cg(1,1+npw*(idat-1))
1653          cfft(2,1,1,1+n6*(idat-1))=zero
1654        end do
1655        npwmin=2
1656      end if
1657 
1658      if (iflag==1) then
1659 !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv) IF (ndat>1)
1660        do idat=1,ndat
1661          do ipw=npwmin,npw
1662            i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1663            i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1664            i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1665            ! Construct the coordinates of -k-G
1666            i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1667 
1668            cfft(1,i1,i2,i3+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1669            cfft(2,i1,i2,i3+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1670            cfft(1,i1inv,i2inv,i3inv+n6*(idat-1))= cg(1,ipw+npw*(idat-1))
1671            cfft(2,i1inv,i2inv,i3inv+n6*(idat-1))=-cg(2,ipw+npw*(idat-1))
1672          end do
1673        end do
1674      end if
1675 
1676      if (iflag==2) then
1677 !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv) IF (ndat>1)
1678        do idat=1,ndat
1679          do ipw=npwmin,npw
1680            i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1681            i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1682            i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1683 
1684            ! Construct the coordinates of -k-G
1685            i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1686 
1687            cfft(1,i1,i3,i2+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1688            cfft(2,i1,i3,i2+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1689            cfft(1,i1inv,i3inv,i2inv+n6*(idat-1))= cg(1,ipw+npw*(idat-1))
1690            cfft(2,i1inv,i3inv,i2inv+n6*(idat-1))=-cg(2,ipw+npw*(idat-1))
1691          end do
1692        end do
1693      end if
1694 
1695    end if
1696 
1697  else if (iflag==-1 .or. iflag==-2) then
1698 
1699    use_symmetry=0
1700    identity(:,:)=0; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1701    if(sum((symm(:,:)-identity(:,:))**2)/=0)use_symmetry=1
1702    if(sum(shiftg(:)**2)/=0)use_symmetry=1
1703 
1704    ! Extract cg from cfft, ignoring components outside range of cg:
1705    if (istwf_k==1) then
1706 
1707      if (use_symmetry==0) then
1708        if (iflag==-1) then
1709 !$OMP PARALLEL DO PRIVATE(i1,i2,i3,ipwdat,i3dat) IF (ndat>1)
1710          do idat=1,ndat
1711            do ipw=1,npw
1712              i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1713              i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1714              i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1715              ipwdat = ipw + (idat-1) * npw
1716              i3dat = i3 + (idat-1) * n6
1717 
1718              cg(1,ipwdat)=cfft(1,i1,i2,i3dat)*xnorm
1719              cg(2,ipwdat)=cfft(2,i1,i2,i3dat)*xnorm
1720            end do
1721          end do
1722        else
1723 !$OMP PARALLEL DO PRIVATE(i1,i2,i3,ipwdat,i2dat) IF (ndat>1)
1724          do idat=1,ndat
1725            do ipw=1,npw
1726              i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1727              i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1728              i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1729 
1730              ipwdat = ipw + (idat-1) * npw
1731              i2dat = i2 + (idat-1) * n6
1732 
1733              cg(1,ipwdat)=cfft(1,i1,i3,i2dat)*xnorm
1734              cg(2,ipwdat)=cfft(2,i1,i3,i2dat)*xnorm
1735            end do
1736          end do
1737        end if
1738      else
1739 !$OMP PARALLEL DO PRIVATE(i1,i2,i3,j1,j2,j3,l1,l2,l3,ipwdat,i3dat) IF (ndat>1)
1740        do idat=1,ndat
1741          do ipw=1,npw
1742            l1=kg_k(1,ipw)+shiftg(1)
1743            l2=kg_k(2,ipw)+shiftg(2)
1744            l3=kg_k(3,ipw)+shiftg(3)
1745            j1=symm(1,1)*l1+symm(1,2)*l2+symm(1,3)*l3
1746            j2=symm(2,1)*l1+symm(2,2)*l2+symm(2,3)*l3
1747            j3=symm(3,1)*l1+symm(3,2)*l2+symm(3,3)*l3
1748            if(j1<0)j1=j1+n1; i1=j1+1
1749            if(j2<0)j2=j2+n2; i2=j2+1
1750            if(j3<0)j3=j3+n3; i3=j3+1
1751 
1752            ipwdat = ipw + (idat-1) * npw
1753            i3dat = i3 + (idat-1)*n6
1754 
1755            cg(1,ipwdat)=cfft(1,i1,i2,i3dat)*xnorm
1756            cg(2,ipwdat)=cfft(2,i1,i2,i3dat)*xnorm
1757          end do
1758        end do
1759      end if
1760 
1761    else if (istwf_k>=2) then
1762 
1763      npwmin=1
1764      if (istwf_k==2 .and. me_g0==1) then
1765        ! Extract cg from cfft, in a way that projects on a
1766        ! wavefunction with time-reversal symmetry
1767        do idat=1,ndat
1768          ipwdat = 1 + (idat-1) * npw
1769          i3dat = 1 + (idat-1)*n6
1770          cg(1,ipwdat)=cfft(1,1,1,i3dat)*xnorm
1771          cg(2,ipwdat)=zero
1772        end do
1773        npwmin=2
1774      end if
1775 
1776      if (use_symmetry==0) then
1777 
1778        if (iflag==-1) then
1779 !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv,ipwdat,i3dat,i3invdat) IF (ndat>1)
1780          do idat=1,ndat
1781            do ipw=npwmin,npw
1782              i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1783              i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1784              i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1785 
1786              ! Construct the coordinates of -k-G
1787              i1inv=i1inver(i1); i2inv=i2inver(i2); i3inv=i3inver(i3)
1788 
1789              ipwdat = ipw + (idat-1) * npw
1790              i3dat = i3 + (idat-1) * n6
1791              i3invdat = i3inv + (idat-1) * n6
1792 
1793              ! Here the time-reversal symmetry is used to project from cfft
1794              cg(1,ipwdat)=(cfft(1,i1,i2,i3dat) + cfft(1,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1795              cg(2,ipwdat)=(cfft(2,i1,i2,i3dat) - cfft(2,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1796            end do
1797          end do
1798 
1799        else
1800 !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv,ipwdat,i2dat,i2invdat) IF (ndat>1)
1801          do idat=1,ndat
1802            do ipw=npwmin,npw
1803              i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1804              i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1805              i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1806 
1807              ! Construct the coordinates of -k-G
1808              i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1809 
1810              ipwdat = ipw + (idat-1) * npw
1811              i2dat = i2 + (idat-1) * n6
1812              i2invdat = i2inv + (idat-1) * n6
1813 
1814              ! Here the time-reversal symmetry is used to project from cfft
1815              cg(1,ipwdat)=(cfft(1,i1,i3,i2dat) + cfft(1,i1inv,i3inv,i2invdat))*0.5d0*xnorm
1816              cg(2,ipwdat)=(cfft(2,i1,i3,i2dat) - cfft(2,i1inv,i3inv,i2invdat))*0.5d0*xnorm
1817            end do
1818          end do
1819        end if
1820 
1821      else ! Use symmetry
1822        id1=n1/2+2
1823        id2=n2/2+2
1824        id3=n3/2+2
1825 
1826 !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv,j1,j2,j3,l1,l2,l3,ipwdat,i3dat,i3invdat) IF (ndat>1)
1827        do idat=1,ndat
1828          do ipw=npwmin,npw
1829 
1830            i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1831            i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1832            i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1833 
1834            i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1835 
1836            l1=kg_k(1,ipw)+shiftg(1)
1837            l2=kg_k(2,ipw)+shiftg(2)
1838            l3=kg_k(3,ipw)+shiftg(3)
1839            j1=symm(1,1)*l1+symm(1,2)*l2+symm(1,3)*l3
1840            j2=symm(2,1)*l1+symm(2,2)*l2+symm(2,3)*l3
1841            j3=symm(3,1)*l1+symm(3,2)*l2+symm(3,3)*l3
1842            if(j1<0)j1=j1+n1 ; i1=j1+1
1843            if(j2<0)j2=j2+n2 ; i2=j2+1
1844            if(j3<0)j3=j3+n3 ; i3=j3+1
1845 
1846            ! Construct the coordinates of -k-G
1847            l1=i1inv-(i1inv/id1)*n1-1+shiftg(1)
1848            l2=i2inv-(i2inv/id2)*n2-1+shiftg(2)
1849            l3=i3inv-(i3inv/id3)*n3-1+shiftg(3)
1850            j1=symm(1,1)*l1+symm(1,2)*l2+symm(1,3)*l3
1851            j2=symm(2,1)*l1+symm(2,2)*l2+symm(2,3)*l3
1852            j3=symm(3,1)*l1+symm(3,2)*l2+symm(3,3)*l3
1853            if(j1<0)j1=j1+n1 ; i1inv=j1+1
1854            if(j2<0)j2=j2+n2 ; i2inv=j2+1
1855            if(j3<0)j3=j3+n3 ; i3inv=j3+1
1856 
1857            ipwdat = ipw + (idat-1) * npw
1858            i3dat = i3 + (idat-1) * n6
1859            i3invdat = i3inv + (idat-1) * n6
1860 
1861            ! Here the time-reversal symmetry is used to project from cfft
1862            cg(1,ipwdat)=(cfft(1,i1,i2,i3dat) + cfft(1,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1863            cg(2,ipwdat)=(cfft(2,i1,i2,i3dat) - cfft(2,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1864          end do
1865        end do
1866      end if
1867 
1868    end if
1869 
1870  else
1871    write(msg,'(a,i0,a)')'  iflag=',iflag,' not acceptable.'
1872    MSG_BUG(msg)
1873  end if
1874 
1875  DBG_EXIT("COLL")
1876 
1877 end subroutine sphere

m_fftcore/sphere_fft [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 sphere_fft

FUNCTION

 Array cg is defined in sphere with npw points. Insert cg inside box
 of n1*n2*n3 points to define array cfft for fft box.
 corresponds to given element in cg.  rest of cfft is filled with 0 s.

 iflag=1==>insert cg into cfft.
 iflag=2==>insert cg into cfft, where the second and third dimension
 have been switched (needed for new 2002 SGoedecker FFT)
 iflag=-1==> extract cg from cfft.
 iflag=-2==> extract cg from cfft, where the second and third dimension
 have been switched (needed for new 2002 SGoedecker FFT)
  (WARNING : iflag=-2 cannot use symmetry operations)

 There is also the possibility to apply a symmetry operation,
 as well as to make a shift in reciprocal space, or to multiply
 by a constant factor, in the case iflag=-1.
 Multiplication by a constant factor is also possible in the case iflag=-2.

INPUTS

 cg(2,npw)= contains values for npw G vectors in basis sphere
 ndat=number of FFT to do in //
 npw=number of G vectors in basis at this k point
 cfft(2,n4,n5,n6) = fft box
 n1,n2,n3=physical dimension of the box (cfft)
 n4,n5,n6=memory dimension of cfft
 kg_k(3,npw)=integer coordinates of G vectors in basis sphere
 mpi_enreg=information about MPI parallelization
 tab_fftwf2_local(n2)=local i2 indices in fourwf
 nd2proc TO BE DESCRIBED SB 090831
 iflag=option parameter. Possible values: -1, -2, 1, 2 ; this is used only in debug option

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output
 iflag=1 and 2, insert cg(input) into cfft(output)
 iflag=-1 and -2, extract cg(output) from cfft(input)

NOTES

 cg and cfft are assumed to be of type COMPLEX, although this routine treats
 them as real of twice the length to avoid nonstandard complex*16.

 WARNING
 NO CHECK is DONE over iflag.

TODO

 Order arguments

PARENTS

      fourwf

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

1943 subroutine sphere_fft(cg,ndat,npw,cfft,n1,n2,n3,n4,n5,kg_k,tab_fftwf2_local,nd2proc)
1944 
1945 
1946 !This section has been created automatically by the script Abilint (TD).
1947 !Do not modify the following lines by hand.
1948 #undef ABI_FUNC
1949 #define ABI_FUNC 'sphere_fft'
1950 !End of the abilint section
1951 
1952  implicit none
1953 
1954 !Arguments ------------------------------------
1955 !scalars
1956  integer,intent(in) :: n1,n2,n3,n4,n5,nd2proc,ndat,npw
1957  integer,intent(in) :: tab_fftwf2_local(n2)
1958 !arrays
1959  integer,intent(in) :: kg_k(3,npw)
1960  real(dp),intent(in) :: cg(2,npw*ndat)
1961  real(dp),intent(out) :: cfft(2,n4,n5,nd2proc*ndat)
1962 
1963 !Local variables-------------------------------
1964 !scalars
1965  integer :: i1,i2,i2_local,i3,idat,ipw
1966 
1967 ! *************************************************************************
1968 
1969 !Insert cg into cfft with extra 0 s around outside:
1970  cfft = zero
1971 
1972 !$OMP PARALLEL DO PRIVATE(i1,i2,i2_local,i3)
1973  do ipw=1,npw
1974    i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1975    i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1976    i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1977    i2_local = tab_fftwf2_local(i2)
1978    do idat=1,ndat
1979      cfft(1,i1,i3,i2_local + nd2proc*(idat-1))=cg(1,ipw+npw*(idat-1))
1980      cfft(2,i1,i3,i2_local + nd2proc*(idat-1))=cg(2,ipw+npw*(idat-1))
1981    end do
1982  end do
1983 
1984 end subroutine sphere_fft

m_fftcore/sphere_fft1 [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 sphere_fft1

FUNCTION

 Array cg is defined in sphere with npw points. Insert cg inside box
 of n1*n2*n3 points to define array cfft for fft box.
 corresponds to given element in cg.  rest of cfft is filled with 0 s.

 iflag=1==>insert cg into cfft.
 iflag=2==>insert cg into cfft, where the second and third dimension
 have been switched (needed for new 2002 SGoedecker FFT)
 iflag=-1==> extract cg from cfft.
 iflag=-2==> extract cg from cfft, where the second and third dimension
 have been switched (needed for new 2002 SGoedecker FFT)
  (WARNING : iflag=-2 cannot use symmetry operations)

 There is also the possibility to apply a symmetry operation,
 as well as to make a shift in reciprocal space, or to multiply
 by a constant factor, in the case iflag=-1.
 Multiplication by a constant factor is also possible in the case iflag=-2.

INPUTS

 cg(2,npw)= contains values for npw G vectors in basis sphere
 ndat=number of FFT to do in //
 npw=number of G vectors in basis at this k point
 cfft(2,n4,n5,n6) = fft box
 n1,n2,n3=physical dimension of the box (cfft)
 n4,n5,n6=memory dimension of cfft
 kg_k(3,npw)=integer coordinates of G vectors in basis sphere
 nd2proc TO BE DESCRIBED SB 090831
 iflag=option parameter. Possible values: -1, -2, 1, 2
 tab_fftwf2_local(n2)=local i2 indices in fourwf

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output
 iflag=1 and 2, insert cg(input) into cfft(output)
 iflag=-1 and -2, extract cg(output) from cfft(input)

NOTES

 cg and cfft are assumed to be of type COMPLEX, although this routine treats
 them as real of twice the length to avoid nonstandard complex*16.

 WARNING
 NO CHECK is DONE over iflag.

TODO

   Order arguments
   sphere_fft1 is similar to sphere_fft, the only difference being that ndat > 1 is not supported.
   Why? Should merge the two APIs.

PARENTS

      m_fft

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

2051 subroutine sphere_fft1(cg,ndat,npw,cfft,n1,n2,n3,n4,n5,n6,kg_k,tab_fftwf2_local)
2052 
2053 
2054 !This section has been created automatically by the script Abilint (TD).
2055 !Do not modify the following lines by hand.
2056 #undef ABI_FUNC
2057 #define ABI_FUNC 'sphere_fft1'
2058 !End of the abilint section
2059 
2060  implicit none
2061 
2062 !Arguments ------------------------------------
2063 !scalars
2064  integer,intent(in) :: n1,n2,n3,n4,n5,n6,ndat,npw
2065 !arrays
2066  integer,intent(in) :: kg_k(3,npw)
2067  integer,intent(in) :: tab_fftwf2_local(n2)
2068  real(dp),intent(in) :: cg(2,npw*ndat)
2069  real(dp),intent(inout) :: cfft(2,n4,n5,n6*ndat)
2070 
2071 !Local variables-------------------------------
2072 !scalars
2073  integer :: i1,i2,i2_local,i3,idat,ipw
2074 
2075 ! *************************************************************************
2076 
2077 !Insert cg into cfft with extra 0 s around outside:
2078 
2079  cfft = zero
2080 !$OMP PARALLEL DO PRIVATE(i1,i2,i2_local,i3)
2081  do idat=1,ndat
2082    do ipw=1,npw
2083      i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
2084      i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
2085      i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
2086      i2_local = tab_fftwf2_local(i2) + n6*(idat-1)
2087      cfft(1,i1,i3,i2_local)=cg(1,ipw+npw*(idat-1))
2088      cfft(2,i1,i3,i2_local)=cg(2,ipw+npw*(idat-1))
2089    end do
2090  end do
2091 
2092 end subroutine sphere_fft1

m_fftcore/sphereboundary [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 sphereboundary

FUNCTION

 Finds the boundary of the basis sphere of G vectors (at a given
 k point) for use in improved zero padding of ffts in 3 dimensions.
 Provides data to be used by subroutine fourwf, in the form of
 an array gbound(2*mgfft+8,2).

 The first component (for use when mod(fftalg,10)==2))
 provides integer values g1min,g1max,g2min,g2max
 and then for g2 in the
 sequence g2=0,1,2,...,g2max,g2min,g2min+1,...,-1, provides g1min, g1max.

 The second component (for use when mod(fftalg,10)==1))
 provides integer values g1min,g1max,g3min,g3max,
 where g3min and g3max have been corrected in case of time-reversal
 and then for g3 in the sequence
 g3=0,1,2,...,g3max,g3min,g3min+1,...,-1, provides g2min, g2max.
 (also corrected in case of time-reversal)

 These values are stored in the above order in array gbound.
 Debug mode, if fftalg is between 000 and 099

INPUTS

  istwf_k=option parameter that describes the storage of wfs
  kg_k(3,npw)=integer coordinates of G vectors in basis sphere
  mgfft=maximum size of 1D FFTs (only for dimensioning purposes)
  npw=number of G vectors in basis at this k point

OUTPUT

  gbound(2*mgfft+8,2)=defined above

PARENTS

      dfpt_eltfrkin,dfpt_mkrho,fock_getghc,m_bandfft_kpt,m_cut3d,m_epjdos
      m_fft,m_fft_prof,m_fock,m_gsphere,m_hamiltonian,m_wfd,mkrho,mlwfovlp
      pawmkaewf,posdoppler,scfcv,spin_current,suscep_stat,susk,tddft,wfconv

CHILDREN

SOURCE

1244 subroutine sphereboundary(gbound,istwf_k,kg_k,mgfft,npw)
1245 
1246 
1247 !This section has been created automatically by the script Abilint (TD).
1248 !Do not modify the following lines by hand.
1249 #undef ABI_FUNC
1250 #define ABI_FUNC 'sphereboundary'
1251 !End of the abilint section
1252 
1253  implicit none
1254 
1255 !Arguments ------------------------------------
1256 !scalars
1257  integer,intent(in) :: istwf_k,mgfft,npw
1258 !arrays
1259  integer,intent(in) :: kg_k(3,npw)
1260  integer,intent(out) :: gbound(2*mgfft+8,2)
1261 
1262 !Local variables-------------------------------
1263 !scalars
1264  integer :: dim_a,dim_b,fftalgc,g_a,gmax_a,gmax_b,gmax_b1,gmax_b2,gmin_a,gmin_b
1265  integer :: gmin_b1,gmin_b2,igb,ii,iloop,ipw,testm,testp,kgk
1266  character(len=500) :: message
1267 !arrays
1268  integer :: gmax(2),gmin(2)
1269 
1270 ! *************************************************************************
1271 !
1272 !DEBUG
1273 !write(std_out,*)' sphereboundary : enter'
1274 !write(std_out, '(a)' )' sphereboundary : list of plane waves coordinates for k point '
1275 !write(std_out, '(a)' )'       ipw       kg_k(1:3,ipw) '
1276 !do ipw=1,npw
1277 !write(std_out, '(i10,a,3i6)' )ipw,'  ',kg_k(1:3,ipw)
1278 !end do
1279 !gbound=-999
1280 !ENDDEBUG
1281 
1282 !Determine cube boundaries
1283  gbound(1,1)=minval(kg_k(1,:))
1284  gbound(2,1)=maxval(kg_k(1,:))
1285  gbound(1:2,2)=gbound(1:2,1)
1286 
1287 !Treat differently the fftalgc cases
1288  do ii=1,2
1289 
1290    fftalgc=3-ii
1291 
1292    if(fftalgc/=2)then
1293      dim_a=3
1294      dim_b=2
1295    else
1296      dim_a=2
1297      dim_b=1
1298    end if
1299 
1300 !  Relevant boundaries
1301    gbound(3,ii)=minval(kg_k(dim_a,:))
1302    gbound(4,ii)=maxval(kg_k(dim_a,:))
1303    gmin_a=gbound(3,ii)
1304    gmax_a=gbound(4,ii)
1305 
1306 !  Must complete the sphere for fftalgc==1 and special storage modes.
1307 !  Explanation : sg_fftpad is not able to take into account
1308 !  the time-reversal symmetry, so that the boundaries will not be delimited
1309 !  by the kg_k set, but by their symmetric also.
1310    if(istwf_k>=2 .and. fftalgc==1)then
1311      if( istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7 )then
1312        gbound(4,2)=max(gmax_a,-gmin_a)
1313        gbound(3,2)=-gbound(4,2)
1314      else if( istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9 )then
1315        gbound(4,2)=max(gmax_a,-gmin_a-1)
1316        gbound(3,2)=-gbound(4,2)-1
1317      end if
1318      gmax_a=gbound(4,2) ; gmin_a=gbound(3,2)
1319    end if
1320 
1321    igb=5
1322 
1323 !  Consider first every g_a in range 0 ... gmax_a, then gmin_a ... -1
1324    gmin(1)=0         ; gmax(1)=gmax_a
1325    gmin(2)=gmin_a    ; gmax(2)=-1
1326 
1327    do iloop=1,2
1328 
1329      if( gmin(iloop) <= gmax(iloop) )then
1330 
1331        do g_a=gmin(iloop),gmax(iloop)
1332 
1333          if(istwf_k==1 .or. fftalgc/=1)then
1334            ! Select the minimal and maximal values, in the selected plane
1335            gmin_b=mgfft+1 ! Initialized with a value larger than all possible ones
1336            gmax_b=-mgfft-1 ! Initialized with a value smaller than all possible ones
1337            do ipw=1,npw
1338              if(kg_k(dim_a,ipw)==g_a)then
1339                kgk=kg_k(dim_b,ipw)
1340                if(kgk<=gmin_b)gmin_b=kgk
1341                if(kgk>=gmax_b)gmax_b=kgk
1342              end if
1343            end do
1344 
1345          else if(istwf_k>=2 .and. fftalgc==1)then
1346 
1347            ! Here, must take into account time-reversal symmetry explicitely
1348 
1349            ! Determine the boundaries for the plane g_a
1350            testp=0
1351            if(g_a<=gmax_a)then
1352              ! Select the minimal and maximal values, in the selected plane
1353              gmin_b1=mgfft+1 ! Initialized with a value larger than all possible ones
1354              gmax_b1=-mgfft-1 ! Initialized with a value smaller than all possible ones
1355              do ipw=1,npw
1356                if(kg_k(dim_a,ipw)==g_a)then
1357                  kgk=kg_k(dim_b,ipw)
1358                  if(kgk<=gmin_b1)gmin_b1=kgk
1359                  if(kgk>=gmax_b1)gmax_b1=kgk
1360                end if
1361              end do
1362 
1363 
1364              testp=1
1365            end if
1366 
1367            ! Determine the boundaries for the plane -g_a or -g_a-1
1368            testm=0
1369            if( istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7 )then
1370 
1371              if(-g_a>=gmin_a)then
1372                ! Select the minimal and maximal values, in the selected plane
1373                ! Warning : there is an inversion of search (might be confusing)
1374                gmax_b2=mgfft+1 ! Initialized with a value larger than all possible ones
1375                gmin_b2=-mgfft-1 ! Initialized with a value smaller than all possible ones
1376                do ipw=1,npw
1377                  if(kg_k(dim_a,ipw)==-g_a)then
1378                    kgk=kg_k(dim_b,ipw)
1379                    if(kgk<=gmax_b2)gmax_b2=kgk
1380                    if(kgk>=gmin_b2)gmin_b2=kgk
1381                  end if
1382                end do
1383                testm=1
1384              end if
1385 
1386            else if( istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9 )then
1387 
1388              if(-g_a-1>=gmin_a)then
1389                ! Select the minimal and maximal values, in the selected plane
1390                ! Warning : there is an inversion of search (might be confusing)
1391                gmax_b2=mgfft+1 ! Initialized with a value larger than all possible ones
1392                gmin_b2=-mgfft-1 ! Initialized with a value smaller than all possible ones
1393                do ipw=1,npw
1394                  if(kg_k(dim_a,ipw)==-g_a-1)then
1395                    kgk=kg_k(dim_b,ipw)
1396                    if(kgk<=gmax_b2)gmax_b2=kgk
1397                    if(kgk>=gmin_b2)gmin_b2=kgk
1398                  end if
1399                end do
1400                testm=1
1401              end if
1402 
1403            end if
1404 
1405            !  Must invert the boundaries, to use them for plane g_a
1406            if(testm==1)then
1407              ! This is needed to avoid border effect
1408              ! if the search did not lead to any element
1409              gmin_b2=max(gmin_b2,-mgfft) ; gmax_b2=min(gmax_b2,mgfft)
1410              if(istwf_k<=5)then
1411                gmax_b2=-gmax_b2 ; gmin_b2=-gmin_b2
1412              else
1413                gmax_b2=-gmax_b2-1 ; gmin_b2=-gmin_b2-1
1414              end if
1415            end if
1416 
1417            if( testp==1 .and. testm==1)then
1418              gmin_b=min(gmin_b1,gmin_b2) ; gmax_b=max(gmax_b1,gmax_b2)
1419            else if( testp==1 )then
1420              gmin_b=gmin_b1 ; gmax_b=gmax_b1
1421            else if( testm==1 )then
1422              gmin_b=gmin_b2 ; gmax_b=gmax_b2
1423            end if
1424 
1425          end if ! Endif take into account time-reversal symmetry
1426 
1427          if (igb+1>2*mgfft+4) then
1428            write(message, '(2a, 4(a,3(i0,1x),a))' )&
1429              "About to overwrite gbound array (FFT mesh too small) ",ch10, &
1430              "   iloop, igb, mgb = ",iloop,igb,2*mgfft+4, ch10, &
1431              "   istwfk, mgfft, npw = ",istwf_k, mgfft, npw, ch10, &
1432              "   minval(kg_k) = ",minval(kg_k, dim=2), ch10, &
1433              "   maxval(kg_k) = ",maxval(kg_k, dim=2), ch10
1434            MSG_BUG(message)
1435          end if
1436 
1437          gbound(igb,ii)=gmin_b
1438          gbound(igb+1,ii)=gmax_b
1439 
1440          if( iloop==1 .and. istwf_k>=2 .and. istwf_k<=5 .and. fftalgc==2 .and. g_a==0)then
1441 !          If k_y=0 , for fftalgc==2, the g_a==0 plane must be completed
1442            if(istwf_k==2 .or. istwf_k==4)then
1443              gbound(igb+1,ii)=max(gmax_b,-gmin_b)
1444              gbound(igb,ii)=-gbound(igb+1,ii)
1445            else if(istwf_k==3 .or. istwf_k==5)then
1446              gbound(igb+1,ii)=max(gmax_b,-gmin_b-1)
1447              gbound(igb,ii)=-gbound(igb+1,ii)-1
1448            end if
1449 
1450          end if
1451 
1452          igb=igb+2
1453 
1454        end do ! g_a
1455      end if
1456    end do  ! iloop
1457  end do ! ii (fftalgc)
1458 
1459 !DEBUG
1460 !write(std_out,'(a)')' sphereoundary : list of plane waves coordinates for 1st k point '
1461 !write(std_out,'(a)')'       ipw       kg_k(1:3,ipw) '
1462 !do ipw=1,npw
1463 !write(std_out, '(i10,a,3i6)' )ipw,'  ',kg_k(1:3,ipw)
1464 !end do
1465 !write(std_out, '(a)' )' sphereboundary : list of boundaries '
1466 !do igb=1,2*mgfft+8
1467 !write(std_out, '(i10,a,2i6)' )igb,'  ',gbound(igb,1),gbound(igb,2)
1468 !end do
1469 !write(std_out,*)' sphereboundary : exit '
1470 !ENDDEBUG
1471 
1472 end subroutine sphereboundary

m_fftcore/switch [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  switch

FUNCTION

INPUTS

OUTPUT

PARENTS

SOURCE

2190 pure subroutine switch(n1dfft,n2,lot,n1,lzt,zt,zw)
2191 
2192 
2193 !This section has been created automatically by the script Abilint (TD).
2194 !Do not modify the following lines by hand.
2195 #undef ABI_FUNC
2196 #define ABI_FUNC 'switch'
2197 !End of the abilint section
2198 
2199  implicit none
2200 
2201 !Arguments ------------------------------------
2202  integer,intent(in) :: n1dfft,n2,lot,n1,lzt
2203  real(dp),intent(in) :: zt(2,lzt,n1)
2204  real(dp),intent(inout) :: zw(2,lot,n2)
2205 
2206 !Local variables-------------------------------
2207  integer :: i,j
2208 ! *************************************************************************
2209 
2210  do j=1,n1dfft
2211    do i=1,n2
2212      zw(1,j,i)=zt(1,i,j)
2213      zw(2,j,i)=zt(2,i,j)
2214    end do
2215  end do
2216 
2217 end subroutine switch

m_fftcore/switch_cent [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 switch_cent

FUNCTION

   Perform the rotation:

     input: I2,i1,j3,(jp3)
     output: i1,I2,j3,(jp3)

   and padd the signal with zeros.

INPUTS

  n1dfft=Number of 1D FFTs to perform
  max2=Max G_y in the small box enclosing the G-sphere.
  m2=Size of the small box enclosing the G-sphere along y
  n2=Dimension of the transform along y
  lot=Cache blocking factor.
  n1=Dimension of the transform along x
  lzt=Second dimension of z
  zt(2,lzt,n1)

OUTPUT

  zw(2,lot,n2)=Cache working array

PARENTS

SOURCE

2251 pure subroutine switch_cent(n1dfft,max2,m2,n2,lot,n1,lzt,zt,zw)
2252 
2253 
2254 !This section has been created automatically by the script Abilint (TD).
2255 !Do not modify the following lines by hand.
2256 #undef ABI_FUNC
2257 #define ABI_FUNC 'switch_cent'
2258 !End of the abilint section
2259 
2260  implicit none
2261 
2262 !Arguments ------------------------------------
2263  integer,intent(in) :: n1dfft,max2,m2,n2,lot,n1,lzt
2264  real(dp),intent(in) :: zt(2,lzt,n1)
2265  real(dp),intent(inout) :: zw(2,lot,n2)
2266 
2267 !Local variables-------------------------------
2268  integer :: i,j
2269 
2270 ! *************************************************************************
2271 
2272  ! Here, zero and positive frequencies
2273  do j=1,n1dfft
2274    do i=1,max2+1
2275      zw(1,j,i)=zt(1,i,j)
2276      zw(2,j,i)=zt(2,i,j)
2277    end do
2278  end do
2279 
2280  ! Fill the center region with zeros
2281  do i=max2+2,n2-m2+max2+1
2282    do j=1,n1dfft
2283      zw(1,j,i)=zero
2284      zw(2,j,i)=zero
2285    end do
2286  end do
2287 
2288  ! Here, negative frequencies
2289  if (m2>=max2+2) then
2290    do j=1,n1dfft
2291      do i=max2+2,m2
2292        zw(1,j,i+n2-m2)=zt(1,i,j)
2293        zw(2,j,i+n2-m2)=zt(2,i,j)
2294      end do
2295    end do
2296  end if
2297 
2298 end subroutine switch_cent

m_fftcore/switchreal [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  switchreal

FUNCTION

   Perform the rotation:

     input: I2,i1,j3,(jp3)
     output: i1,I2,j3,(jp3)

   and padd the signal with zeros.
   Used for real wavefunctions.

INPUTS

  includelast
  n1dfft=Number of 1D FFTs to perform
  n2=Dimension of the transform along y
  n2eff
  lot=Cache blocking factor.
  n1zt
  lzt
  zt(2,lzt,n1zt)

OUTPUT

  zw(2,lot,n2)

PARENTS

SOURCE

2333 pure subroutine switchreal(includelast,n1dfft,n2,n2eff,lot,n1zt,lzt,zt,zw)
2334 
2335 
2336 !This section has been created automatically by the script Abilint (TD).
2337 !Do not modify the following lines by hand.
2338 #undef ABI_FUNC
2339 #define ABI_FUNC 'switchreal'
2340 !End of the abilint section
2341 
2342  implicit none
2343 
2344 !Arguments ------------------------------------
2345  integer,intent(in) :: includelast,n1dfft,n2,n2eff,lot,n1zt,lzt
2346  real(dp),intent(in) :: zt(2,lzt,n1zt)
2347  real(dp),intent(inout) :: zw(2,lot,n2)
2348 
2349 !Local variables-------------------------------
2350  integer :: i,j
2351 ! *************************************************************************
2352 
2353  if (includelast==1) then
2354 
2355    ! Compute symmetric and antisymmetric combinations
2356    do j=1,n1dfft
2357      zw(1,j,1)=zt(1,1,2*j-1)
2358      zw(2,j,1)=zt(1,1,2*j  )
2359    end do
2360    do i=2,n2eff
2361      do j=1,n1dfft
2362        zw(1,j,i)=      zt(1,i,2*j-1)-zt(2,i,2*j)
2363        zw(2,j,i)=      zt(2,i,2*j-1)+zt(1,i,2*j)
2364        zw(1,j,n2+2-i)= zt(1,i,2*j-1)+zt(2,i,2*j)
2365        zw(2,j,n2+2-i)=-zt(2,i,2*j-1)+zt(1,i,2*j)
2366      end do
2367    end do
2368 
2369  else
2370 
2371    ! An odd number of FFTs
2372    ! Compute symmetric and antisymmetric combinations
2373    do j=1,n1dfft-1
2374      zw(1,j,1)=zt(1,1,2*j-1)
2375      zw(2,j,1)=zt(1,1,2*j  )
2376    end do
2377    zw(1,n1dfft,1)=zt(1,1,2*n1dfft-1)
2378    zw(2,n1dfft,1)=zero
2379 
2380    do i=2,n2eff
2381      do j=1,n1dfft-1
2382        zw(1,j,i)=      zt(1,i,2*j-1)-zt(2,i,2*j)
2383        zw(2,j,i)=      zt(2,i,2*j-1)+zt(1,i,2*j)
2384        zw(1,j,n2+2-i)= zt(1,i,2*j-1)+zt(2,i,2*j)
2385        zw(2,j,n2+2-i)=-zt(2,i,2*j-1)+zt(1,i,2*j)
2386      end do
2387      zw(1,n1dfft,i)=      zt(1,i,2*n1dfft-1)
2388      zw(2,n1dfft,i)=      zt(2,i,2*n1dfft-1)
2389      zw(1,n1dfft,n2+2-i)= zt(1,i,2*n1dfft-1)
2390      zw(2,n1dfft,n2+2-i)=-zt(2,i,2*n1dfft-1)
2391    end do
2392  end if
2393 
2394 end subroutine switchreal

m_fftcore/switchreal_cent [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  switchreal_cent

FUNCTION

   Perform the rotation:

     input: I2,i1,j3,(jp3)
     output: i1,I2,j3,(jp3)

   and padd the signal with zeros.
   Used for the Fourier transform of real wavefunctions.

INPUTS

  includelast
  n1dfft=Number of 1D FFTs to perform
  max2
  n2=Dimension of the transform along y
  lot=Cache blocking factor.
  n1zt, lzt=Dimensions of zt.
  zt(2,lzt,n1zt)

OUTPUT

  zw(2,lot,n2)

PARENTS

SOURCE

2428 pure subroutine switchreal_cent(includelast,n1dfft,max2,n2,lot,n1zt,lzt,zt,zw)
2429 
2430 
2431 !This section has been created automatically by the script Abilint (TD).
2432 !Do not modify the following lines by hand.
2433 #undef ABI_FUNC
2434 #define ABI_FUNC 'switchreal_cent'
2435 !End of the abilint section
2436 
2437  implicit none
2438 
2439 !Arguments ------------------------------------
2440  integer,intent(in) :: includelast,n1dfft,max2,n2,lot,n1zt,lzt
2441  real(dp),intent(in) :: zt(2,lzt,n1zt)
2442  real(dp),intent(inout) :: zw(2,lot,n2)
2443 
2444 !Local variables-------------------------------
2445  integer :: i,j
2446 ! *************************************************************************
2447 
2448  if (includelast==1) then
2449 
2450    ! Compute symmetric and antisymmetric combinations
2451    do j=1,n1dfft
2452     zw(1,j,1)=zt(1,1,2*j-1)
2453     zw(2,j,1)=zt(1,1,2*j  )
2454    end do
2455 
2456    do i=2,max2+1
2457     do j=1,n1dfft
2458      zw(1,j,i)=      zt(1,i,2*j-1)-zt(2,i,2*j)
2459      zw(2,j,i)=      zt(2,i,2*j-1)+zt(1,i,2*j)
2460      zw(1,j,n2+2-i)= zt(1,i,2*j-1)+zt(2,i,2*j)
2461      zw(2,j,n2+2-i)=-zt(2,i,2*j-1)+zt(1,i,2*j)
2462     end do
2463    end do
2464 
2465    if(max2+1<n2-max2)then
2466     do i=max2+2,n2-max2
2467      do j=1,n1dfft
2468       zw(1,j,i)=zero
2469       zw(2,j,i)=zero
2470      end do
2471     end do
2472    end if
2473 
2474  else
2475    ! Compute symmetric and antisymmetric combinations
2476    do j=1,n1dfft-1
2477      zw(1,j,1)=zt(1,1,2*j-1)
2478      zw(2,j,1)=zt(1,1,2*j  )
2479    end do
2480 
2481    zw(1,n1dfft,1)=zt(1,1,2*n1dfft-1)
2482    zw(2,n1dfft,1)=zero
2483    do i=2,max2+1
2484      do j=1,n1dfft-1
2485        zw(1,j,i)=      zt(1,i,2*j-1)-zt(2,i,2*j)
2486        zw(2,j,i)=      zt(2,i,2*j-1)+zt(1,i,2*j)
2487        zw(1,j,n2+2-i)= zt(1,i,2*j-1)+zt(2,i,2*j)
2488        zw(2,j,n2+2-i)=-zt(2,i,2*j-1)+zt(1,i,2*j)
2489      end do
2490      zw(1,n1dfft,i)=      zt(1,i,2*n1dfft-1)
2491      zw(2,n1dfft,i)=      zt(2,i,2*n1dfft-1)
2492      zw(1,n1dfft,n2+2-i)= zt(1,i,2*n1dfft-1)
2493      zw(2,n1dfft,n2+2-i)=-zt(2,i,2*n1dfft-1)
2494    end do
2495 
2496    if(max2+1<n2-max2)then
2497      do i=max2+2,n2-max2
2498        do j=1,n1dfft
2499         zw(1,j,i)=zero
2500         zw(2,j,i)=zero
2501        end do
2502      end do
2503    end if
2504  end if
2505 
2506 end subroutine switchreal_cent

m_fftcore/unfill [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  unfill

FUNCTION

  Move data from the cache work array to zf

INPUTS

  nd1,nd3=Dimensions of the input array zf.
  lot=Cache blocking factor.
  n1dfft=Number of 1D FFTs to perform
  n3=Dimension of the transform along z
  zw(2,lot,n3)=Cache work array with the z-lines.

OUTPUT

  zf(2,nd1,nd3)= zf(:,1:n1dfft,:1:n3) is filled with the results stored in zw

PARENTS

SOURCE

2724 pure subroutine unfill(nd1,nd3,lot,n1dfft,n3,zw,zf)
2725 
2726 
2727 !This section has been created automatically by the script Abilint (TD).
2728 !Do not modify the following lines by hand.
2729 #undef ABI_FUNC
2730 #define ABI_FUNC 'unfill'
2731 !End of the abilint section
2732 
2733  implicit none
2734 
2735 !Arguments ------------------------------------
2736  integer,intent(in) :: nd1,nd3,lot,n1dfft,n3
2737  real(dp),intent(in) :: zw(2,lot,n3)
2738  real(dp),intent(inout) :: zf(2,nd1,nd3)
2739 
2740 !Local variables-------------------------------
2741  integer :: i1,i3
2742 ! *************************************************************************
2743 
2744  do i3=1,n3
2745    do i1=1,n1dfft
2746      zf(1,i1,i3)=zw(1,i1,i3)
2747      zf(2,i1,i3)=zw(2,i1,i3)
2748    end do
2749  end do
2750 
2751 end subroutine unfill

m_fftcore/unfill_cent [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  unfill_cent

FUNCTION

  Transfer data from the cache work array to zf.
  Takes into account zero padding (only the non-zero entries are moved)

INPUTS

  md1,md3=Leading dimension of zf along x and z
  lot=Cache blocking factor.
  n1dfft=Number of 1d transforms performed along z.
  max3=Max index of G_z the small box enclosing the G-sphere.
  m3=Number of points in the *small* box enclosing the G-sphere
  n3=Dimension of the FFT transform along z
  zw(2,lot,n3)=Cache work array

OUTPUT

  zf(2,md1,md3)= zf(:,1:n1dfft,1:m3) is filled with the non-zero components.

PARENTS

SOURCE

2780 pure subroutine unfill_cent(md1,md3,lot,n1dfft,max3,m3,n3,zw,zf)
2781 
2782 
2783 !This section has been created automatically by the script Abilint (TD).
2784 !Do not modify the following lines by hand.
2785 #undef ABI_FUNC
2786 #define ABI_FUNC 'unfill_cent'
2787 !End of the abilint section
2788 
2789  implicit none
2790 
2791 !Arguments ------------------------------------
2792  integer,intent(in) :: md1,md3,lot,n1dfft,max3,m3,n3
2793  real(dp),intent(in) :: zw(2,lot,n3)
2794  real(dp),intent(inout) :: zf(2,md1,md3)
2795 
2796 !Local variables-------------------------------
2797  integer :: i1,i3
2798 ! *************************************************************************
2799 
2800  ! Here, zero and positive frequencies
2801  do i3=1,max3+1
2802    do i1=1,n1dfft
2803      zf(1,i1,i3)=zw(1,i1,i3)
2804      zf(2,i1,i3)=zw(2,i1,i3)
2805    end do
2806  end do
2807 
2808  ! Here, negative frequencies
2809  do i3=max3+2,m3
2810    do i1=1,n1dfft
2811      zf(1,i1,i3)=zw(1,i1,i3+n3-m3)
2812      zf(2,i1,i3)=zw(2,i1,i3+n3-m3)
2813    end do
2814  end do
2815 
2816 end subroutine unfill_cent

m_fftcore/unmpiswitch [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  unmpiswitch

FUNCTION

INPUTS

OUTPUT

PARENTS

      m_fftw3,m_sg2002

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

2839 pure subroutine unmpiswitch(j3,n1dfft,Jp2st,J2st,lot,n1,nd2proc,nd3proc,nproc,ioption,zw,zmpi1)
2840 
2841 
2842 !This section has been created automatically by the script Abilint (TD).
2843 !Do not modify the following lines by hand.
2844 #undef ABI_FUNC
2845 #define ABI_FUNC 'unmpiswitch'
2846 !End of the abilint section
2847 
2848  implicit none
2849 
2850 !Arguments ------------------------------------
2851  integer,intent(in) :: j3,n1dfft,lot,n1,nd2proc,nd3proc,nproc,ioption
2852  integer,intent(inout) :: Jp2st,J2st
2853  real(dp),intent(in) :: zw(2,lot,n1)
2854  real(dp),intent(inout) :: zmpi1(2,n1,nd2proc,nd3proc,nproc)
2855 
2856 !Local variables-------------------------------
2857  integer :: i1,jp2,j2,ind,jjp2,mfft,jj2
2858 
2859 ! *************************************************************************
2860 
2861  mfft=0
2862  if (ioption == 2) then
2863    do Jp2=Jp2st,nproc
2864      do J2=J2st,nd2proc
2865        mfft=mfft+1
2866        if (mfft.gt.n1dfft) then
2867          Jp2st=Jp2
2868          J2st=J2
2869          return
2870        end if
2871        do I1=1,n1
2872          zmpi1(1,I1,J2,j3,Jp2)=zw(1,mfft,I1)
2873          zmpi1(2,I1,J2,j3,Jp2)=zw(2,mfft,I1)
2874        end do
2875      end do
2876      J2st=1
2877    end do
2878 
2879  else
2880    do Jp2=Jp2st,nproc
2881      do J2=J2st,nd2proc
2882        mfft=mfft+1
2883        if (mfft.gt.n1dfft) then
2884          Jp2st=Jp2
2885          J2st=J2
2886          return
2887        end if
2888        ind=(Jp2-1) * nd2proc + J2
2889        jj2=(ind-1)/nproc +1
2890 
2891        !jjp2=modulo(ind,nproc) +1
2892        jjp2=modulo(ind-1,nproc)+1
2893 
2894        do I1=1,n1
2895          zmpi1(1,I1,jj2,j3,jjp2)=zw(1,mfft,I1)
2896          zmpi1(2,I1,jj2,j3,jjp2)=zw(2,mfft,I1)
2897        end do
2898      end do
2899      J2st=1
2900    end do
2901  end if
2902 
2903 end subroutine unmpiswitch

m_fftcore/unmpiswitch_cent [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  unmpiswitch_cent

FUNCTION

INPUTS

OUTPUT

PARENTS

      m_fftw3,m_sg2002

CHILDREN

      xmpi_sum,xmpi_sum_master

SOURCE

2926 pure subroutine unmpiswitch_cent(j3,n1dfft,Jp2stf,J2stf,lot,max1,md1,m1,n1,md2proc,nd3proc,nproc,ioption,zw,zmpi1)
2927 
2928 
2929 !This section has been created automatically by the script Abilint (TD).
2930 !Do not modify the following lines by hand.
2931 #undef ABI_FUNC
2932 #define ABI_FUNC 'unmpiswitch_cent'
2933 !End of the abilint section
2934 
2935  implicit none
2936 
2937 !Arguments ------------------------------------
2938  integer,intent(in) :: j3,n1dfft,lot,max1,md1,m1,n1,md2proc,nd3proc,nproc,ioption
2939  integer,intent(inout) :: Jp2stf,J2stf
2940  real(dp),intent(inout) :: zmpi1(2,md1,md2proc,nd3proc,nproc)
2941  real(dp),intent(in) :: zw(2,lot,n1)
2942 
2943 !Local variables-------------------------------
2944  integer :: mfft,Jp2,J2,I1,ind,jj2,jjp2
2945 
2946 ! *************************************************************************
2947 
2948  mfft=0
2949 
2950  if (ioption == 2) then
2951    do Jp2=Jp2stf,nproc
2952      do J2=J2stf,md2proc
2953        mfft=mfft+1
2954 
2955        if (mfft.gt.n1dfft) then
2956          Jp2stf=Jp2
2957          J2stf=J2
2958          return
2959        end if
2960 
2961        ! Here, zero and positive frequencies
2962        do I1=1,max1+1
2963          zmpi1(1,I1,J2,j3,Jp2)=zw(1,mfft,I1)
2964          zmpi1(2,I1,J2,j3,Jp2)=zw(2,mfft,I1)
2965        end do
2966 
2967        ! Here, negative frequencies
2968        do I1=max1+2,m1
2969          zmpi1(1,I1,J2,j3,Jp2)=zw(1,mfft,I1+n1-m1)
2970          zmpi1(2,I1,J2,j3,Jp2)=zw(2,mfft,I1+n1-m1)
2971        end do
2972 
2973      end do
2974      J2stf=1
2975    end do
2976 
2977  else
2978    do Jp2=Jp2stf,nproc
2979      do J2=J2stf,md2proc
2980        mfft=mfft+1
2981        if (mfft.gt.n1dfft) then
2982          Jp2stf=Jp2
2983          J2stf=J2
2984          return
2985        end if
2986        ind=(Jp2-1) * md2proc + J2
2987        jj2=(ind-1)/nproc +1
2988 
2989        !jjp2=modulo(ind,nproc) +1
2990        jjp2=modulo(ind-1,nproc)+1
2991 
2992        ! Here, zero and positive frequencies
2993        do I1=1,max1+1
2994          zmpi1(1,I1,Jj2,j3,Jjp2)=zw(1,mfft,I1)
2995          zmpi1(2,I1,Jj2,j3,Jjp2)=zw(2,mfft,I1)
2996        end do
2997 
2998        ! Here, negative frequencies
2999        do I1=max1+2,m1
3000          zmpi1(1,I1,Jj2,j3,Jjp2)=zw(1,mfft,I1+n1-m1)
3001          zmpi1(2,I1,Jj2,j3,Jjp2)=zw(2,mfft,I1+n1-m1)
3002        end do
3003      end do
3004      J2stf=1
3005    end do
3006  end if
3007 
3008 end subroutine unmpiswitch_cent

m_fftcore/unscramble [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  unscramble

FUNCTION

INPUTS

  i1
  j2
  lot=Cache blocking factor.
  n1dfft=Number of 1D FFTs to perform
  md1,n3,md2proc,nnd3
  zmpi2(2,md1,md2proc,nnd3)

OUTPUT

  zw(2,lot,n3)= cache work array

PARENTS

SOURCE

3034 pure subroutine unscramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zmpi2,zw)
3035 
3036 
3037 !This section has been created automatically by the script Abilint (TD).
3038 !Do not modify the following lines by hand.
3039 #undef ABI_FUNC
3040 #define ABI_FUNC 'unscramble'
3041 !End of the abilint section
3042 
3043  implicit none
3044 
3045 !Arguments ------------------------------------
3046  integer,intent(in) :: i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3
3047  real(dp),intent(in) :: zmpi2(2,md1,md2proc,nnd3)
3048  real(dp),intent(inout) :: zw(2,lot,n3)
3049 
3050 !Local variables-------------------------------
3051 !scalars
3052  integer :: i,i3
3053 
3054 ! *************************************************************************
3055 
3056  do i3=1,n3
3057    do i=0,n1dfft-1
3058      zw(1,i+1,i3)=zmpi2(1,i1+i,j2,i3)
3059      zw(2,i+1,i3)=zmpi2(2,i1+i,j2,i3)
3060    end do
3061  end do
3062 
3063 end subroutine unscramble

m_fftcore/unswitch [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  unswitch

FUNCTION

INPUTS

  n1dfft=Number of 1D FFTs
  n2=Dimension of the transform along y
  lot=Cache blocking factor.
  n1=Dimension of the transform along x.
  lzt
  zw(2,lot,n2)=Cache work array

OUTPUT

  zt(2,lzt,n1)

PARENTS

SOURCE

3089 pure subroutine unswitch(n1dfft,n2,lot,n1,lzt,zw,zt)
3090 
3091 
3092 !This section has been created automatically by the script Abilint (TD).
3093 !Do not modify the following lines by hand.
3094 #undef ABI_FUNC
3095 #define ABI_FUNC 'unswitch'
3096 !End of the abilint section
3097 
3098  implicit none
3099 
3100 !Arguments ------------------------------------
3101  integer,intent(in) :: n1dfft,n2,lot,n1,lzt
3102  real(dp),intent(in) :: zw(2,lot,n2)
3103  real(dp),intent(inout) :: zt(2,lzt,n1)
3104 
3105 !Local variables-------------------------------
3106  integer :: i,j
3107 ! *************************************************************************
3108 
3109  do j=1,n1dfft
3110    do i=1,n2
3111      zt(1,i,j)=zw(1,j,i)
3112      zt(2,i,j)=zw(2,j,i)
3113    end do
3114  end do
3115 
3116 end subroutine unswitch

m_fftcore/unswitch_cent [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  unswitch_cent

FUNCTION

INPUTS

  n1dfft=Number of 1D FFTs to perform
  max2=Max G_y in the small box enclosing the G-sphere.
  m2=Size of the small box enclosing the G-sphere along y
  n2=Dimension of the transform along y
  lot=Cache blocking factor.
  n1=Dimension of the transform along x
  lzt
  zw(2,lot,n2)=Cache working array

OUTPUT

  zt(2,lzt,n1)

PARENTS

SOURCE

3144 pure subroutine unswitch_cent(n1dfft,max2,m2,n2,lot,n1,lzt,zw,zt)
3145 
3146 
3147 !This section has been created automatically by the script Abilint (TD).
3148 !Do not modify the following lines by hand.
3149 #undef ABI_FUNC
3150 #define ABI_FUNC 'unswitch_cent'
3151 !End of the abilint section
3152 
3153  implicit none
3154 
3155 !Arguments ------------------------------------
3156  integer,intent(in) :: n1dfft,max2,m2,n2,lot,n1,lzt
3157  real(dp),intent(in) :: zw(2,lot,n2)
3158  real(dp),intent(inout) :: zt(2,lzt,n1)
3159 
3160 !Local variables-------------------------------
3161  integer :: i,j
3162 ! *************************************************************************
3163 
3164 ! Here, zero and positive frequencies
3165  do j=1,n1dfft
3166    do i=1,max2+1
3167      zt(1,i,j)=zw(1,j,i)
3168      zt(2,i,j)=zw(2,j,i)
3169    end do
3170  end do
3171 
3172 ! Here, negative frequencies
3173  if(m2>=max2+2)then
3174    do j=1,n1dfft
3175      do i=max2+2,m2
3176        zt(1,i,j)=zw(1,j,i+n2-m2)
3177        zt(2,i,j)=zw(2,j,i+n2-m2)
3178      end do
3179    end do
3180  end if
3181 
3182 end subroutine unswitch_cent

m_fftcore/unswitchreal [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

  unswitchreal

FUNCTION

INPUTS

  n1dfft=Number of 1D FFTs to perform
  n2=Dimension of the transform along y
  n2eff=
  lot=Cache blocking factor.
  n1zt
  lzt
  zw(2,lot,n2)=Cache working array

OUTPUT

  zt(2,lzt,n1)

PARENTS

SOURCE

3209 pure subroutine unswitchreal(n1dfft,n2,n2eff,lot,n1zt,lzt,zw,zt)
3210 
3211 
3212 !This section has been created automatically by the script Abilint (TD).
3213 !Do not modify the following lines by hand.
3214 #undef ABI_FUNC
3215 #define ABI_FUNC 'unswitchreal'
3216 !End of the abilint section
3217 
3218  implicit none
3219 
3220 !Arguments ------------------------------------
3221  integer,intent(in) :: n1dfft,n2,n2eff,lot,n1zt,lzt
3222  real(dp),intent(in) :: zw(2,lot,n2)
3223  real(dp),intent(inout) :: zt(2,lzt,n1zt)
3224 
3225 !Local variables-------------------------------
3226  integer :: i,j
3227 ! *************************************************************************
3228 
3229 ! Decompose symmetric and antisymmetric parts
3230  do j=1,n1dfft
3231    zt(1,1,2*j-1)=zw(1,j,1)
3232    zt(2,1,2*j-1)=zero
3233    zt(1,1,2*j)  =zw(2,j,1)
3234    zt(2,1,2*j)  =zero
3235  end do
3236 
3237  do i=2,n2eff
3238    do j=1,n1dfft
3239      zt(1,i,2*j-1)= (zw(1,j,i)+zw(1,j,n2+2-i))*half
3240      zt(2,i,2*j-1)= (zw(2,j,i)-zw(2,j,n2+2-i))*half
3241      zt(1,i,2*j)  = (zw(2,j,i)+zw(2,j,n2+2-i))*half
3242      zt(2,i,2*j)  =-(zw(1,j,i)-zw(1,j,n2+2-i))*half
3243    end do
3244  end do
3245 
3246 end subroutine unswitchreal

m_fftcore/unswitchreal_cent [ Functions ]

[ Top ] [ m_fftcore ] [ Functions ]

NAME

 unswitchreal_cent

FUNCTION

INPUTS

  n1dfft=Number of 1D FFTs to perform
  max2=Max G_y in the small box enclosing the G-sphere.
  n2=Dimension of the transform along y
  lot=Cache blocking factor.
  n1zt
  lzt
  zw(2,lot,n2)=Cache working array

OUTPUT

  zt(2,lzt,n1)

PARENTS

SOURCE

3273 pure subroutine unswitchreal_cent(n1dfft,max2,n2,lot,n1zt,lzt,zw,zt)
3274 
3275 
3276 !This section has been created automatically by the script Abilint (TD).
3277 !Do not modify the following lines by hand.
3278 #undef ABI_FUNC
3279 #define ABI_FUNC 'unswitchreal_cent'
3280 !End of the abilint section
3281 
3282  implicit none
3283 
3284 !Arguments ------------------------------------
3285  integer,intent(in) :: n1dfft,max2,n2,lot,n1zt,lzt
3286  real(dp),intent(in) :: zw(2,lot,n2)
3287  real(dp),intent(inout) :: zt(2,lzt,n1zt)
3288 
3289 !Local variables-------------------------------
3290  integer :: i,j
3291 ! *************************************************************************
3292 
3293  do j=1,n1dfft
3294    zt(1,1,2*j-1)=zw(1,j,1)
3295    zt(2,1,2*j-1)=zero
3296    zt(1,1,2*j)  =zw(2,j,1)
3297    zt(2,1,2*j)  =zero
3298  end do
3299 
3300  do i=2,max2+1
3301    do j=1,n1dfft
3302      zt(1,i,2*j-1)= (zw(1,j,i)+zw(1,j,n2+2-i))*half
3303      zt(2,i,2*j-1)= (zw(2,j,i)-zw(2,j,n2+2-i))*half
3304      zt(1,i,2*j)  = (zw(2,j,i)+zw(2,j,n2+2-i))*half
3305      zt(2,i,2*j)  =-(zw(1,j,i)-zw(1,j,n2+2-i))*half
3306    end do
3307  end do
3308 
3309 !       Here, zero and positive frequencies
3310 !        do 90,j=1,n1dfft
3311 !        do 90,i=1,max2+1
3312 !        zt(1,i,j)=zw(1,j,i)
3313 !        zt(2,i,j)=zw(2,j,i)
3314 !90      continue
3315 
3316 !       Here, negative frequencies
3317 !        if(m2>=max2+2)then
3318 !         do 110,j=1,n1dfft
3319 !         do 110,i=max2+2,m2
3320 !         zt(1,i,j)=zw(1,j,i+n2-m2)
3321 !         zt(2,i,j)=zw(2,j,i+n2-m2)
3322 !110      continue
3323 !        end if
3324 
3325 end subroutine unswitchreal_cent