TABLE OF CONTENTS

ABINIT/m_cplxtools [ Modules ]

[ Top ] [ Modules ]

NAME

  m_cplxtools

FUNCTION

 This module defines helper functions to operate on complex arrays (mainly used in the GW code)

COPYRIGHT

 Copyright (C) 1992-2024 ABINIT group (MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

NOTES

 1) The convention about names of interfaced routine is: cplx_<name>,
    where <name> is equal to the name of the standard BLAS routine

m_cplxtools/cplx_addtorho_dpc [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_addtorho_dpc

FUNCTION

  Add |ur|**2 to the ground-states density rho.
    rho = rho + weight_r * |ur|**2

INPUTS

  nx,ny,nz=physical dimension of the FFT box.
  ldx,ldy,ldz=leading dimensions of the arrays.
  ndat=number of contributions to accumulate.
  weight_r=weight used for the accumulation of the density in real space
  ur(ldx,ldy,ldz*ndat)=wavefunctions in real space

SIDE EFFECTS

  rho(ldx,ldy,ldz) = contains the input density at input,
                     modified in input with the contribution gived by ur.

SOURCE

1046 subroutine cplx_addtorho_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,weight_r,ur,rho)
1047 
1048 !Arguments ------------------------------------
1049 !scalars
1050  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
1051  real(dp),intent(in) :: weight_r
1052 !arrays
1053  complex(dpc),intent(in) :: ur(ldx*ldy*ldz*ndat)
1054  real(dp),intent(inout) :: rho(ldx*ldy*ldz)
1055 
1056 !Local variables-------------------------------
1057 !scalars
1058  integer :: ix,iy,iz,dat,ifft,ldxyz,pad_box,padz,pady
1059 
1060 ! *************************************************************************
1061 
1062  ldxyz = ldx*ldy*ldz
1063 
1064  if (ndat==1) then
1065 !$OMP PARALLEL DO PRIVATE(padz, pady, ifft)
1066    do iz=1,nz
1067      padz = (iz-1)*ldx*ldy
1068      do iy=1,ny
1069        pady = (iy-1)*ldx
1070        do ix=1,nx
1071          ifft = ix + pady + padz
1072          rho(ifft) = rho(ifft) + weight_r * (REAL(ur(ifft))**2 + AIMAG(ur(ifft))**2)
1073        end do
1074      end do
1075    end do
1076 
1077  else
1078 ! It would be nice to use $OMP PARALLEL DO REDUCTION(+:rho)
1079 ! but it's risky as the private rho is allocated on the stack of the thread.
1080 !$OMP PARALLEL PRIVATE(pad_box, padz, pady, ifft)
1081    do dat=1,ndat
1082      pad_box = (dat-1)*ldxyz
1083 !$OMP DO
1084      do iz=1,nz
1085        padz = (iz-1)*ldx*ldy
1086        do iy=1,ny
1087          pady = (iy-1)*ldx
1088          do ix=1,nx
1089            ifft = ix + pady + padz
1090            rho(ifft) = rho(ifft) + weight_r * (REAL(ur(ifft+pad_box)**2 + AIMAG(ur(ifft+pad_box))**2))
1091          end do
1092        end do
1093      end do
1094 !$OMP END DO NOWAIT
1095    end do
1096 !$OMP END PARALLEL
1097  end if
1098 
1099 end subroutine cplx_addtorho_dpc

