<commit_before>!
!  Ab-initio Real-Time Electron Dynamics Simulator, ARTED
!  Copyright (C) 2016  ARTED developers
!
!  This file is part of CG.f90.
!
!  CG.f90 is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!
!  CG.f90 is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with CG.f90.  If not, see <http://www.gnu.org/licenses/>.
!
!This file is "CG.f90"
!This file contain a subroutine.
!Subroutine CG
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine CG_omp(iter_cg_max)
  use Global_Variables
  use timelog
  implicit none
  real(8),parameter :: delta_cg=1.d-15
  integer iter,ik,ib,ibt
  integer :: iter_cg_max
  real(8) :: xkHxk,gkgk,pkHpk,xkTxk
  real(8) :: uk,s,ev
!  complex(8) :: xk(NL),hxk(NL),gk(NL),pk(NL),pko(NL),txk(NL)
  complex(8) :: cx,cp,xkHpk
  real(8) :: esp_var_l(1:NB,1:NK)
  complex(8) :: zs
! sato
  integer :: ia,j,i,ix,iy,iz
  real(8) :: kr
! omp
  integer :: thr_id,omp_get_thread_num
  thr_id=0

  call timelog_begin(LOG_CG)
  esp_var_l(:,:)=0.d0
!$omp parallel private(thr_id)
!$  thr_id=omp_get_thread_num()
!$omp do private(ia,j,i,ix,iy,iz,kr,ib,ibt,s,xkHxk,xkTxk,iter,uk,gkgk,xkHpk,pkHpk,ev,cx,cp,zs)
  do ik=NK_s,NK_e
! sato -----------------------------------------------------------------------------------------
!    k2=sum(kAc(ik,:)**2)
!Constructing nonlocal part
  do ia=1,NI
  do j=1,Mps(ia)
    i=Jxyz(j,ia); ix=Jxx(j,ia); iy=Jyy(j,ia); iz=Jzz(j,ia)
    kr=kAc(ik,1)*(Lx(i)*Hx-ix*aLx)+kAc(ik,2)*(Ly(i)*Hy-iy*aLy)+kAc(ik,3)*(Lz(i)*Hz-iz*aLz)
    ekr_omp(j,ia,ik)=exp(zI*kr)
  enddo
  enddo
! sato -----------------------------------------------------------------------------------------
  do ib=1,NB
    do ibt=1,ib-1
      s=sum(conjg(zu_GS(:,ibt,ik))*zu_GS(:,ib,ik))*Hxyz
      zu_GS(1:NL,ib,ik)=zu_GS(1:NL,ib,ik)-zu_GS(1:NL,ibt,ik)*s
    end do
    s=1.0d0/sqrt(sum(abs(zu_GS(:,ib,ik))**2)*Hxyz)
    xk_omp(1:NL,thr_id)=zu_GS(1:NL,ib,ik)*s
    call hpsi_omp_KB(ik,xk_omp(:,thr_id),txk_omp(:,thr_id),hxk_omp(:,thr_id))
    xkHxk=sum(conjg(xk_omp(:,thr_id))*hxk_omp(:,thr_id))*Hxyz
    xkTxk=sum(conjg(xk_omp(:,thr_id))*txk_omp(:,thr_id))*Hxyz

    do iter=1,iter_cg_max
      gk_omp(1:NL,thr_id)=(hxk_omp(1:NL,thr_id)-xkHxk*xk_omp(1:NL,thr_id))
      do ibt=1,ib-1
        zs=sum(conjg(zu_GS(:,ibt,ik))*gk_omp(:,thr_id))*Hxyz
        gk_omp(1:NL,thr_id)=gk_omp(1:NL,thr_id)-zu_GS(1:NL,ibt,ik)*zs
      end do
      s=sum(abs(gk_omp(:,thr_id))**2)*Hxyz

      select case (iter)
      case(1)
        pk_omp(1:NL,thr_id)=gk_omp(1:NL,thr_id)
      case default
        uk=s/gkgk
        pk_omp(1:NL,thr_id)=gk_omp(1:NL,thr_id)+uk*pk_omp(1:NL,thr_id)
      end select
      gkgk=s

      zs=sum(conjg(xk_omp(:,thr_id))*pk_omp(:,thr_id))*Hxyz
      pko_omp(1:NL,thr_id)=pk_omp(1:NL,thr_id)-xk_omp(1:NL,thr_id)*zs
      s=1.0d0/sqrt(sum(abs(pko_omp(:,thr_id))**2)*Hxyz)
      pko_omp(1:NL,thr_id)=pko_omp(1:NL,thr_id)*s
      call hpsi_omp_KB(ik,pko_omp(:,thr_id),ttpsi_omp(:,thr_id),htpsi_omp(:,thr_id))
      xkHpk=sum(conjg( xk_omp(:,thr_id))*htpsi_omp(:,thr_id))*Hxyz
      pkHpk=sum(conjg(pko_omp(:,thr_id))*htpsi_omp(:,thr_id))*Hxyz
      ev=0.5d0*((xkHxk+pkHpk)-sqrt((xkHxk-pkHpk)**2+4*abs(xkHpk)**2))
      cx=xkHpk/(ev-xkHxk)
      cp=1.d0/sqrt(1.d0+abs(cx)**2)
      cx=cx*cp
      if(abs(ev-xkHxk)<delta_cg) exit
       xk_omp(1:NL,thr_id)=cx* xk_omp(1:NL,thr_id)+cp*  pko_omp(1:NL,thr_id)
      hxk_omp(1:NL,thr_id)=cx*hxk_omp(1:NL,thr_id)+cp*htpsi_omp(1:NL,thr_id)
      txk_omp(1:NL,thr_id)=cx*txk_omp(1:NL,thr_id)+cp*ttpsi_omp(1:NL,thr_id)
      xkHxk=sum(conjg(xk_omp(:,thr_id))*hxk_omp(:,thr_id))*Hxyz
      xkTxk=sum(conjg(xk_omp(:,thr_id))*txk_omp(:,thr_id))*Hxyz
    enddo

    s=1.0d0/sqrt(sum(abs(xk_omp(:,thr_id))**2)*Hxyz)
    zu_GS(1:NL,ib,ik)=xk_omp(1:NL,thr_id)*s
    call hpsi_omp_KB(ik,zu_GS(:,ib,ik),ttpsi_omp(:,thr_id),htpsi_omp(:,thr_id))
    xkHxk=sum(conjg(tpsi_omp(:,thr_id))*htpsi_omp(:,thr_id))*Hxyz
    esp_var_l(ib,ik)=sqrt(sum(abs(htpsi_omp(:,thr_id)-xkHxk*tpsi_omp(:,thr_id))**2)*Hxyz)*occ(ib,ik)
  enddo
  enddo

