LCOV - code coverage report
Current view: top level - eigen - hsmt_sph.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 73 77 94.8 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_hsmt_sph
       8             :    USE m_juDFT
       9             :    IMPLICIT NONE
      10             : 
      11             :    INTERFACE hsmt_sph
      12             :       module procedure hsmt_sph_cpu
      13             : #ifdef CPP_GPU
      14             :       module procedure hsmt_sph_gpu
      15             : #endif
      16             :    END INTERFACE
      17             : 
      18             : CONTAINS
      19             : 
      20        4396 : SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_shift,usdus,fj,gj,smat,hmat)
      21             :    USE m_constants, ONLY : fpi_const,tpi_const
      22             :    USE m_types
      23             :    IMPLICIT NONE
      24             :    TYPE(t_input),INTENT(IN)      :: input
      25             :    TYPE(t_mpi),INTENT(IN)        :: mpi
      26             :    TYPE(t_noco),INTENT(IN)       :: noco
      27             :    TYPE(t_atoms),INTENT(IN)      :: atoms
      28             :    TYPE(t_lapw),INTENT(IN)       :: lapw
      29             :    TYPE(t_usdus),INTENT(IN)      :: usdus
      30             :    CLASS(t_mat),INTENT(INOUT)    :: smat,hmat
      31             :    !     ..
      32             :    !     .. Scalar Arguments ..
      33             :    INTEGER, INTENT (IN) :: n,isp,iintsp,jintsp
      34             :    COMPLEX, INTENT(IN)  :: chi
      35             :    !     ..
      36             :    !     .. Array Arguments ..
      37             :    REAL,    INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,input%jspins)
      38             :    REAL,    INTENT (IN) :: e_shift!(atoms%ntype,input%jspins)
      39             :    REAL,    INTENT (IN) :: fj(:,0:,:),gj(:,0:,:)
      40             :    !     ..
      41             :    !     .. Local Scalars ..
      42             :    REAL tnn(3), elall,fjkiln,gjkiln,ddnln,ski(3)
      43             :    REAL apw_lo1,apw_lo2,w1
      44             : 
      45             :    INTEGER kii,ki,kj,l,nn,kj_end,l3,jv,kj_off,kj_vec
      46             : 
      47             :    !     ..
      48             :    !     .. Local Arrays ..
      49        8792 :    REAL fleg1(0:atoms%lmaxd),fleg2(0:atoms%lmaxd),fl2p1(0:atoms%lmaxd)
      50             :    REAL qssbti(3),qssbtj(3)
      51        4396 :    REAL, ALLOCATABLE :: plegend(:,:)
      52        4396 :    REAL, ALLOCATABLE :: xlegend(:)
      53        4396 :    REAL, ALLOCATABLE :: VecHelpS(:),VecHelpH(:)
      54        4396 :    REAL, ALLOCATABLE :: cph_re(:), cph_im(:)
      55        4396 :    REAL, ALLOCATABLE :: dot(:), fct(:), fct2(:)
      56             :    INTEGER, PARAMETER :: NVEC = 128 
      57             :    INTEGER :: NVEC_rem  !remainder
      58             : 
      59        4396 :    CALL timestart("spherical setup")
      60             : 
      61       45232 :    DO l = 0,atoms%lmaxd
      62       40836 :       fleg1(l) = REAL(l+l+1)/REAL(l+1)
      63       40836 :       fleg2(l) = REAL(l)/REAL(l+1)
      64       45232 :       fl2p1(l) = REAL(l+l+1)/fpi_const
      65             :    END DO ! l
      66             : !$OMP     PARALLEL DEFAULT(NONE)&
      67             : !$OMP     SHARED(lapw,atoms,noco,mpi,input,usdus,smat,hmat)&
      68             : !$OMP     SHARED(jintsp,iintsp,n,fleg1,fleg2,fj,gj,isp,fl2p1,el,e_shift,chi)&
      69             : !$OMP     PRIVATE(kii,ki,ski,kj,kj_off,kj_vec,plegend,xlegend,l,l3,kj_end,qssbti,qssbtj,fct2)&
      70             : !$OMP     PRIVATE(cph_re,cph_im,dot,nn,tnn,fjkiln,gjkiln)&
      71             : !$OMP     PRIVATE(w1,apw_lo1,apw_lo2,ddnln,elall,fct)&
      72       13188 : !$OMP     PRIVATE(VecHelpS,VecHelpH,NVEC_rem)
      73        4396 :    ALLOCATE(cph_re(NVEC),cph_im(NVEC))
      74        4396 :    ALLOCATE(dot(NVEC),fct(NVEC),fct2(NVEC))
      75        4396 :    ALLOCATE(plegend(NVEC,0:2))
      76        4396 :    ALLOCATE(xlegend(NVEC))
      77        4396 :    ALLOCATE(VecHelpS(NVEC),VecHelpH(NVEC))
      78       17584 :    qssbti=MERGE(- noco%qss/2,+ noco%qss/2,jintsp.EQ.1)
      79       17584 :    qssbtj=MERGE(- noco%qss/2,+ noco%qss/2,iintsp.EQ.1)
      80        4396 : !$OMP      DO SCHEDULE(DYNAMIC,1)
      81             :    DO  ki =  mpi%n_rank+1, lapw%nv(jintsp), mpi%n_size
      82      516944 :       kii=(ki-1)/mpi%n_size+1
      83     2067776 :       ski(1:3) = lapw%gvec(1:3,ki,jintsp) + qssbti(1:3)
      84             : 
      85      516944 :       IF (smat%l_real) THEN
      86             :          kj_end = ki
      87             :       ELSE
      88      209256 :          kj_end = MIN(ki,lapw%nv(iintsp))
      89             :       ENDIF
      90      516944 :       NVEC_rem = NVEC
      91     1588178 :       DO  kj_off = 1, kj_end, NVEC
      92     1071234 :          kj_vec = kj_off - 1 + NVEC
      93     1071234 :          IF (kj_vec > kj_end) THEN
      94      514662 :              kj_vec = kj_end
      95      514662 :              NVEC_rem = kj_end - kj_off + 1
      96             :          ENDIF
      97             : 
      98             :          !--->          update overlap and l-diagonal hamiltonian matrix
      99   275307138 :          VecHelpS = 0.0
     100   275307138 :          VecHelpH = 0.0
     101             :          
     102             :          !--->       x for legendre polynomials
     103   202634210 :          DO jv = 0, NVEC_rem-1
     104   100781488 :             kj = jv + kj_off
     105   404197186 :             xlegend(jv+1) = DOT_PRODUCT(lapw%gk(1:3,kj,iintsp), lapw%gk(1:3,ki,jintsp))
     106             :          END DO ! kj
     107             :          
     108    10837908 :          DO  l = 0,atoms%lmax(n)
     109             :          
     110     9766674 :             fjkiln = fj(ki,l,jintsp)
     111     9766674 :             gjkiln = gj(ki,l,jintsp)
     112             : 
     113     9766674 :             IF (input%l_useapw) THEN
     114           0 :                w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + usdus%us(l,n,isp)*usdus%duds(l,n,isp) )
     115           0 :                apw_lo1 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( gjkiln * w1 + fjkiln * usdus%us(l,n,isp)  * usdus%dus(l,n,isp) )
     116           0 :                apw_lo2 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( fjkiln * w1 + gjkiln * usdus%uds(l,n,isp) * usdus%duds(l,n,isp) )
     117             :             ENDIF ! useapw
     118             : 
     119     9766674 :             ddnln = usdus%ddn(l,n,isp)
     120     9766674 :             elall = el(l,n,isp)
     121    17313756 :             IF (l<=atoms%lnonsph(n)) elall=elall-e_shift!(isp)
     122             : 
     123             :             !--->       legendre polynomials
     124     9766674 :             l3 = modulo(l, 3)
     125     9766674 :             IF (l == 0) THEN
     126   138189186 :                plegend(:,0) = 1.0
     127     8695440 :             ELSE IF (l == 1) THEN
     128     1071234 :                plegend(:NVEC_REM,1) = xlegend(:NVEC_REM)
     129             :             ELSE
     130     7624206 :                plegend(:NVEC_REM,l3) = fleg1(l-1)*xlegend(:NVEC_REM)*plegend(:NVEC_REM,modulo(l-1,3)) - fleg2(l-1)*plegend(:NVEC_REM,modulo(l-2,3))
     131             :             END IF ! l
     132             : 
     133   924056002 :             fct(:NVEC_REM)  = plegend(:NVEC_REM,l3)*fl2p1(l)       * ( fjkiln*fj(kj_off:kj_vec,l,iintsp) + gjkiln*gj(kj_off:kj_vec,l,iintsp)*ddnln )
     134  1838345330 :             fct2(:NVEC_REM) = plegend(:NVEC_REM,l3)*fl2p1(l) * 0.5 * ( gjkiln*fj(kj_off:kj_vec,l,iintsp) + fjkiln*gj(kj_off:kj_vec,l,iintsp) )
     135             : 
     136  1838345330 :             VecHelpS(:NVEC_REM) = VecHelpS(:NVEC_REM) + fct(:NVEC_REM)
     137  1838345330 :             VecHelpH(:NVEC_REM) = VecHelpH(:NVEC_REM) + fct(:NVEC_REM)*elall + fct2(:NVEC_REM)
     138             : 
     139    10837908 :             IF (input%l_useapw) THEN
     140           0 :                VecHelpH(:NVEC_REM) = VecHelpH(:NVEC_REM) + plegend(:NVEC_REM,l3) * ( apw_lo1*fj(kj_off:kj_vec,l,iintsp) + apw_lo2*gj(kj_off:kj_vec,l,iintsp) )
     141             :             ENDIF ! useapw
     142             : 
     143             :             !--->          end loop over l
     144             :          ENDDO ! l
     145             : 
     146             :          !--->             set up phase factors
     147   275307138 :          cph_re = 0.0
     148   275307138 :          cph_im = 0.0
     149     5519166 :          DO nn = SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
     150     6838344 :             tnn(1:3) = tpi_const*atoms%taual(1:3,nn)
     151   336355928 :             DO jv = 0, NVEC_rem-1
     152   167323171 :                kj = jv + kj_off
     153   671002270 :                dot(jv+1) = DOT_PRODUCT(ski(1:3) - lapw%gvec(1:3,kj,iintsp) - qssbtj(1:3), tnn(1:3))
     154             :             END DO ! kj
     155   169032757 :             cph_re(:NVEC_REM) = cph_re(:NVEC_REM) + COS(dot(:NVEC_REM))
     156   337427162 :             cph_im(:NVEC_REM) = cph_im(:NVEC_REM) - SIN(dot(:NVEC_REM))
     157             :             ! IF (iintsp.NE.jintsp) cph_im=-cph_im
     158             :          END DO ! nn
     159             :          
     160     1588178 :          IF (smat%l_real) THEN
     161             :             smat%data_r(kj_off:kj_vec,kii) = &
     162    83198264 :             smat%data_r(kj_off:kj_vec,kii) + cph_re(:NVEC_REM) * VecHelpS(:NVEC_REM)
     163             :             hmat%data_r(kj_off:kj_vec,kii) = &
     164      805974 :             hmat%data_r(kj_off:kj_vec,kii) + cph_re(:NVEC_REM) * VecHelpH(:NVEC_REM)
     165             :          ELSE  ! real
     166             :             smat%data_c(kj_off:kj_vec,kii) = &
     167    18654458 :             smat%data_c(kj_off:kj_vec,kii) + chi*cmplx(cph_re(:NVEC_REM),cph_im(:NVEC_REM)) * VecHelpS(:NVEC_REM)
     168             :             hmat%data_c(kj_off:kj_vec,kii) = &
     169      265260 :             hmat%data_c(kj_off:kj_vec,kii) + chi*cmplx(cph_re(:NVEC_REM),cph_im(:NVEC_REM)) * VecHelpH(:NVEC_REM)
     170             :          ENDIF ! real
     171             : 
     172             :       END DO ! kj_off
     173             : 
     174             :       !--->    end loop over ki
     175             :    ENDDO
     176             : !$OMP     END DO
     177        4396 :    DEALLOCATE(plegend)
     178        4396 :    DEALLOCATE(VecHelpS,VecHelpH)
     179             : !$OMP     END PARALLEL
     180        4396 :    CALL timestop("spherical setup")
     181             : 
     182        4396 :    RETURN
     183        8792 : END SUBROUTINE hsmt_sph_cpu
     184             : 
     185             : 
     186             : #ifdef CPP_GPU
     187             :    ATTRIBUTES(global)&
     188             :       SUBROUTINE HsmtSphGpuKernel_real(loop_size,iintsp,jintsp,nv,lmaxd,lmax,ki_start,ki_end,ki_step,nn_start,nn_end,&
     189             :                                        lnonsph,qssbti,qssbtj,gvec,gk,fleg1,fleg2,fl2p1,fl2p1bt,fj,gj,taual,ddn,el,e_shift,&
     190             :                                        smat_data,hmat_data,&
     191             :                                        uds,dus,us,duds,rmt)
     192             :    INTEGER, VALUE, INTENT(IN) :: loop_size
     193             :    INTEGER,VALUE,INTENT(IN) :: iintsp,jintsp,lmaxd,lmax,ki_start,ki_end,ki_step,nn_start,nn_end,lnonsph
     194             :    REAL,INTENT(IN)    :: qssbti(3),qssbtj(3)
     195             :    INTEGER,INTENT(IN) :: gvec(:,:,:),nv(2)
     196             :    REAL,INTENT(IN)    :: gk(:,:,:)
     197             :    REAL,INTENT(IN)    :: fleg1(0:lmaxd),fleg2(0:lmaxd),fl2p1(0:lmaxd)
     198             :    REAL,INTENT(IN)    :: fl2p1bt(0:lmaxd)
     199             :    REAL,INTENT(IN)    :: fj(:,0:,:),gj(:,0:,:)
     200             :    REAL,INTENT(IN)    :: taual(:,:)
     201             :    REAL,INTENT(IN)    :: ddn(0:lmaxd)
     202             :    REAL,INTENT(IN)    :: el(0:lmaxd)
     203             :    REAL,VALUE,INTENT(IN)    :: e_shift
     204             :    REAL,INTENT(INOUT) :: smat_data(:,:),hmat_data(:,:)
     205             :    !+APW
     206             :    REAL,INTENT(IN),OPTIONAL   :: uds(0:lmaxd),dus(0:lmaxd),us(0:lmaxd),duds(0:lmaxd)
     207             :    REAL,INTENT(IN),OPTIONAL   :: rmt
     208             :    !-APW
     209             : 
     210             :    REAL,   PARAMETER :: tpi_const=2.*3.1415926535897932
     211             :    REAL, ALLOCATABLE :: plegend(:)
     212             :    REAL cph
     213             :    REAL tnn(3), elall,fct,fct2,fjkiln,gjkiln,ddnln,ski(3)
     214             :    REAL apw_lo1,apw_lo2,apw1,w1
     215             :    INTEGER kii,ki,kj,l,nn,k
     216             :    INTEGER :: loop_start, loop_end, i
     217             : 
     218             :    ALLOCATE(plegend(0:lmaxd))
     219             :    plegend=0.0
     220             :    plegend(0)=1.0
     221             : 
     222             :    k = (blockidx%x-1)*blockdim%x + threadidx%x
     223             : 
     224             :    !TODO!!!
     225             :    !for seq, i.e. ki_start = 1, ki_step = 1
     226             :    loop_start = (k-1) * loop_size + 1
     227             :    loop_end = loop_start + loop_size - 1
     228             :    if (loop_end > ki_end ) loop_end = ki_end
     229             : 
     230             :    DO ki = loop_start,loop_end,ki_step
     231             :       !DO  ki =  ki_start,ki_end,ki_step
     232             :       DO kj = 1,ki
     233             :          kii=(ki-1)/ki_step+1
     234             :          ski = gvec(:,ki,jintsp) + qssbti
     235             :          !--->             set up phase factors
     236             :          cph = 0.0
     237             :          DO nn = nn_start,nn_end
     238             :             tnn = tpi_const*taual(:,nn)
     239             :             cph = cph + COS(DOT_PRODUCT(ski-gvec(:,kj,iintsp)-qssbtj,tnn))
     240             :             ! IF (iintsp.NE.jintsp) cph(kj)=CONJG(cph(kj))
     241             :          ENDDO
     242             :          !--->       legendre polynomials
     243             :          plegend(1) = DOT_PRODUCT(gk(:,kj,iintsp),gk(:,ki,jintsp))
     244             :          DO l = 1,lmax - 1
     245             :             plegend(l+1) = fleg1(l)*plegend(1)*plegend(l) - fleg2(l)*plegend(l-1)
     246             :          END DO
     247             :          DO  l = 0,lmax
     248             :             fjkiln = fj(ki,l,jintsp)
     249             :             gjkiln = gj(ki,l,jintsp)
     250             :             !+APW
     251             :             IF (PRESENT(uds)) THEN
     252             :                w1 = 0.5 * ( uds(l)*dus(l) + us(l)*duds(l) )
     253             :                apw_lo1 = fl2p1(l) * 0.5 * rmt**2 * ( gjkiln * w1 +&
     254             :                                                     fjkiln * us(l) * dus(l) )
     255             :                apw_lo2 = fl2p1(l) * 0.5 * rmt**2 * ( fjkiln * w1 +&
     256             :                                                     gjkiln * uds(l) * duds(l) )
     257             :             ENDIF
     258             :             !-APW
     259             :             ddnln =  ddn(l)
     260             :             elall = el(l)
     261             :             IF (l<=lnonsph) elall=elall-e_shift!(isp)
     262             :             !DO kj = 1,ki
     263             :             fct  = plegend(l)*fl2p1(l)*&
     264             :                    ( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln )
     265             :             fct2 = plegend(l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) )
     266             : 
     267             :             smat_data(kj,kii)=smat_data(kj,kii)+cph*fct
     268             :             hmat_data(kj,kii)=hmat_data(kj,kii) + cph * ( fct * elall + fct2)
     269             :             !+APW
     270             :             IF (PRESENT(uds)) THEN
     271             :                apw1 = cph * plegend(l)  * &
     272             :                       ( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) )
     273             :                hmat_data(kj,kii)=hmat_data(kj,kii) + apw1
     274             :             ENDIF
     275             :             !-APW
     276             :             !ENDDO
     277             :             !--->          end loop over l
     278             :          ENDDO
     279             :       ENDDO
     280             :       !--->    end loop over ki
     281             :    ENDDO
     282             :    DEALLOCATE(plegend)
     283             : END SUBROUTINE HsmtSphGpuKernel_real
     284             : 
     285             : ATTRIBUTES(global)&
     286             :    SUBROUTINE HsmtSphGpuKernel_cmplx(loop_size,iintsp,jintsp,nv,lmaxd,lmax,ki_start,ki_end,ki_step,nn_start,nn_end,&
     287             :                                      lnonsph,chi,qssbti,qssbtj,gvec,gk,fleg1,fleg2,fl2p1,fl2p1bt,fj,gj,taual,ddn,el,e_shift,&
     288             :                                      smat_data,hmat_data,&
     289             :                                      uds,dus,us,duds,rmt)
     290             : INTEGER, VALUE, INTENT(IN) :: loop_size
     291             : INTEGER, VALUE, INTENT(IN) :: iintsp,jintsp,lmaxd,lmax,ki_start,ki_end,ki_step,nn_start,nn_end,lnonsph
     292             : COMPLEX, VALUE, INTENT(IN)  :: chi
     293             : REAL,INTENT(IN)    :: qssbti(3),qssbtj(3)
     294             : INTEGER,INTENT(IN) :: gvec(:,:,:),nv(2)
     295             : REAL,INTENT(IN)    :: gk(:,:,:)
     296             : REAL,INTENT(IN)    :: fleg1(0:lmaxd),fleg2(0:lmaxd),fl2p1(0:lmaxd)
     297             : REAL,INTENT(IN)    :: fl2p1bt(0:lmaxd)
     298             : REAL,INTENT(IN)    :: fj(:,0:,:),gj(:,0:,:)
     299             : REAL,INTENT(IN)    :: taual(:,:)
     300             : REAL,INTENT(IN)    :: ddn(0:lmaxd)
     301             : REAL,INTENT(IN)    :: el(0:lmaxd)
     302             : REAL, VALUE,INTENT(IN)    :: e_shift
     303             : COMPLEX,INTENT(INOUT) :: smat_data(:,:),hmat_data(:,:)
     304             : !+APW
     305             : REAL,INTENT(IN),OPTIONAL   :: uds(0:lmaxd),dus(0:lmaxd),us(0:lmaxd),duds(0:lmaxd)
     306             : REAL,INTENT(IN),OPTIONAL   :: rmt
     307             : !-APW
     308             : 
     309             : REAL,   PARAMETER :: tpi_const=2.*3.1415926535897932
     310             : REAL, ALLOCATABLE :: plegend(:)
     311             : COMPLEX :: cph
     312             : REAL apw_lo1,apw_lo2,w1
     313             : COMPLEX capw1
     314             : REAL tnn(3), elall,fct,fct2,fjkiln,gjkiln,ddnln,ski(3)
     315             : INTEGER kii,ki,kj,kjj,l,nn,kj_end,k
     316             : INTEGER :: loop_start, loop_end, i
     317             : 
     318             : ALLOCATE(plegend(0:lmaxd))
     319             : plegend=0.0
     320             : plegend(0)=1.0
     321             : 
     322             : k = (blockidx%x-1)*blockdim%x + threadidx%x
     323             : 
     324             : !TODO!!!
     325             : !for seq, i.e. ki_start = 1, ki_step = 1
     326             : loop_start = (k-1) * loop_size + 1
     327             : loop_end = loop_start + loop_size - 1
     328             : if (loop_end > ki_end ) loop_end = ki_end
     329             : 
     330             : DO ki = loop_start,loop_end,ki_step
     331             :    !DO  ki =  ki_start,ki_end,ki_step
     332             :    kj_end = MIN(ki,nv(iintsp))
     333             :    DO kj = 1,kj_end
     334             :       kii=(ki-1)/ki_step+1
     335             :       ski = gvec(:,ki,jintsp) + qssbti
     336             : 
     337             :       !--->             set up phase factors
     338             :       cph = 0.0
     339             :       DO nn = nn_start,nn_end
     340             :          tnn = tpi_const*taual(:,nn)
     341             :          cph = cph +&
     342             :                CMPLX(COS(DOT_PRODUCT(ski-gvec(:,kj,iintsp)-qssbtj,tnn)),&
     343             :                      SIN(DOT_PRODUCT(gvec(:,kj,iintsp)+qssbtj-ski,tnn)))
     344             :          ! IF (iintsp.NE.jintsp) cph(kj)=CONJG(cph(kj))
     345             :       ENDDO
     346             :       !--->       legendre polynomials
     347             :       plegend(1) = DOT_PRODUCT(gk(:,kj,iintsp),gk(:,ki,jintsp))
     348             :       DO l = 1,lmax - 1
     349             :          plegend(l+1) = fleg1(l)*plegend(1)*plegend(l) - fleg2(l)*plegend(l-1)
     350             :       END DO
     351             :       DO  l = 0,lmax
     352             :          fjkiln = fj(ki,l,jintsp)
     353             :          gjkiln = gj(ki,l,jintsp)
     354             :          !+APW
     355             :          IF (PRESENT(uds)) THEN
     356             :             w1 = 0.5 * ( uds(l)*dus(l) + us(l)*duds(l) )
     357             :             apw_lo1 = fl2p1(l) * 0.5 * rmt**2 * ( gjkiln * w1 +&
     358             :                                                  fjkiln * us(l) * dus(l) )
     359             :             apw_lo2 = fl2p1(l) * 0.5 * rmt**2 * ( fjkiln * w1 +&
     360             :                                                  gjkiln * uds(l) * duds(l) )
     361             :          ENDIF
     362             :          !-APW
     363             :          ddnln =  ddn(l)
     364             :          elall = el(l)
     365             :          IF (l<=lnonsph) elall=elall-e_shift!(isp)
     366             :          fct  = plegend(l)*fl2p1(l)*&
     367             :                 ( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln )
     368             :          fct2 = plegend(l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) )
     369             : 
     370             :          smat_data(kj,kii)=smat_data(kj,kii) + chi * cph * fct
     371             :          hmat_data(kj,kii)=hmat_data(kj,kii) + chi * cph * ( fct * elall + fct2)
     372             :          !+APW
     373             :          IF (PRESENT(uds)) THEN
     374             :             capw1 = cph*plegend(l)&
     375             :                     * ( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) )
     376             :             hmat_data(kj,kii)=hmat_data(kj,kii) + capw1
     377             :          ENDIF
     378             :          !-APW
     379             :          !--->          end loop over l
     380             :       ENDDO
     381             :    ENDDO
     382             :    !--->    end loop over ki
     383             : ENDDO
     384             : DEALLOCATE(plegend)
     385             : END SUBROUTINE HsmtSphGpuKernel_cmplx
     386             : 
     387             : SUBROUTINE hsmt_sph_gpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_shift,usdus,fj,gj,smat,hmat)
     388             :    USE m_constants, ONLY : fpi_const,tpi_const
     389             :    USE m_types
     390             :    USE nvtx
     391             :    IMPLICIT NONE
     392             :    TYPE(t_input),INTENT(IN)      :: input
     393             :    TYPE(t_mpi),INTENT(IN)        :: mpi
     394             :    TYPE(t_noco),INTENT(IN)       :: noco
     395             :    TYPE(t_atoms),INTENT(IN)      :: atoms
     396             :    TYPE(t_lapw),INTENT(IN)       :: lapw
     397             :    TYPE(t_usdus),INTENT(IN)      :: usdus
     398             :    CLASS(t_mat),INTENT(INOUT)     :: smat,hmat
     399             :    !     ..
     400             :    !     .. Scalar Arguments ..
     401             :    INTEGER, INTENT (IN) :: n,isp,iintsp,jintsp
     402             :    COMPLEX, INTENT(IN)  :: chi
     403             :    !     ..
     404             :    !     .. Array Arguments ..
     405             :    REAL,MANAGED,INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,input%jspins)
     406             :    REAL,    INTENT (IN) :: e_shift!(atoms%ntype,input%jspins)
     407             :    REAL,MANAGED,INTENT(IN)    :: fj(:,0:,:),gj(:,0:,:)
     408             :    !     ..
     409             :    !     .. Local Scalars ..
     410             :    INTEGER l
     411             :    INTEGER :: grid, block, loop_size
     412             : 
     413             :    !     ..
     414             :    !     .. Local Arrays ..
     415             :    REAL,MANAGED :: fleg1(0:atoms%lmaxd),fleg2(0:atoms%lmaxd),fl2p1(0:atoms%lmaxd)
     416             :    REAL,MANAGED :: fl2p1bt(0:atoms%lmaxd)
     417             :    REAL,MANAGED :: qssbti(3),qssbtj(3)
     418             :    INTEGER, DEVICE :: nv_dev(2)
     419             : 
     420             :    call nvtxStartRange("hsmt_sph",2)
     421             :    CALL timestart("spherical setup")
     422             :    print*, "HsmtSph_GPU"
     423             : 
     424             :    DO l = 0,atoms%lmaxd
     425             :       fleg1(l) = REAL(l+l+1)/REAL(l+1)
     426             :       fleg2(l) = REAL(l)/REAL(l+1)
     427             :       fl2p1(l) = REAL(l+l+1)/fpi_const
     428             :       fl2p1bt(l) = fl2p1(l)*0.5
     429             :    END DO
     430             :    qssbti=MERGE(- noco%qss/2,+ noco%qss/2,jintsp.EQ.1)
     431             :    qssbtj=MERGE(- noco%qss/2,+ noco%qss/2,iintsp.EQ.1)
     432             : 
     433             :    ! pretty ugly solution
     434             :    nv_dev = lapw%nv
     435             :    loop_size = 1
     436             :    block = 32   ! number of threads in a block
     437             :    grid = ceiling(real(lapw%nv(jintsp))/(loop_SIZE*block))
     438             :    !loop_size = max(lapw%nv(jintsp)/(grid*block),1)   !number of iterations performed by each thread
     439             :    !if (loop_size * grid*block < lapw%nv(jintsp)) loop_size = loop_size + 1
     440             :    IF (input%l_useapw) THEN
     441             :       !TODO!!!!
     442             :       ! APW case is not testet
     443             :       IF (smat%l_real) THEN
     444             :          CALL HsmtSphGpuKernel_real<<<grid,block>>>(loop_size,iintsp,jintsp,nv_dev,&
     445             :                                                     atoms%lmaxd,atoms%lmax(n),mpi%n_rank+1,&
     446             :                                             lapw%nv(jintsp), mpi%n_size,SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n)),atoms%lnonsph(n),&
     447             :                                                     qssbti,qssbtj,lapw%gvec,lapw%gk,fleg1,fleg2,fl2p1,fl2p1bt,fj,gj,atoms%taual,&
     448             :                                                     usdus%ddn(:,n,isp),el(:,n,isp),e_shift,&
     449             :                                                     smat%data_r,hmat%data_r,&
     450             :                                            usdus%uds(:,n,isp),usdus%dus(:,n,isp),usdus%us(:,n,isp),usdus%duds(:,n,isp),atoms%rmt(n))
     451             :       ELSE
     452             :          CALL HsmtSphGpuKernel_cmplx<<<grid,block>>>(loop_size,iintsp,jintsp,nv_dev,&
     453             :                                                      atoms%lmaxd,atoms%lmax(n),mpi%n_rank+1,&
     454             :                                             lapw%nv(jintsp), mpi%n_size,SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n)),atoms%lnonsph(n),&
     455             :                                                    chi,qssbti,qssbtj,lapw%gvec,lapw%gk,fleg1,fleg2,fl2p1,fl2p1bt,fj,gj,atoms%taual,&
     456             :                                                      usdus%ddn(:,n,isp),el(:,n,isp),e_shift,&
     457             :                                                      smat%data_c,hmat%data_c,&
     458             :                                            usdus%uds(:,n,isp),usdus%dus(:,n,isp),usdus%us(:,n,isp),usdus%duds(:,n,isp),atoms%rmt(n))
     459             :       ENDIF
     460             :    ELSE
     461             :       IF (smat%l_real) THEN
     462             :          CALL HsmtSphGpuKernel_real<<<grid,block>>>(loop_size,iintsp,jintsp,nv_dev,&
     463             :                                                     atoms%lmaxd,atoms%lmax(n),mpi%n_rank+1,&
     464             :                                             lapw%nv(jintsp), mpi%n_size,SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n)),atoms%lnonsph(n),&
     465             :                                                     qssbti,qssbtj,lapw%gvec,lapw%gk,fleg1,fleg2,fl2p1,fl2p1bt,fj,gj,atoms%taual,&
     466             :                                                     usdus%ddn(:,n,isp),el(:,n,isp),e_shift,smat%data_r,hmat%data_r)
     467             :       ELSE
     468             :          CALL HsmtSphGpuKernel_cmplx<<<grid,block>>>(loop_size,iintsp,jintsp,nv_dev,&
     469             :                                                      atoms%lmaxd,atoms%lmax(n),mpi%n_rank+1,&
     470             :                                             lapw%nv(jintsp), mpi%n_size,SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n)),atoms%lnonsph(n),&
     471             :                                                    chi,qssbti,qssbtj,lapw%gvec,lapw%gk,fleg1,fleg2,fl2p1,fl2p1bt,fj,gj,atoms%taual,&
     472             :                                                      usdus%ddn(:,n,isp),el(:,n,isp),e_shift,smat%data_c,hmat%data_c)
     473             :       ENDIF
     474             :    ENDIF
     475             :    CALL timestop("spherical setup")
     476             : 
     477             :    call nvtxEndRange
     478             :    RETURN
     479             : END SUBROUTINE hsmt_sph_gpu
     480             : #endif
     481             : 
     482             : END MODULE m_hsmt_sph

Generated by: LCOV version 1.13