m_cplxtools/cplx_box2gsph_dpc [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_box2gsph_dpc

FUNCTION

   Transfer data from the FFT box to the G-sphere. Target DPC complex array.

INPUTS

  nx,ny,nz=physical dimension of the FFT box.
  ldx,ldy,ldz=Logical dimensions of the arrays.
  ndat=number of data in iarrbox
  npw_k=Number of planewaves in the G-sphere.
  kg_k(3,npw_k)=Reduced coordinates of the G-vectoes.
  iarrbox(ldx*ldy*ldz*ndat)=Complex Input arrays on the FFT box.
  [rscal] = Scaling factor

OUTPUT

  oarrsph(npw_k*ndat)=Complex Data defined on the G-sphere.

SOURCE

506 subroutine cplx_box2gsph_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,npw_k,kg_k,iarrbox,oarrsph,rscal)
507 
508 !Arguments ------------------------------------
509 !scalars
510  integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat
511  real(dp),optional,intent(in) :: rscal
512 !arrays
513  integer,intent(in) :: kg_k(3,npw_k)
514  complex(dpc),intent(in) :: iarrbox(ldx*ldy*ldz*ndat)
515  complex(dpc),intent(out) :: oarrsph(npw_k*ndat)
516 
517 !Local variables-------------------------------
518 !scalars
519  integer :: ig,ix,iy,iz,dat,pad_sph,pad_box,ifft,ldxyz
520 
521 ! *************************************************************************
522 
523  ldxyz = ldx*ldy*ldz
524  if (.not. PRESENT(rscal)) then
525    !
526    if (ndat==1) then
527 !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
528      do ig=1,npw_k
529        ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
530        iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
531        iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
532        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
533        oarrsph(ig) = iarrbox(ifft)
534      end do
535    else
536 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
537      do dat=1,ndat
538        pad_sph = (dat-1)*npw_k
539        pad_box = (dat-1)*ldxyz
540        do ig=1,npw_k
541          ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
542          iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
543          iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
544          ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
545          oarrsph(ig+pad_sph) = iarrbox(ifft+pad_box)
546        end do
547      end do
548    end if
549    !
550  else
551    if (ndat==1) then
552 !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
553      do ig=1,npw_k
554        ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
555        iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
556        iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
557        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
558        oarrsph(ig) = iarrbox(ifft) * rscal
559      end do
560    else
561 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
562      do dat=1,ndat
563        pad_sph = (dat-1)*npw_k
564        pad_box = (dat-1)*ldxyz
565        do ig=1,npw_k
566          ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
567          iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
568          iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
569          ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
570          oarrsph(ig+pad_sph) = iarrbox(ifft+pad_box) * rscal
571        end do
572      end do
573    end if
574  end if
575 
576 end subroutine cplx_box2gsph_dpc