!$omp end parallel

  call timelog_begin(LOG_ALLREDUCE)
  call MPI_ALLREDUCE(esp_var_l,esp_var,NB*NK,MPI_REAL8,MPI_SUM,NEW_COMM_WORLD,ierr)
  call timelog_end(LOG_ALLREDUCE)

  call timelog_end(LOG_CG)

  return
End Subroutine CG_OMP
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130<commit_msg>A bug in CG was fixed.<commit_after>!
!  Ab-initio Real-Time Electron Dynamics Simulator, ARTED
!  Copyright (C) 2016  ARTED developers
!
!  This file is part of CG.f90.
!
!  CG.f90 is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!
!  CG.f90 is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with CG.f90.  If not, see <http://www.gnu.org/licenses/>.
!
!This file is "CG.f90"
!This file contain a subroutine.
!Subroutine CG
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine CG_omp(iter_cg_max)
  use Global_Variables
  use timelog
  implicit none
  real(8),parameter :: delta_cg=1.d-15
  integer iter,ik,ib,ibt
  integer :: iter_cg_max
  real(8) :: xkHxk,gkgk,pkHpk,xkTxk
  real(8) :: uk,s,ev
!  complex(8) :: xk(NL),hxk(NL),gk(NL),pk(NL),pko(NL),txk(NL)
  complex(8) :: cx,cp,xkHpk
  real(8) :: esp_var_l(1:NB,1:NK)
  complex(8) :: zs
! sato
  integer :: ia,j,i,ix,iy,iz
  real(8) :: kr
! omp
  integer :: thr_id,omp_get_thread_num
  thr_id=0

  call timelog_begin(LOG_CG)
  esp_var_l(:,:)=0.d0
!$omp parallel private(thr_id)
!$  thr_id=omp_get_thread_num()
!$omp do private(ia,j,i,ix,iy,iz,kr,ib,ibt,s,xkHxk,xkTxk,iter,uk,gkgk,xkHpk,pkHpk,ev,cx,cp,zs)
  do ik=NK_s,NK_e
! sato -----------------------------------------------------------------------------------------
!    k2=sum(kAc(ik,:)**2)
!Constructing nonlocal part
  do ia=1,NI
  do j=1,Mps(ia)
    i=Jxyz(j,ia); ix=Jxx(j,ia); iy=Jyy(j,ia); iz=Jzz(j,ia)
    kr=kAc(ik,1)*(Lx(i)*Hx-ix*aLx)+kAc(ik,2)*(Ly(i)*Hy-iy*aLy)+kAc(ik,3)*(Lz(i)*Hz-iz*aLz)
    ekr_omp(j,ia,ik)=exp(zI*kr)
  enddo
  enddo
! sato -----------------------------------------------------------------------------------------
  do ib=1,NB
    do ibt=1,ib-1
      s=sum(conjg(zu_GS(:,ibt,ik))*zu_GS(:,ib,ik))*Hxyz
      zu_GS(1:NL,ib,ik)=zu_GS(1:NL,ib,ik)-zu_GS(1:NL,ibt,ik)*s
    end do
    s=1.0d0/sqrt(sum(abs(zu_GS(:,ib,ik))**2)*Hxyz)
    xk_omp(1:NL,thr_id)=zu_GS(1:NL,ib,ik)*s
    call hpsi_omp_KB(ik,xk_omp(:,thr_id),txk_omp(:,thr_id),hxk_omp(:,thr_id))
    xkHxk=sum(conjg(xk_omp(:,thr_id))*hxk_omp(:,thr_id))*Hxyz
    xkTxk=sum(conjg(xk_omp(:,thr_id))*txk_omp(:,thr_id))*Hxyz

    do iter=1,iter_cg_max
      gk_omp(1:NL,thr_id)=(hxk_omp(1:NL,thr_id)-xkHxk*xk_omp(1:NL,thr_id))
      do ibt=1,ib-1
        zs=sum(conjg(zu_GS(:,ibt,ik))*gk_omp(:,thr_id))*Hxyz
        gk_omp(1:NL,thr_id)=gk_omp(1:NL,thr_id)-zu_GS(1:NL,ibt,ik)*zs
      end do
      s=sum(abs(gk_omp(:,thr_id))**2)*Hxyz

      select case (iter)
      case(1)
        pk_omp(1:NL,thr_id)=gk_omp(1:NL,thr_id)
      case default
        uk=s/gkgk
        pk_omp(1:NL,thr_id)=gk_omp(1:NL,thr_id)+uk*pk_omp(1:NL,thr_id)
      end select
      gkgk=s

      zs=sum(conjg(xk_omp(:,thr_id))*pk_omp(:,thr_id))*Hxyz
      pko_omp(1:NL,thr_id)=pk_omp(1:NL,thr_id)-xk_omp(1:NL,thr_id)*zs
      s=1.0d0/sqrt(sum(abs(pko_omp(:,thr_id))**2)*Hxyz)
      pko_omp(1:NL,thr_id)=pko_omp(1:NL,thr_id)*s
      call hpsi_omp_KB(ik,pko_omp(:,thr_id),ttpsi_omp(:,thr_id),htpsi_omp(:,thr_id))
      xkHpk=sum(conjg( xk_omp(:,thr_id))*htpsi_omp(:,thr_id))*Hxyz
      pkHpk=sum(conjg(pko_omp(:,thr_id))*htpsi_omp(:,thr_id))*Hxyz
      ev=0.5d0*((xkHxk+pkHpk)-sqrt((xkHxk-pkHpk)**2+4*abs(xkHpk)**2))
      cx=xkHpk/(ev-xkHxk)
      cp=1.d0/sqrt(1.d0+abs(cx)**2)
      cx=cx*cp
      if(abs(ev-xkHxk)<delta_cg) exit
       xk_omp(1:NL,thr_id)=cx* xk_omp(1:NL,thr_id)+cp*  pko_omp(1:NL,thr_id)
      hxk_omp(1:NL,thr_id)=cx*hxk_omp(1:NL,thr_id)+cp*htpsi_omp(1:NL,thr_id)
      txk_omp(1:NL,thr_id)=cx*txk_omp(1:NL,thr_id)+cp*ttpsi_omp(1:NL,thr_id)
      xkHxk=sum(conjg(xk_omp(:,thr_id))*hxk_omp(:,thr_id))*Hxyz
      xkTxk=sum(conjg(xk_omp(:,thr_id))*txk_omp(:,thr_id))*Hxyz
    enddo

    s=1.0d0/sqrt(sum(abs(xk_omp(:,thr_id))**2)*Hxyz)
    zu_GS(1:NL,ib,ik)=xk_omp(1:NL,thr_id)*s
    call hpsi_omp_KB(ik,zu_GS(:,ib,ik),ttpsi_omp(:,thr_id),htpsi_omp(:,thr_id))
    xkHxk=sum(conjg(zu_GS(1:NL,ib,ik))*htpsi_omp(:,thr_id))*Hxyz
    esp_var_l(ib,ik)=sqrt(sum(abs(htpsi_omp(:,thr_id)-xkHxk*zu_GS(1:NL,ib,ik))**2)*Hxyz)*occ(ib,ik)
  enddo
  enddo

!$omp end parallel

  call timelog_begin(LOG_ALLREDUCE)
  call MPI_ALLREDUCE(esp_var_l,esp_var,NB*NK,MPI_REAL8,MPI_SUM,NEW_COMM_WORLD,ierr)
  call timelog_end(LOG_ALLREDUCE)

  call timelog_end(LOG_CG)

  return
End Subroutine CG_OMP
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
{'fraction_non_alphanumeric': 0.24458456218854877, 'fraction_numerical': 0.03335706295128649, 'mean_word_length': 3.411843876177658, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5536491', 'n_tokens_mistral': 4939, 'n_tokens_neox': 4701, 'n_words': 552}