m_cplxtools/cplx_box2gsph_spc [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_box2gsph_spc

FUNCTION

   Transfer data from the FFT box to the G-sphere. Target SPC complex array.

INPUTS

  nx,ny,nz=physical dimension of the FFT box.
  ldx,ldy,ldz=Logical dimensions of the arrays.
  ndat=number of data in iarrbox
  npw_k=Number of planewaves in the G-sphere.
  kg_k(3,npw_k)=Reduced coordinates of the G-vectoes.
  iarrbox(ldx*ldy*ldz*ndat)=Complex Input arrays on the FFT box.
  [rscal] = Scaling factor

OUTPUT

  oarrsph(npw_k*ndat)=Complex Data defined on the G-sphere.

SOURCE

408 subroutine cplx_box2gsph_spc(nx,ny,nz,ldx,ldy,ldz,ndat,npw_k,kg_k,iarrbox,oarrsph,rscal)
409 
410 !Arguments ------------------------------------
411 !scalars
412  integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat
413  real(sp),optional,intent(in) :: rscal
414 !arrays
415  integer,intent(in) :: kg_k(3,npw_k)
416  complex(spc),intent(in) :: iarrbox(ldx*ldy*ldz*ndat)
417  complex(spc),intent(out) :: oarrsph(npw_k*ndat)
418 
419 !Local variables-------------------------------
420 !scalars
421  integer :: ig,ix,iy,iz,dat,pad_sph,pad_box,ifft,ldxyz
422 
423 ! *************************************************************************
424 
425  ldxyz = ldx*ldy*ldz
426  if (.not. PRESENT(rscal)) then
427    !
428    if (ndat==1) then
429 !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
430      do ig=1,npw_k
431        ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
432        iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
433        iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
434        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
435        oarrsph(ig) = iarrbox(ifft)
436      end do
437    else
438 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
439      do dat=1,ndat
440        pad_sph = (dat-1)*npw_k
441        pad_box = (dat-1)*ldxyz
442        do ig=1,npw_k
443          ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
444          iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
445          iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
446          ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
447          oarrsph(ig+pad_sph) = iarrbox(ifft+pad_box)
448        end do
449      end do
450    end if
451    !
452  else
453    if (ndat==1) then
454 !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
455      do ig=1,npw_k
456        ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
457        iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
458        iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
459        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
460        oarrsph(ig) = iarrbox(ifft) * rscal
461      end do
462    else
463 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
464      do dat=1,ndat
465        pad_sph = (dat-1)*npw_k
466        pad_box = (dat-1)*ldxyz
467        do ig=1,npw_k
468          ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
469          iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
470          iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
471          ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
472          oarrsph(ig+pad_sph) = iarrbox(ifft+pad_box) * rscal
473        end do
474      end do
475    end if
476  end if
477 
478 end subroutine cplx_box2gsph_spc

m_cplxtools/cplx_filter [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_filter

FUNCTION

  Set all the elements of x to zero where mask is .TRUE.

INPUTS

  n=Specifies the number of elements in vectors x and y.
  mask(n)=Logical array.

SIDE EFFECTS

  x(n)=See description.

SOURCE

153 subroutine cplx_filter(n, x, mask)
154 
155 !Arguments ------------------------------------
156 !scalars
157  integer,intent(in) :: n
158 !arrays
159  complex(dpc),intent(inout) :: x(n)
160  logical,intent(in) :: mask(n)
161 
162 ! *************************************************************************
163 
164  where (mask)
165    x = czero
166  end where
167 
168 end subroutine cplx_filter

m_cplxtools/cplx_fromreal [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_fromreal

FUNCTION

  Convert a real array with (real,imag) part to a complex array

INPUTS

  n = Specifies the number of elements in ocplx
  ireal(2*n)=Input real array.

OUTPUT

  ocplx(n)=Output complex array

SOURCE

112 subroutine cplx_fromreal(n,ireal,ocplx)
113 
114 !Arguments ------------------------------------
115 !scalars
116  integer,intent(in) :: n
117 !arrays
118  real(dp),intent(in) :: ireal(2,n)
119  complex(dpc),intent(out) :: ocplx(n)
120 
121 !Local variables ------------------------------
122 !scalars
123  integer :: ii
124 
125 ! *************************************************************************
126 
127 !$OMP PARALLEL DO PRIVATE(ii)
128  do ii=1,n
129    ocplx(ii) = DCMPLX(ireal(1,ii),ireal(2,ii))
130  end do
131 
132 end subroutine cplx_fromreal

m_cplxtools/cplx_gsph2box_dpc [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

 cplx_gsph2box_dpc

FUNCTION

 Array iarrsph is defined in sphere with npw points. Insert iarrsph inside box
 of nx*ny*nz points to define array oarrbox for fft box. rest of oarrbox is filled with 0 s.
 targer: DPC complex arrays

INPUTS

 iarrsph(2,npw*ndat)= contains values for npw G vectors in basis sphere
 ndat=number of FFT to perform.
 npw=number of G vectors in basis at this k point
 oarrbox(2,ldx*ldy*ldz*ndat) = fft box
 nx,ny,nz=physical dimension of the box (oarrbox)
 ldx,ldy,ldz=memory dimension of oarrbox
 kg_k(3,npw)=integer coordinates of G vectors in basis sphere
 istwf_k=option parameter that describes the storage of wfs

OUTPUT

   oarrbox(ldx*ldy*ldz*ndat)

NOTES

 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.

SOURCE

774 subroutine cplx_gsph2box_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,npw,istwf_k,kg_k,iarrsph,oarrbox)
775 
776 !Arguments ------------------------------------
777 !scalars
778  integer,intent(in) :: istwf_k,nx,ny,nz,ldx,ldy,ldz,ndat,npw
779 !arrays
780  integer,intent(in) :: kg_k(3,npw)
781  complex(dpc),intent(in) :: iarrsph(npw*ndat)
782  complex(dpc),intent(out) :: oarrbox(ldx*ldy*ldz*ndat)
783 
784 !Local variables-------------------------------
785 !scalars
786  integer,parameter :: me_g0=1
787  integer :: ix,ixinv,iy,iyinv,iz,izinv,dat,ipw,npwmin,pad_box,pad_sph,ifft,ifft_inv,ldxyz
788  !character(len=500) :: msg
789 !arrays
790  integer,allocatable :: ixinver(:),iyinver(:),izinver(:)
791 
792 ! *************************************************************************
793 
794 !In the case of special k-points, invariant under time-reversal,
795 !but not Gamma, initialize the inverse coordinates
796 !Remember indeed that
797 !u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^*
798 !and therefore:
799 !u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
800  if (istwf_k>=2) then
801    ABI_MALLOC(ixinver,(nx))
802    ABI_MALLOC(iyinver,(ny))
803    ABI_MALLOC(izinver,(nz))
804    if ( ANY(istwf_k==(/2,4,6,8/)) ) then
805      ixinver(1)=1
806      do ix=2,nx
807        ixinver(ix)=nx+2-ix
808      end do
809    else
810      do ix=1,nx
811        ixinver(ix)=nx+1-ix
812      end do
813    end if
814    if (istwf_k>=2 .and. istwf_k<=5) then
815      iyinver(1)=1
816      do iy=2,ny
817        iyinver(iy)=ny+2-iy
818      end do
819    else
820      do iy=1,ny
821        iyinver(iy)=ny+1-iy
822      end do
823    end if
824    if ( ANY(istwf_k==(/2,3,6,7/)) ) then
825      izinver(1)=1
826      do iz=2,nz
827        izinver(iz)=nz+2-iz
828      end do
829    else
830      do iz=1,nz
831        izinver(iz)=nz+1-iz
832      end do
833    end if
834  end if
835 
836  ldxyz = ldx*ldy*ldz
837 
838  if (istwf_k==1) then
839 
840 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
841    do dat=1,ndat
842      pad_sph = (dat-1)*npw
843      pad_box = (dat-1)*ldxyz
844      oarrbox(1+pad_box:ldxyz+pad_box) = czero_dpc ! zero the sub-array
845      do ipw=1,npw
846        ix=kg_k(1,ipw); if (ix<0) ix=ix+nx; ix=ix+1
847        iy=kg_k(2,ipw); if (iy<0) iy=iy+ny; iy=iy+1
848        iz=kg_k(3,ipw); if (iz<0) iz=iz+nz; iz=iz+1
849        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
850 #if defined __INTEL_COMPILER && defined HAVE_OPENMP
851        if (ifft==0) then
852          ABI_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
853        end if
854 #endif
855        oarrbox(ifft+pad_box) = iarrsph(ipw+pad_sph)
856      end do
857    end do
858 
859  else if (istwf_k>=2) then
860    !
861    npwmin=1
862    if(istwf_k==2 .and. me_g0==1) then ! If gamma point, then oarrbox must be completed
863      do dat=1,ndat
864        pad_sph = (dat-1)*npw
865        pad_box = (dat-1)*ldxyz
866        oarrbox(1+pad_box) = REAL(iarrsph(1+pad_sph))
867      end do
868      npwmin=2
869    end if
870 
871 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,ixinv,iy,iyinv,iz,izinv,ifft)
872    do dat=1,ndat
873      pad_sph = (dat-1)*npw
874      pad_box = (dat-1)*ldxyz
875      oarrbox(npwmin+pad_box:ldxyz+pad_box) = czero_dpc
876      do ipw=npwmin,npw
877        ix=kg_k(1,ipw); if(ix<0)ix=ix+nx; ix=ix+1
878        iy=kg_k(2,ipw); if(iy<0)iy=iy+ny; iy=iy+1
879        iz=kg_k(3,ipw); if(iz<0)iz=iz+nz; iz=iz+1
880        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
881        ! Construct the coordinates of -k-G
882        ixinv=ixinver(ix); iyinv=iyinver(iy); izinv=izinver(iz)
883        ifft_inv = ixinv + (iyinv-1)*ldx + (izinv-1)*ldx*ldy
884 #if defined __INTEL_COMPILER && defined HAVE_OPENMP
885        if (ifft==0 .or. ifft_inv==0) then
886          ABI_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
887        end if
888 #endif
889        oarrbox(ifft    +pad_box) =        iarrsph(ipw+pad_sph)
890        oarrbox(ifft_inv+pad_box) = DCONJG(iarrsph(ipw+pad_sph))
891      end do
892    end do
893    !
894  else
895    ABI_ERROR("Wrong istwfk")
896  end if
897 
898  if (istwf_k>=2) then
899    ABI_FREE(ixinver)
900    ABI_FREE(iyinver)
901    ABI_FREE(izinver)
902  end if
903 
904 end subroutine cplx_gsph2box_dpc

m_cplxtools/cplx_gsph2box_spc [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

 cplx_gsph2box_spc

FUNCTION

 Array iarrsph is defined in sphere with npw points. Insert iarrsph inside box
 of nx*ny*nz points to define array oarrbox for fft box. rest of oarrbox is filled with 0 s.
 targer: SPC complex arrays

INPUTS

 iarrsph(2,npw*ndat)= contains values for npw G vectors in basis sphere
 ndat=number of FFT to perform.
 npw=number of G vectors in basis at this k point
 oarrbox(2,ldx*ldy*ldz*ndat) = fft box
 nx,ny,nz=physical dimension of the box (oarrbox)
 ldx,ldy,ldz=memory dimension of oarrbox
 kg_k(3,npw)=integer coordinates of G vectors in basis sphere
 istwf_k=option parameter that describes the storage of wfs

OUTPUT

   oarrbox(ldx*ldy*ldz*ndat)

NOTES

 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.

SOURCE

610 subroutine cplx_gsph2box_spc(nx,ny,nz,ldx,ldy,ldz,ndat,npw,istwf_k,kg_k,iarrsph,oarrbox)
611 
612 !Arguments ------------------------------------
613 !scalars
614  integer,intent(in) :: istwf_k,nx,ny,nz,ldx,ldy,ldz,ndat,npw
615 !arrays
616  integer,intent(in) :: kg_k(3,npw)
617  complex(spc),intent(in) :: iarrsph(npw*ndat)
618  complex(spc),intent(out) :: oarrbox(ldx*ldy*ldz*ndat)
619 
620 !Local variables-------------------------------
621 !scalars
622  integer,parameter :: me_g0=1
623  integer :: ix,ixinv,iy,iyinv,iz,izinv,dat,ipw,npwmin,pad_box,pad_sph,ifft,ifft_inv,ldxyz
624  !character(len=500) :: msg
625 !arrays
626  integer,allocatable :: ixinver(:),iyinver(:),izinver(:)
627 
628 ! *************************************************************************
629 
630 !In the case of special k-points, invariant under time-reversal,
631 !but not Gamma, initialize the inverse coordinates
632 !Remember indeed that
633 !u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^*
634 !and therefore:
635 !u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
636  if (istwf_k>=2) then
637    ABI_MALLOC(ixinver,(nx))
638    ABI_MALLOC(iyinver,(ny))
639    ABI_MALLOC(izinver,(nz))
640    if ( ANY(istwf_k==(/2,4,6,8/)) ) then
641      ixinver(1)=1
642      do ix=2,nx
643        ixinver(ix)=nx+2-ix
644      end do
645    else
646      do ix=1,nx
647        ixinver(ix)=nx+1-ix
648      end do
649    end if
650    if (istwf_k>=2 .and. istwf_k<=5) then
651      iyinver(1)=1
652      do iy=2,ny
653        iyinver(iy)=ny+2-iy
654      end do
655    else
656      do iy=1,ny
657        iyinver(iy)=ny+1-iy
658      end do
659    end if
660    if ( ANY(istwf_k==(/2,3,6,7/)) ) then
661      izinver(1)=1
662      do iz=2,nz
663        izinver(iz)=nz+2-iz
664      end do
665    else
666      do iz=1,nz
667        izinver(iz)=nz+1-iz
668      end do
669    end if
670  end if
671 
672  ldxyz = ldx*ldy*ldz
673 
674  if (istwf_k==1) then
675 
676 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
677    do dat=1,ndat
678      pad_sph = (dat-1)*npw
679      pad_box = (dat-1)*ldxyz
680      oarrbox(1+pad_box:ldxyz+pad_box) = czero_spc ! zero the sub-array
681      do ipw=1,npw
682        ix=kg_k(1,ipw); if (ix<0) ix=ix+nx; ix=ix+1
683        iy=kg_k(2,ipw); if (iy<0) iy=iy+ny; iy=iy+1
684        iz=kg_k(3,ipw); if (iz<0) iz=iz+nz; iz=iz+1
685        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
686 #if defined __INTEL_COMPILER && defined HAVE_OPENMP
687        if (ifft==0) then
688          ABI_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
689        end if
690 #endif
691        oarrbox(ifft+pad_box) = iarrsph(ipw+pad_sph)
692      end do
693    end do
694 
695  else if (istwf_k>=2) then
696    !
697    npwmin=1
698    if(istwf_k==2 .and. me_g0==1) then ! If gamma point, then oarrbox must be completed
699      do dat=1,ndat
700        pad_sph = (dat-1)*npw
701        pad_box = (dat-1)*ldxyz
702        oarrbox(1+pad_box) = REAL(iarrsph(1+pad_sph))
703      end do
704      npwmin=2
705    end if
706 
707 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,ixinv,iy,iyinv,iz,izinv,ifft)
708    do dat=1,ndat
709      pad_sph = (dat-1)*npw
710      pad_box = (dat-1)*ldxyz
711      oarrbox(npwmin+pad_box:ldxyz+pad_box) = czero_spc
712      do ipw=npwmin,npw
713        ix=kg_k(1,ipw); if(ix<0)ix=ix+nx; ix=ix+1
714        iy=kg_k(2,ipw); if(iy<0)iy=iy+ny; iy=iy+1
715        iz=kg_k(3,ipw); if(iz<0)iz=iz+nz; iz=iz+1
716        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
717        ! Construct the coordinates of -k-G
718        ixinv=ixinver(ix); iyinv=iyinver(iy); izinv=izinver(iz)
719        ifft_inv = ixinv + (iyinv-1)*ldx + (izinv-1)*ldx*ldy
720 #if defined __INTEL_COMPILER && defined HAVE_OPENMP
721        if (ifft==0 .or. ifft_inv==0) then
722          ABI_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
723        end if
724 #endif
725        oarrbox(ifft    +pad_box) =       iarrsph(ipw+pad_sph)
726        oarrbox(ifft_inv+pad_box) = CONJG(iarrsph(ipw+pad_sph))
727      end do
728    end do
729    !
730  else
731    ABI_ERROR("Wrong istwfk")
732  end if
733 
734  if (istwf_k>=2) then
735    ABI_FREE(ixinver)
736    ABI_FREE(iyinver)
737    ABI_FREE(izinver)
738  end if
739 
740 end subroutine cplx_gsph2box_spc

m_cplxtools/cplx_real_zdotc [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_real_zdotc

FUNCTION

   Perform a vector-vector operation defined as res = REAL (\Sigma (conjg(x)*y)) where x and y are n-element vectors.

INPUTS

  n = Specifies the number of elements in vector x and y
  x,y = Input arrays.

OUTPUT

  res=Real part of the scalar product.

SOURCE

189 function cplx_real_zdotc(n,x,y) result(res)
190 
191 !Arguments ------------------------------------
192 !scalars
193  integer,intent(in) :: n
194 !arrays
195  complex(dpc),intent(in) :: x(n)
196  complex(dpc),intent(in) :: y(n)
197  real(dp) :: res
198 
199 !Local variables-------------------------------
200  real(dp),external :: ddot
201 
202 ! *************************************************************************
203 
204  res = ddot(2*n,x,1,y,1)
205 
206 end function cplx_real_zdotc

m_cplxtools/cplx_setaug_zero_dpc [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_setaug_zero_dpc

FUNCTION

  Set to zero all elements of the array that are not in the FFT box.

INPUTS

 nx,ny,nz=physical dimensions of the FFT box
 ldx,ldy,ldx=memory dimension of arr
 ndat=number of FFTs

 SIDE EFFECT
  arr(ldx,ldy,ldz*ndat)= all entries in the augmented region are set to zero

SOURCE

 984 subroutine cplx_setaug_zero_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,arr)
 985 
 986 !Arguments ------------------------------------
 987 !scalars
 988  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
 989 !arrays
 990  complex(dpc),intent(inout) :: arr(ldx,ldy,ldz*ndat)
 991 
 992 !Local variables-------------------------------
 993  integer :: iy,iz,dat,padat
 994 
 995 ! *************************************************************************
 996 
 997  if (nx /= ldx) then
 998    do iz=1,ldz*ndat
 999      do iy=1,ldy
1000        arr(nx+1:ldx,iy,iz) = czero_dpc
1001      end do
1002    end do
1003  end if
1004 
1005  if (ny /= ldy) then
1006    do iz=1,ldz*ndat
1007      arr(:,ny+1:ldy,iz) = czero_dpc
1008    end do
1009  end if
1010 
1011  if (nz /= ldz) then
1012    do dat=1,ndat
1013      padat = ldz*(dat-1)
1014      do iz=nz+1,ldz
1015        arr(:,:,iz+padat) = czero_dpc
1016      end do
1017    end do
1018  end if
1019 
1020 end subroutine cplx_setaug_zero_dpc

m_cplxtools/cplx_setaug_zero_spc [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_setaug_zero_spc

FUNCTION

  Set to zero all elements of the array that are not in the FFT box.

INPUTS

 nx,ny,nz=physical dimensions of the FFT box
 ldx,ldy,ldx=memory dimension of arr
 ndat=number of FFTs

 SIDE EFFECT
  arr(ldx,ldy,ldz*ndat)= all entries in the augmented region are set to zero

SOURCE

926 subroutine cplx_setaug_zero_spc(nx,ny,nz,ldx,ldy,ldz,ndat,arr)
927 
928 !Arguments ------------------------------------
929 !scalars
930  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
931 !arrays
932  complex(spc),intent(inout) :: arr(ldx,ldy,ldz*ndat)
933 
934 !Local variables-------------------------------
935  integer :: iy,iz,dat,padat
936 
937 ! *************************************************************************
938 
939  if (nx /= ldx) then
940    do iz=1,ldz*ndat
941      do iy=1,ldy
942        arr(nx+1:ldx,iy,iz) = czero_spc
943      end do
944    end do
945  end if
946 
947  if (ny /= ldy) then
948    do iz=1,ldz*ndat
949      arr(:,ny+1:ldy,iz) = czero_spc
950    end do
951  end if
952 
953  if (nz /= ldz) then
954    do dat=1,ndat
955      padat = ldz*(dat-1)
956      do iz=nz+1,ldz
957        arr(:,:,iz+padat) = czero_spc
958      end do
959    end do
960  end if
961 
962 end subroutine cplx_setaug_zero_spc

m_cplxtools/cplx_zaxpby [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_zaxpby

FUNCTION

  Scales two vectors, adds them to one another and stores result in the vector.
  y := a*x + b*y

INPUTS

 n = the number of elements in vectors x and y.
 a = Specifies the scalar a.
 x = Array.
 b = Specifies the scalar b.
 y = Array

OUTPUT

 y Contains the updated vector y.

SOURCE

231 subroutine cplx_zaxpby(n,a,x,b,y)
232 
233 !Arguments ------------------------------------
234 !scalars
235  integer,intent(in) :: n
236  complex(dpc),intent(in) :: a,b
237 !arrays
238  complex(dpc),intent(in) :: x(n)
239  complex(dpc),intent(inout) :: y(n)
240 
241 ! *************************************************************************
242 
243 #ifdef HAVE_LINALG_AXPBY
244  call zaxpby(n, a, x, 1, b, y, 1)
245 #else
246  call zscal(n, b, y, 1)
247  call zaxpy(n, a, x, 1, y,1)
248 #endif
249 
250 end subroutine cplx_zaxpby

m_cplxtools/cplx_zgemm [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_zgemm

FUNCTION

  The ?gemm routines perform a matrix-matrix operation with general matrices.
  The operation is defined as C := alpha*op(A)*op(B) + beta*C,
  where:

  op(x) is one of op(x) = x, or op(x) = x', or op(x) = conjg(x'),

  alpha and beta are scalars,
  A, B and C are matrices:
  op(A) is an m-by-k matrix,
  op(B) is a k-by-n matrix,
  C is an m-by-n matrix.

INPUTS

OUTPUT

SOURCE

343 subroutine cplx_zgemm(transa,transb,npws,ncola,ncolb,amat,bmat,cmat,alpha,beta)
344 
345 !Arguments ------------------------------------
346 !scalars
347  integer,intent(in) :: npws,ncola,ncolb
348  complex(dpc),optional,intent(in) :: alpha,beta
349  character(len=1),intent(in) :: transa,transb
350 !arrays
351  complex(dpc),intent(in) :: amat(npws*ncola)
352  complex(dpc),intent(in) :: bmat(npws*ncolb)
353  complex(dpc),intent(inout) :: cmat(*)
354 
355 !Local variables-------------------------------
356 !scalars
357  integer :: mm,nn,kk,lda,ldb,ldc
358  complex(dpc) :: my_alpha,my_beta
359 
360 ! *************************************************************************
361 
362  lda = npws
363  ldb = npws
364 
365  mm  = npws
366  nn  = ncolb
367  kk  = ncola
368 
369  if (toupper(transa) /= 'N') then
370    mm = ncola
371    kk = npws
372  end if
373  if (toupper(transb) /= 'N') nn = npws
374 
375  ldc = mm
376 
377  my_alpha = cone_dpc;  if (PRESENT(alpha)) my_alpha = alpha
378  my_beta  = czero_dpc; if (PRESENT(beta))  my_beta  = beta
379 
380  call ZGEMM(transa,transb,mm,nn,kk,my_alpha,amat,lda,bmat,ldb,my_beta,cmat,ldc)
381 
382 end subroutine cplx_zgemm

m_cplxtools/cplx_zgemv [ Functions ]

[ Top ] [ m_cplxtools ] [ Functions ]

NAME

  cplx_zgemv

FUNCTION

 The ?gemv routines perform a matrix-vector operation defined as

 y := alpha*A*x + beta*y,
 or
 y := alpha*A'*x + beta*y,
 or
 y := alpha*conjg(A')*x + beta*y,

 where: alpha and beta are scalars, x and y are vectors, A is an m-by-n matrix.

INPUTS

OUTPUT

SOURCE

276 subroutine cplx_zgemv(trans,nrows,ncols,mat,vec,matvec,alpha,beta)
277 
278 !Arguments ------------------------------------
279 !scalars
280  integer,intent(in) :: nrows,ncols
281  complex(dpc),optional,intent(in) :: alpha,beta
282  character(len=1),intent(in) :: trans
283 !arrays
284  complex(dpc),intent(in) :: mat(nrows*ncols)
285  complex(dpc),intent(in) :: vec(*)
286  complex(dpc),intent(inout) :: matvec(*)
287 
288 !Local variables-------------------------------
289 !scalars
290  integer :: mm,nn,kk,lda,ldb,ldc
291  complex(dpc) :: my_alpha,my_beta
292 
293 ! *************************************************************************
294 
295  lda = nrows
296  mm  = nrows
297  nn  = 1
298  kk  = ncols
299 
300  if (toupper(trans) /= 'N') then
301    mm = ncols
302    kk = nrows
303  end if
304 
305  ldb = kk
306  ldc = mm
307 
308  my_alpha = cone_dpc;  if (PRESENT(alpha)) my_alpha = alpha
309  my_beta  = czero_dpc; if (PRESENT(beta))  my_beta  = beta
310 
311  call ZGEMM(trans,"N",mm,nn,kk,my_alpha,mat,lda,vec,ldb,my_beta,matvec,ldc)
312 
313  !call ZGEMV(trans,mm,nn,my_alpha,mat,lda,vec,1,my_beta,matvec,1)
314 
315 end subroutine cplx_zgemv