LCOV - code coverage report
Current view: top level - src/79_seqpar_mpi - m_lobpcgwf_old.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 70.1 % 622 436
Test Date: 2026-11-08 08:00:28 Functions: 100.0 % 1 1
Branches: 57.9 % 1686 976

             Branch data     Line data    Source code
       1                 :             : !!****m* ABINIT/m_lobpcgwf_old
       2                 :             : !! NAME
       3                 :             : !!   m_lobpcgwf_old
       4                 :             : !!
       5                 :             : !! FUNCTION
       6                 :             : !!
       7                 :             : !!
       8                 :             : !! COPYRIGHT
       9                 :             : !!  Copyright (C) 2008-2026 ABINIT group ()
      10                 :             : !!  This file is distributed under the terms of the
      11                 :             : !!  GNU General Public License, see ~abinit/COPYING
      12                 :             : !!  or http://www.gnu.org/copyleft/gpl.txt .
      13                 :             : !!
      14                 :             : !! SOURCE
      15                 :             : 
      16                 :             : #if defined HAVE_CONFIG_H
      17                 :             : #include "config.h"
      18                 :             : #endif
      19                 :             : 
      20                 :             : #include "abi_common.h"
      21                 :             : 
      22                 :             : module m_lobpcgwf_old
      23                 :             : 
      24                 :             :  implicit none
      25                 :             : 
      26                 :             :  private
      27                 :             : !!***
      28                 :             : 
      29                 :             :  public :: lobpcgwf
      30                 :             : !!***
      31                 :             : 
      32                 :             : contains
      33                 :             : !!***
      34                 :             : 
      35                 :             : !!****f* ABINIT/lobpcgwf
      36                 :             : !! NAME
      37                 :             : !! lobpcgwf
      38                 :             : !!
      39                 :             : !! FUNCTION
      40                 :             : !! this routine updates the whole wave functions at a given k-point,
      41                 :             : !! using the lobpcg method
      42                 :             : !! for a given spin-polarization, from a fixed hamiltonian
      43                 :             : !! but might also simply compute eigenvectors and eigenvalues at this k point.
      44                 :             : !! it will also update the matrix elements of the hamiltonian.
      45                 :             : !!
      46                 :             : !! COPYRIGHT
      47                 :             : !! Copyright (C) 1998-2026 ABINIT group (FBottin,GZ,AR,MT,FDahm)
      48                 :             : !! this file is distributed under the terms of the
      49                 :             : !! gnu general public license, see ~abinit/COPYING
      50                 :             : !! or http://www.gnu.org/copyleft/gpl.txt .
      51                 :             : !! for the initials of contributors, see ~abinit/doc/developers/contributors.txt .
      52                 :             : !!
      53                 :             : !! INPUTS
      54                 :             : !!  dtset <type(dataset_type)>=all input variales for this dataset
      55                 :             : !!  gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k
      56                 :             : !!  icg=shift to be applied on the location of data in the array cg
      57                 :             : !!  igsc=shift to be applied on the location of data in the array gsc
      58                 :             : !!  kinpw(npw)=(modified) kinetic energy for each plane wave (hartree)
      59                 :             : !!  kinpw_kp(npw)=GBT kinetic energy for the second spinor component
      60                 :             : !!  mcg=second dimension of the cg array
      61                 :             : !!  mgsc=second dimension of the gsc array
      62                 :             : !!  mpi_enreg=information about MPI parallelization
      63                 :             : !!  nband_k=number of bands at this k point for that spin polarization
      64                 :             : !!  nbdblock : number of blocks
      65                 :             : !!  npw_k=number of plane waves at this k point
      66                 :             : !!  prtvol=control print volume and debugging output
      67                 :             : !!  use_totvnlx=1 if one has to compute totvnlx
      68                 :             : !!
      69                 :             : !! OUTPUT
      70                 :             : !!  resid_k(nband_k)=residuals for each states
      71                 :             : !!  subham(nband_k*(nband_k+1))=the matrix elements of h
      72                 :             : !!  If gs_hamk%usepaw==0:
      73                 :             : !!    gsc(2,mgsc)=<g|s|c> matrix elements (s=overlap)
      74                 :             : !!    totvnlx(nband_k*use_totvnlx,nband_k*use_totvnlx)=the matrix elements of vnl+vfockACE
      75                 :             : !!
      76                 :             : !! SIDE EFFECTS
      77                 :             : !!  cg(2,mcg)=updated wavefunctions
      78                 :             : !!
      79                 :             : !! SOURCE
      80                 :             : 
      81         [ -  + ]:        2250 : subroutine lobpcgwf(cg,dtset,gs_hamk,gsc,icg,igsc,kinpw,kinpw_kp,mcg,mgsc,mpi_enreg,&
      82                 :        2250 : &                   nband_k,nbdblock,npw_k,prtvol,resid_k,subham,totvnlx,use_totvnlx)
      83                 :             : 
      84                 :             : 
      85                 :             :  use, intrinsic :: iso_c_binding
      86                 :             :  use defs_basis
      87                 :             :  use m_abicore
      88                 :             :  use m_lobpcg
      89                 :             :  use m_gputk
      90                 :             :  use m_abi_linalg
      91                 :             :  use m_wfutils
      92                 :             :  use m_xmpi
      93                 :             :  use m_errors
      94                 :             :  use m_dtset
      95                 :             : 
      96                 :             :  use defs_abitypes, only : mpi_type
      97                 :             :  use m_time,        only : timab
      98                 :             :  use m_hamiltonian, only : gs_hamiltonian_type
      99                 :             :  use m_pawcprj,     only : pawcprj_type
     100                 :             :  use m_getghc,      only : getghc
     101                 :             :  use m_prep_kgb,    only : prep_getghc
     102                 :             : #ifdef HAVE_GPU
     103                 :             :  use m_gputk
     104                 :             : #endif
     105                 :             : 
     106                 :             : !Arguments ------------------------------------
     107                 :             :  integer,intent(in) :: icg,igsc,mcg,mgsc,nband_k,nbdblock,npw_k,prtvol,use_totvnlx
     108                 :             :  type(gs_hamiltonian_type),intent(inout) :: gs_hamk
     109                 :             :  type(dataset_type),intent(in) :: dtset
     110                 :             :  type(mpi_type),intent(in) :: mpi_enreg
     111                 :             :  real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc)
     112                 :             :  real(dp),intent(in) :: kinpw(npw_k),kinpw_kp(:)
     113                 :             :  real(dp),intent(out) :: resid_k(nband_k)
     114                 :             :  real(dp),intent(inout) :: subham(nband_k*(nband_k+1))
     115                 :             :  real(dp),intent(inout) :: totvnlx((3-gs_hamk%istwf_k)*nband_k*use_totvnlx,nband_k*use_totvnlx)
     116                 :             : 
     117                 :             : !Local variables-------------------------------
     118                 :             :  integer, parameter :: tim_getghc=5
     119                 :             :  integer :: activepsize,activersize,bblocksize,bigorder,blocksize,cpopt
     120                 :             :  integer :: cond_try
     121                 :             :  integer :: iblocksize,iblock,ierr,ii,info,istwf_k,isubh
     122                 :             :  integer :: iterationnumber
     123                 :             :  integer :: iwavef,i1,i2,i3,i4,maxiterations,my_nspinor
     124                 :             :  integer :: nrestart,optekin,optpcon,restart
     125                 :             :  integer :: sij_opt,timopt,tim_wfcopy,tim_xeigen
     126                 :             :  integer :: tim_xortho,tim_xprecon,use_lapack_gpu,use_linalg_gpu,vectsize
     127                 :             :  logical :: gen_eigenpb
     128                 :             :  integer :: cplx
     129                 :             :  real(dp) :: condestgramb,deltae,deold,dum
     130                 :             :  complex(dp) :: cminusone
     131                 :             :  real(dp) :: zvar(2)
     132                 :             :  logical :: havetoprecon
     133                 :             :  real(dp) :: tsec(2)
     134                 :        2250 :  real(dp), allocatable :: gwavef(:,:),cwavef(:,:),gvnlxc(:,:)
     135                 :        2250 :  real(dp), allocatable :: swavef(:,:)
     136                 :        2250 :  real(dp), allocatable :: residualnorms(:),eigen(:)
     137                 :        2250 :  real(dp), allocatable :: tmpeigen(:)
     138                 :        2250 :  real(dp), allocatable :: pcon(:,:),pcon_kp(:,:)
     139                 :        2250 :  real(dp), allocatable, target :: blockvectorx(:,:),blockvectorvx(:,:),blockvectorax(:,:),blockvectorbx(:,:)
     140                 :        2250 :  real(dp), allocatable, target :: blockvectorr(:,:),blockvectorvr(:,:),blockvectorar(:,:),blockvectorbr(:,:)
     141                 :        2250 :  real(dp), allocatable, target :: blockvectorp(:,:),blockvectorvp(:,:),blockvectorap(:,:),blockvectorbp(:,:),blockvectordumm(:,:)
     142                 :        2250 :  real(dp), allocatable, target :: blockvectory(:,:),blockvectorby(:,:),blockvectorz(:,:)
     143                 :        2250 :  real(dp), allocatable, target :: gramxax(:,:),gramxar(:,:),gramxap(:,:),gramrar(:,:),gramrap(:,:),grampap(:,:)
     144                 :        2250 :  real(dp), allocatable, target :: gramxbx(:,:),gramxbr(:,:),gramxbp(:,:),gramrbr(:,:),gramrbp(:,:),grampbp(:,:)
     145                 :        2250 :  real(dp), allocatable, target :: coordx1(:,:),coordx2(:,:),coordx3(:,:),lambda(:,:),grama(:,:),gramb(:,:),gramyx(:,:)
     146                 :        2250 :  real(dp), allocatable :: tmpgramb(:,:),transf3(:,:,:),transf5(:,:,:)
     147                 :        2250 :  real(dp), allocatable :: tsubham(:,:)
     148                 :        2250 :  type(pawcprj_type), allocatable :: cprj_dum(:,:)
     149                 :             :  character(len=500) :: message
     150                 :             :  character, dimension(2) :: cparam
     151                 :             :  type(c_ptr) :: A_gpu,C_gpu,coordx2_gpu,coordx3_gpu,bblockvector_gpu,gram_gpu
     152                 :             :  type(c_ptr) :: blockvectorr_gpu,blockvectorar_gpu,blockvectorbr_gpu
     153                 :             : 
     154                 :             : !Index of a given band
     155                 :             : !gramindex(iblocksize)=(iblocksize-1)*cplx+1
     156                 :             : 
     157                 :             : ! *********************************************************************
     158                 :             : 
     159                 :             :  DBG_ENTER("COLL")
     160                 :             : 
     161                 :        2250 :  call timab(530,1,tsec)
     162         [ -  + ]:        2250 :  if(abs(dtset%timopt)==4) then
     163                 :           0 :    call timab(520,1,tsec)
     164                 :             :  end if
     165                 :             : 
     166                 :             : !###########################################################################
     167                 :             : !################ INITIALISATION  ##########################################
     168                 :             : !###########################################################################
     169                 :             : 
     170                 :             : !For timing
     171                 :        2250 :  timopt=dtset%timopt
     172                 :        2250 :  tim_wfcopy=584
     173                 :             :  !tim_xcopy=584
     174                 :        2250 :  tim_xeigen=587
     175                 :             :  !tim_xgemm=532
     176                 :        2250 :  tim_xortho=535
     177                 :        2250 :  tim_xprecon=536
     178                 :             :  !tim_xtrsm=535
     179                 :             : 
     180                 :             : !Variables
     181                 :        2250 :  maxiterations=dtset%nline
     182                 :        2250 :  gen_eigenpb=(gs_hamk%usepaw==1)
     183                 :        2250 :  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
     184                 :        2250 :  cminusone=-cone
     185                 :        2250 :  istwf_k=gs_hamk%istwf_k
     186                 :        2250 :  info = 0
     187                 :        2250 :  cparam(1)='t'
     188                 :        2250 :  cparam(2)='c'
     189                 :             : 
     190                 :             : !Depends on istwfk
     191         [ +  + ]:        2250 :  if ( istwf_k == 2 ) then
     192                 :         182 :    cplx=1
     193         [ +  + ]:         182 :    if (mpi_enreg%me_g0 == 1) then
     194                 :          83 :      vectsize=2*npw_k*my_nspinor-1
     195                 :             :    else
     196                 :          99 :      vectsize=2*npw_k*my_nspinor
     197                 :             :    end if
     198                 :             :  else
     199                 :        2068 :    cplx=2
     200                 :        2068 :    vectsize=npw_k*my_nspinor
     201                 :             :  end if
     202                 :             : 
     203                 :             : !For preconditionning
     204                 :        2250 :  optekin=0;if (dtset%wfoptalg>10) optekin=0
     205         [ +  + ]:        2250 :  optpcon=1;if (dtset%wfoptalg>10) optpcon=0
     206                 :             : 
     207                 :             : !For communication
     208                 :             :  !blocksize=mpi_enreg%nproc_fft
     209                 :             :  !if(mpi_enreg%paral_kgb==1) blocksize=mpi_enreg%nproc_band*mpi_enreg%bandpp
     210                 :        2250 :  blocksize=nband_k/dtset%nblock_lobpcg
     211                 :             :  !IF you want to compare with new lobpcg in sequential uncomment the following
     212                 :             :  !line
     213                 :             :  !blocksize=mpi_enreg%nproc_band*mpi_enreg%bandpp
     214                 :             : 
     215                 :             : !Iniitializations/allocations of GPU parallelism
     216                 :        2250 :  use_linalg_gpu=0;use_lapack_gpu=0
     217   [ +  -  -  - ]:        2250 :  if ((dtset%gpu_option==ABI_GPU_LEGACY).and. &
     218                 :        2250 : & (vectsize*blocksize*blocksize>dtset%gpu_linalg_limit)) use_linalg_gpu=1
     219         [ +  - ]:        2250 :  if (dtset%gpu_option==ABI_GPU_OPENMP) use_linalg_gpu=1
     220                 :             : #ifdef HAVE_GPU_HIP
     221                 :             :  use_linalg_gpu=0
     222                 :             : #endif
     223                 :             : #if defined HAVE_LINALG_MAGMA
     224                 :             :  use_lapack_gpu=use_linalg_gpu
     225                 :             : #endif
     226         [ -  + ]:        2250 :  if(use_linalg_gpu==1) then
     227                 :           0 :    call alloc_on_gpu(A_gpu,             INT(cplx, c_size_t)*dp*vectsize*blocksize)
     228                 :           0 :    call alloc_on_gpu(C_gpu,             INT(cplx, c_size_t)*dp*vectsize*blocksize)
     229                 :           0 :    call alloc_on_gpu(blockvectorr_gpu,  INT(cplx, c_size_t)*dp*vectsize*blocksize)
     230                 :           0 :    call alloc_on_gpu(blockvectorar_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     231                 :           0 :    call alloc_on_gpu(blockvectorbr_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     232                 :           0 :    call alloc_on_gpu(coordx2_gpu,       INT(cplx, c_size_t)*dp*blocksize*blocksize)
     233                 :           0 :    call alloc_on_gpu(coordx3_gpu,       INT(cplx, c_size_t)*dp*blocksize*blocksize)
     234                 :             :  end if
     235                 :             : 
     236   [ -  +  +  -  :       13658 :  ABI_MALLOC(cprj_dum, (gs_hamk%natom, 1))
          +  -  +  -  -  
          +  -  +  -  +  
          +  +  +  +  +  
                      - ]
     237                 :             : 
     238                 :             :  ! Work arrays eventually mapped on GPU via OpenMP
     239   [ -  +  +  -  :        6750 :  ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     240   [ -  +  +  -  :        4500 :  ABI_MALLOC(gwavef,(2,npw_k*my_nspinor*blocksize))
          -  +  -  +  -  
                      + ]
     241   [ -  +  +  -  :        4500 :  ABI_MALLOC(gvnlxc,(2,npw_k*my_nspinor*blocksize))
          -  +  -  +  -  
                      + ]
     242   [ -  +  +  -  :        4500 :  ABI_MALLOC(swavef,(2,npw_k*my_nspinor*blocksize))
          -  +  -  +  -  
                      + ]
     243                 :             : #ifdef HAVE_OPENMP_OFFLOAD
     244                 :             :  !$OMP TARGET ENTER DATA MAP(alloc:cwavef,gwavef,gvnlxc,swavef) IF(dtset%gpu_option==ABI_GPU_OPENMP)
     245                 :             : #endif
     246                 :             : 
     247         [ -  + ]:        2250 :  if(abs(dtset%timopt)==4) then
     248                 :           0 :    call timab(520,2,tsec)
     249                 :             :  end if
     250                 :             : 
     251                 :             : !###########################################################################
     252                 :             : !################ BIG LOOP OVER BLOCKS  ####################################
     253                 :             : !###########################################################################
     254                 :             : 
     255         [ +  + ]:        7886 :  do iblock=1,nbdblock
     256                 :             : 
     257         [ -  + ]:        5636 :    if(abs(dtset%timopt)==4) then
     258                 :           0 :      call timab(521,1,tsec)
     259                 :             :    end if
     260                 :             : 
     261                 :        5636 :    havetoprecon=.true.
     262                 :        5636 :    nrestart=0
     263                 :        5636 :    bblocksize=(iblock-1)*blocksize
     264                 :             : 
     265                 :             : !  allocations
     266   [ -  +  +  -  :       22544 :    ABI_MALLOC(pcon,(npw_k,blocksize))
          -  +  +  -  +  
          -  +  -  -  +  
                   -  + ]
     267         [ -  + ]:        5636 :    if (gs_hamk%use_gbt /= 0) then
     268   [ #  #  #  #  :           0 :      ABI_MALLOC(pcon_kp,(npw_k,blocksize))
          #  #  #  #  #  
             #  #  #  #  
                      # ]
     269                 :             :    end if
     270   [ -  +  +  -  :       22544 :    ABI_MALLOC(blockvectorx,(cplx*vectsize,blocksize))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     271   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorax,(cplx*vectsize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     272   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorbx,(cplx*vectsize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     273   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorr,(cplx*vectsize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     274   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorar,(cplx*vectsize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     275   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorbr,(cplx*vectsize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     276   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorp,(cplx*vectsize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     277   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorap,(cplx*vectsize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     278   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorbp,(cplx*vectsize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     279   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectordumm,(cplx*vectsize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     280   [ -  +  +  -  :       22544 :    ABI_MALLOC(gramxax,(cplx*blocksize,blocksize))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     281   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramxar,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     282   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramxap,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     283   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramrar,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     284   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramrap,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     285   [ -  +  +  -  :       16908 :    ABI_MALLOC(grampap,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     286   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramxbx,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     287   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramxbr,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     288   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramxbp,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     289   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramrbr,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     290   [ -  +  +  -  :       16908 :    ABI_MALLOC(gramrbp,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     291   [ -  +  +  -  :       16908 :    ABI_MALLOC(grampbp,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     292   [ -  +  +  -  :       28180 :    ABI_MALLOC(transf3,(cplx*blocksize,blocksize,3))
          -  +  +  -  +  
          -  +  -  +  -  
          -  +  -  +  -  
                      + ]
     293   [ -  +  +  -  :       28180 :    ABI_MALLOC(transf5,(cplx*blocksize,blocksize,5))
          -  +  +  -  +  
          -  +  -  +  -  
          -  +  -  +  -  
                      + ]
     294   [ -  +  +  -  :       16908 :    ABI_MALLOC(lambda,(cplx*blocksize,blocksize))
          -  +  +  -  -  
                +  -  + ]
     295   [ -  +  +  -  :       16908 :    ABI_MALLOC(residualnorms,(blocksize))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     296                 :             : 
     297   [ -  +  +  -  :       22544 :    ABI_MALLOC(blockvectory,(cplx*vectsize,bblocksize))
          +  +  +  -  +  
          -  +  +  -  +  
             -  +  -  + ]
     298   [ -  +  +  -  :       16908 :    ABI_MALLOC(blockvectorby,(cplx*vectsize,bblocksize))
          +  +  +  -  -  
             +  -  +  -  
                      + ]
     299   [ +  +  +  -  :       22544 :    ABI_MALLOC(gramyx,(cplx*bblocksize,blocksize))
          -  +  +  -  +  
          -  +  +  -  +  
             -  +  -  + ]
     300         [ +  + ]:        5636 :    if (gs_hamk%usepaw==0) then
     301   [ -  +  +  -  :        2496 :      ABI_MALLOC(blockvectorvx,(cplx*vectsize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     302   [ -  +  +  -  :        2496 :      ABI_MALLOC(blockvectorvr,(cplx*vectsize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     303   [ -  +  +  -  :        2496 :      ABI_MALLOC(blockvectorvp,(cplx*vectsize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     304                 :             :    end if
     305                 :             : 
     306         [ -  + ]:        5636 :    if(use_linalg_gpu==1) then
     307         [ #  # ]:           0 :      if(iblock/=1) then
     308                 :           0 :        call alloc_on_gpu(bblockvector_gpu, INT(cplx, c_size_t)*dp*vectsize*bblocksize)
     309                 :           0 :        call alloc_on_gpu(gram_gpu,         INT(cplx, c_size_t)*dp*bblocksize*blocksize)
     310                 :             :      else
     311                 :           0 :        call alloc_on_gpu(bblockvector_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     312                 :           0 :        call alloc_on_gpu(gram_gpu,         INT(cplx, c_size_t)*dp*blocksize*blocksize)
     313                 :             :      end if
     314                 :             :    end if
     315                 :             : 
     316                 :             : ! Initialize global variables in m_wfutils.
     317                 :        5636 :    call setWFParameter(cplx,mpi_enreg%me_g0,npw_k,my_nspinor,icg,igsc,blocksize)
     318                 :             : 
     319                 :             : !  transfer array of wf coeff in iblock to blockvectorx
     320                 :             :    call wfcopy('D',blocksize*vectsize,cg,1,blockvectorx,1,blocksize,iblock,'C',withbbloc=.true.,&
     321                 :        5636 : &   timopt=timopt,tim_wfcopy=tim_wfcopy)
     322                 :             : 
     323                 :             : !  !!!!!!!!!!!!!!!!!!!!!!!! Begin if iblock /=1 !!!!!!!!!!!!!!!!!!!!!!!!!!
     324                 :             : !  transfer array of wf coeff less than iblock to blockvectory
     325         [ +  + ]:        5636 :    if(iblock /=1) then
     326                 :             :      call wfcopy('D',bblocksize*vectsize,cg,1,blockvectory,1,bblocksize,iblock,'C',withbbloc=.false.,&
     327                 :        3386 : &     timopt=timopt,tim_wfcopy=tim_wfcopy)
     328                 :             : 
     329         [ +  + ]:        3386 :      if(gen_eigenpb) then
     330                 :             :        call wfcopy('D',bblocksize*vectsize,gsc,1,blockvectorby,1,bblocksize,iblock,'S',withbbloc=.false.,&
     331                 :        3258 : &       timopt=timopt,tim_wfcopy=tim_wfcopy)
     332                 :             :      else
     333   [ +  -  +  -  :       62824 :        blockvectorby = blockvectory
          -  +  -  -  -  
          -  -  -  +  +  
                   +  + ]
     334                 :             :      end if
     335                 :             : 
     336                 :             : !    b-orthogonalize x to the constraint y (supposed b-orthonormal)
     337                 :             : !    blockvectorx=blockvectorx-matmul(blockvectory,matmul((blockvectorby)^T,blockvectorx))
     338                 :             : 
     339                 :             :      call abi_xgemm(cparam(cplx),'n',bblocksize,blocksize,vectsize,cone,blockvectorby,&
     340                 :        3386 : &     vectsize,blockvectorx,vectsize,czero,gramyx,bblocksize,x_cplx=x_cplx)
     341                 :             : 
     342         [ -  + ]:        3386 :      if(abs(dtset%timopt)==3) then
     343                 :           0 :        call timab(533,1,tsec)
     344                 :             :      end if
     345                 :        3386 :      call xmpi_sum(gramyx,mpi_enreg%comm_bandspinorfft,ierr)
     346         [ -  + ]:        3386 :      if(abs(dtset%timopt)==3) then
     347                 :           0 :        call timab(533,2,tsec)
     348                 :             :      end if
     349                 :             : 
     350                 :             :      call abi_xgemm('n','n',vectsize,blocksize,bblocksize,cminusone,blockvectory,&
     351                 :        3386 : &     vectsize,gramyx,bblocksize,cone,blockvectorx,vectsize,x_cplx=x_cplx)
     352                 :             : 
     353                 :             :    end if
     354                 :             : !  !!!!!!!!!!!!!!!!!!!!!!!! End if iblock /=1 !!!!!!!!!!!!!!!!!!!!!!!!!!!
     355                 :             : 
     356                 :             :    call wfcopy('I',vectsize*blocksize,blockvectorx,1,cwavef,1,blocksize,iblock,'W',withbbloc=.false.,&
     357                 :        5636 : &   timopt=timopt,tim_wfcopy=tim_wfcopy)
     358                 :             : 
     359         [ -  + ]:        5636 :    if(abs(dtset%timopt)==4) then
     360                 :           0 :      call timab(521,2,tsec)
     361                 :             :    end if
     362         [ -  + ]:        5636 :    if(abs(dtset%timopt)==4) then
     363                 :           0 :      call timab(526,1,tsec)
     364                 :             :    end if
     365                 :             : 
     366         [ +  + ]:        5636 :    cpopt=-1;sij_opt=0;if (gen_eigenpb) sij_opt=1
     367                 :             : 
     368                 :             : #ifdef HAVE_OPENMP_OFFLOAD
     369                 :             :    !$OMP TARGET UPDATE TO(cwavef) IF(dtset%gpu_option==ABI_GPU_OPENMP)
     370                 :             : #endif
     371         [ +  + ]:        5636 :    if (mpi_enreg%paral_kgb==0) then
     372                 :             :      call getghc(cpopt,cwavef,cprj_dum,gwavef,swavef,gs_hamk,gvnlxc,dum,&
     373                 :         836 : &     mpi_enreg,blocksize,prtvol,sij_opt,tim_getghc,0)
     374                 :             :    else
     375                 :             :      call prep_getghc(cwavef,gs_hamk,gvnlxc,gwavef,swavef,dum,blocksize,mpi_enreg,&
     376                 :        4800 : &     prtvol,sij_opt,cpopt,cprj_dum,already_transposed=.false.)
     377                 :             :    end if
     378                 :             : #ifdef HAVE_OPENMP_OFFLOAD
     379                 :             :    !$OMP TARGET UPDATE FROM(gwavef,gvnlxc,swavef) IF(dtset%gpu_option==ABI_GPU_OPENMP)
     380                 :             : #endif
     381         [ -  + ]:        5636 :    if(abs(dtset%timopt)==4) then
     382                 :           0 :      call timab(526,2,tsec)
     383                 :             :    end if
     384         [ -  + ]:        5636 :    if(abs(dtset%timopt)==4) then
     385                 :           0 :      call timab(522,1,tsec)
     386                 :             :    end if
     387                 :             : 
     388         [ +  + ]:        5636 :    if ( gen_eigenpb ) then
     389                 :             :      call wfcopy('D',vectsize*blocksize,swavef,1,blockvectorbx,1,blocksize,iblock,'W',withbbloc=.false.,&
     390                 :        4804 : &     timopt=timopt,tim_wfcopy=tim_wfcopy)
     391                 :             :    else
     392                 :             :      call wfcopy('D',vectsize*blocksize,gvnlxc,1,blockvectorvx,1,blocksize,iblock,'W',withbbloc=.false.,&
     393                 :         832 : &     timopt=timopt,tim_wfcopy=tim_wfcopy)
     394   [ +  -  +  -  :      626448 :      blockvectorbx = blockvectorx
          -  +  -  -  -  
          -  -  -  +  +  
                   +  + ]
     395                 :             :    end if
     396                 :             : 
     397                 :             :    call wfcopy('D',vectsize*blocksize,gwavef,1,blockvectorax,1,blocksize,iblock,'W',withbbloc=.false.,&
     398                 :        5636 : &   timopt=timopt,tim_wfcopy=tim_wfcopy)
     399                 :             : 
     400                 :             :    call abi_xorthonormalize(blockvectorx,blockvectorbx,blocksize,mpi_enreg%comm_bandspinorfft,gramxbx,vectsize,&
     401                 :        5636 : &   x_cplx,timopt=timopt,tim_xortho=tim_xortho)
     402                 :             : 
     403                 :        5636 :    call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,gramxbx,blocksize,blockvectorbx,vectsize,x_cplx=x_cplx)
     404                 :        5636 :    call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,gramxbx,blocksize,blockvectorax,vectsize,x_cplx=x_cplx)
     405                 :             : 
     406         [ +  + ]:        5636 :    if (gs_hamk%usepaw==0) then
     407                 :         832 :      call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,gramxbx,blocksize,blockvectorvx,vectsize,x_cplx=x_cplx)
     408                 :             :    end if
     409                 :             : 
     410                 :             : !  Do rayleigh ritz on a in space x
     411                 :             : !  gramxax=matmul(transpose(blockvectorx),blockvectorax)
     412                 :             :    call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorx,&
     413                 :        5636 : &   vectsize,blockvectorax,vectsize,czero,gramxax,blocksize,x_cplx=x_cplx)
     414                 :             : 
     415         [ -  + ]:        5636 :    if(abs(dtset%timopt)==3) then
     416                 :           0 :      call timab(533,1,tsec)
     417                 :             :    end if
     418                 :        5636 :    call xmpi_sum(gramxax,mpi_enreg%comm_bandspinorfft,ierr)
     419         [ -  + ]:        5636 :    if(abs(dtset%timopt)==3) then
     420                 :           0 :      call timab(533,2,tsec)
     421                 :             :    end if
     422   [ -  +  +  -  :       11272 :    ABI_MALLOC(eigen,(blocksize))
             -  +  -  + ]
     423                 :             : 
     424                 :             :    call abi_xheev('v','u',blocksize,gramxax,blocksize,eigen,x_cplx=cplx,istwf_k=istwf_k, &
     425                 :        5636 :    timopt=timopt,tim_xeigen=tim_xeigen,use_slk=dtset%use_slk,use_gpu_magma=use_lapack_gpu)
     426                 :             : 
     427                 :             : !  blockvectorx=matmul(blockvectorx,gramxax)
     428                 :             :    call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorx,&
     429                 :        5636 : &   vectsize,gramxax,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
     430   [ +  -  +  -  :     6334194 :    blockvectorx = blockvectordumm
          -  +  -  -  +  
                +  +  + ]
     431                 :             : 
     432                 :             : !  blockvectorax=matmul(blockvectorax,gramxax)
     433                 :             :    call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorax,&
     434                 :        5636 : &   vectsize,gramxax,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
     435   [ +  -  +  -  :     6334194 :    blockvectorax = blockvectordumm
          -  +  -  -  +  
                +  +  + ]
     436                 :             : 
     437                 :             : !  blockvectorvx=matmul(blockvectorvx,gramxax)
     438         [ +  + ]:        5636 :    if (gs_hamk%usepaw==0) then
     439                 :             :      call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorvx,&
     440                 :         832 : &     vectsize,gramxax,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
     441   [ +  -  +  -  :      626448 :      blockvectorvx = blockvectordumm
          -  +  -  -  +  
                +  +  + ]
     442                 :             :    end if
     443                 :             : 
     444                 :             : !  blockvectorbx=matmul(blockvectorbx,gramxax)
     445                 :             :    call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorbx,&
     446                 :        5636 : &   vectsize,gramxax,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
     447   [ +  -  +  -  :     6334194 :    blockvectorbx = blockvectordumm
          -  +  -  -  +  
                +  +  + ]
     448                 :             : 
     449                 :             : #if FC_CRAY
     450                 :             :    lambda(:,:) = zero
     451                 :             :    do iblocksize=1,blocksize
     452                 :             :      lambda(cplx*(iblocksize-1)+1,iblocksize) = eigen(iblocksize)
     453                 :             :    end do
     454                 :             : #else
     455         [ +  + ]:       32336 :    do iblocksize=1,blocksize
     456         [ +  + ]:       80100 :      zvar=(/eigen(iblocksize),zero/)
     457                 :       32336 :      call abi_xcopy(1,zvar,1,lambda(cplx*(iblocksize-1)+1:cplx*iblocksize,iblocksize),1,x_cplx=x_cplx)
     458                 :             :    end do
     459                 :             : #endif
     460                 :             : 
     461         [ -  + ]:        5636 :    ABI_FREE(eigen)
     462                 :             : 
     463         [ -  + ]:        5636 :    if(abs(dtset%timopt)==4) then
     464                 :           0 :      call timab(522,2,tsec)
     465                 :             :    end if
     466                 :             : 
     467                 :             : !  ###########################################################################
     468                 :             : !  ################ PERFORM LOOP ON NLINE ####################################
     469                 :             : !  ###########################################################################
     470                 :             : !  now the main alogrithm
     471         [ +  + ]:       33474 :    iter: do iterationnumber=1,maxiterations
     472                 :             : 
     473         [ -  + ]:       28923 :      if(abs(dtset%timopt)==4) then
     474                 :           0 :        call timab(523,1,tsec)
     475                 :             :      end if
     476                 :             : 
     477                 :             : !    Build residual
     478                 :             : !    blockvectorr=blockvectorax-matmul(blockvectorx,lambda)
     479         [ +  - ]:       28923 :      if (gs_hamk%use_gbt == 0) then
     480                 :             :        call xprecon(blockvectorbx,lambda,blocksize,&
     481                 :             : &       iterationnumber,kinpw,mpi_enreg,npw_k,my_nspinor,&
     482                 :       28923 : &       optekin,optpcon,pcon,blockvectorax,blockvectorr,vectsize,timopt=timopt,tim_xprecon=tim_xprecon)
     483                 :             :      else
     484                 :             :        call xprecon(blockvectorbx,lambda,blocksize,&
     485                 :             : &       iterationnumber,kinpw,mpi_enreg,npw_k,my_nspinor,&
     486                 :             : &       optekin,optpcon,pcon,blockvectorax,blockvectorr,vectsize,timopt=timopt,tim_xprecon=tim_xprecon,&
     487   [ #  #  #  #  :           0 : &       kinpw_kp=kinpw_kp,pcon_kp=pcon_kp)
             #  #  #  # ]
     488                 :             :      end if
     489                 :             : 
     490   [ +  -  -  +  :    33624786 :      residualnorms=sum(blockvectorr**2,dim=1)
          -  -  +  +  +  
                      + ]
     491                 :             : 
     492         [ -  + ]:       28923 :      if(abs(dtset%timopt)==3) then
     493                 :           0 :        call timab(533,1,tsec)
     494                 :             :      end if
     495                 :       28923 :      call xmpi_sum(residualnorms,mpi_enreg%comm_bandspinorfft,ierr)
     496         [ -  + ]:       28923 :      if(abs(dtset%timopt)==3) then
     497                 :           0 :        call timab(533,2,tsec)
     498                 :             :      end if
     499                 :             : 
     500         [ +  + ]:      192104 :      resid_k(bblocksize+1:bblocksize+blocksize)=residualnorms(1:blocksize)
     501                 :             : 
     502                 :             : !    If residual sufficiently small stop line minimizations
     503   [ +  -  -  +  :      221027 :      if (abs(maxval(residualnorms(1:blocksize)))<dtset%tolwfr_diago) then
          -  -  +  +  +  
                      + ]
     504         [ -  + ]:           8 :        if (prtvol > 0) then
     505                 :             :          write(message, '(a,i0,a,i0,a,es12.4)' ) &
     506                 :           0 : &         ' lobpcgwf: block ',iblock,' converged after ',iterationnumber,&
     507   [ #  #  #  #  :           0 : &         ' line minimizations: maxval(resid(1:blocksize)) =',maxval(residualnorms(1:blocksize))
             #  #  #  # ]
     508                 :           0 :          call wrtout(std_out,message,'PERS')
     509                 :             :        end if
     510                 :             :        havetoprecon=.false.
     511                 :             :        exit
     512                 :             :      end if
     513                 :             : 
     514         [ -  + ]:       28915 :      if(use_linalg_gpu==1) then
     515                 :           0 :        call copy_on_gpu(blockvectorr, blockvectorr_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     516                 :             :      end if
     517                 :             : 
     518         [ +  + ]:       28915 :      if(iblock /=1) then
     519                 :             : !      Residuals orthogonal to blockvectorby
     520                 :             : !      blockvectorr=blockvectorr-matmul(blockvectory,matmul((blockvectorby)^T,blockvectorr))
     521                 :             : 
     522         [ -  + ]:       13700 :        if(use_linalg_gpu==1) then
     523                 :           0 :          call copy_on_gpu(blockvectorby, bblockvector_gpu, INT(cplx, c_size_t)*dp*vectsize*bblocksize)
     524                 :             :          call gpu_xgemm(cplx,cparam(cplx),'n',bblocksize,blocksize,vectsize,cone,bblockvector_gpu,&
     525                 :           0 : &         vectsize,blockvectorr_gpu,vectsize,czero,gram_gpu,bblocksize)
     526                 :           0 :          call copy_from_gpu(gramyx, gram_gpu, INT(cplx, c_size_t)*dp*bblocksize*blocksize)
     527                 :             :        else
     528                 :             :          call abi_xgemm(cparam(cplx),'n',bblocksize,blocksize,vectsize,cone,blockvectorby,&
     529                 :       13700 : &         vectsize,blockvectorr,vectsize,czero,gramyx,bblocksize,x_cplx=x_cplx)
     530                 :             :        end if
     531                 :             : 
     532         [ -  + ]:       13700 :        if(abs(dtset%timopt)==3) then
     533                 :           0 :          call timab(533,1,tsec)
     534                 :             :        end if
     535                 :       13700 :        call xmpi_sum(gramyx,mpi_enreg%comm_bandspinorfft,ierr)
     536         [ -  + ]:       13700 :        if(abs(dtset%timopt)==3) then
     537                 :           0 :          call timab(533,2,tsec)
     538                 :             :        end if
     539                 :             : 
     540         [ -  + ]:       13700 :        if(use_linalg_gpu==1) then
     541                 :           0 :          call copy_on_gpu(gramyx,       gram_gpu,         INT(cplx, c_size_t)*dp*bblocksize*blocksize)
     542                 :           0 :          call copy_on_gpu(blockvectory, bblockvector_gpu, INT(cplx, c_size_t)*dp*vectsize*bblocksize)
     543                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,bblocksize,cminusone,bblockvector_gpu,&
     544                 :           0 : &         vectsize,gram_gpu,bblocksize,cone,blockvectorr_gpu,vectsize)
     545                 :             :        else
     546                 :             :          call abi_xgemm('n','n',vectsize,blocksize,bblocksize,cminusone,blockvectory,&
     547                 :       13700 : &         vectsize,gramyx,bblocksize,cone,blockvectorr,vectsize,x_cplx=x_cplx)
     548                 :             :        end if
     549                 :             : 
     550                 :             :      end if
     551                 :             : 
     552                 :             : !    Residuals orthogonal to blockvectorx
     553                 :             : !    blockvectorr=blockvectorr-matmul(blockvectorx,matmul((blockvectorbx)^T,blockvectorr))
     554         [ -  + ]:       28915 :      if(use_linalg_gpu==1) then
     555                 :           0 :        call copy_on_gpu(blockvectorbx, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     556                 :             :        call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,cone,C_gpu,&
     557                 :           0 : &       vectsize,blockvectorr_gpu,vectsize,czero,gram_gpu,blocksize)
     558                 :           0 :        call copy_from_gpu(gramxax, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     559                 :             :      else
     560                 :             :        call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorbx,&
     561                 :       28915 : &       vectsize,blockvectorr,vectsize,czero,gramxax,blocksize,x_cplx=x_cplx)
     562                 :             :      end if
     563                 :             : 
     564         [ -  + ]:       28915 :      if(abs(dtset%timopt)==3) then
     565                 :           0 :        call timab(533,1,tsec)
     566                 :             :      end if
     567                 :       28915 :      call xmpi_sum(gramxax,mpi_enreg%comm_bandspinorfft,ierr)
     568         [ -  + ]:       28915 :      if(abs(dtset%timopt)==3) then
     569                 :           0 :        call timab(533,2,tsec)
     570                 :             :      end if
     571                 :             : 
     572         [ -  + ]:       28915 :      if(use_linalg_gpu==1) then
     573                 :           0 :        call copy_on_gpu(gramxax,      gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     574                 :           0 :        call copy_on_gpu(blockvectorx, C_gpu,    INT(cplx, c_size_t)*dp*vectsize*blocksize)
     575                 :             :        call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cminusone,C_gpu,&
     576                 :           0 : &       vectsize,gram_gpu,blocksize,cone,blockvectorr_gpu,vectsize)
     577                 :           0 :        call copy_from_gpu(blockvectorr, blockvectorr_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     578                 :             :      else
     579                 :             :        call abi_xgemm('n','n',vectsize,blocksize,blocksize,cminusone,blockvectorx,&
     580                 :       28915 : &       vectsize,gramxax,blocksize,cone,blockvectorr,vectsize,x_cplx=x_cplx)
     581                 :             :      end if
     582                 :             : 
     583                 :             :      call wfcopy('I',vectsize*blocksize,blockvectorr,1,cwavef,1,blocksize,iblock,'W',withbbloc=.false.,&
     584                 :       28915 : &     timopt=timopt,tim_wfcopy=tim_wfcopy)
     585                 :             : 
     586         [ +  + ]:       28915 :      cpopt=-1;sij_opt=0;if (gen_eigenpb) sij_opt=1
     587                 :             : 
     588         [ -  + ]:       28915 :      if(abs(dtset%timopt)==4) then
     589                 :           0 :        call timab(523,2,tsec)
     590                 :             :      end if
     591         [ -  + ]:       28915 :      if(abs(dtset%timopt)==4) then
     592                 :           0 :        call timab(526,1,tsec)
     593                 :             :      end if
     594                 :             : 
     595                 :             : #ifdef HAVE_OPENMP_OFFLOAD
     596                 :             :      !$OMP TARGET UPDATE TO(cwavef) IF(dtset%gpu_option==ABI_GPU_OPENMP)
     597                 :             : #endif
     598         [ +  + ]:       28915 :      if (mpi_enreg%paral_kgb==0) then
     599                 :             :        call getghc(cpopt,cwavef,cprj_dum,gwavef,swavef,gs_hamk,gvnlxc,dum,&
     600                 :        5363 : &       mpi_enreg,blocksize,prtvol,sij_opt,tim_getghc,0)
     601                 :             :      else
     602                 :             :        call prep_getghc(cwavef,gs_hamk,gvnlxc,gwavef,swavef,dum,blocksize,mpi_enreg,&
     603                 :       23552 : &       prtvol,sij_opt,cpopt,cprj_dum,already_transposed=.false.)
     604                 :             :      end if
     605                 :             : #ifdef HAVE_OPENMP_OFFLOAD
     606                 :             :      !$OMP TARGET UPDATE FROM(gwavef,gvnlxc,swavef) IF(dtset%gpu_option==ABI_GPU_OPENMP)
     607                 :             : #endif
     608                 :             : 
     609         [ -  + ]:       28915 :      if(abs(dtset%timopt)==4) then
     610                 :           0 :        call timab(526,2,tsec)
     611                 :             :      end if
     612         [ -  + ]:       28915 :      if(abs(dtset%timopt)==4) then
     613                 :           0 :        call timab(524,1,tsec)
     614                 :             :      end if
     615                 :             : 
     616         [ +  + ]:       28915 :      if (gen_eigenpb) then
     617                 :             :        call wfcopy('D',vectsize*blocksize,swavef,1,blockvectorbr,1,blocksize,iblock,'W',withbbloc=.false.,&
     618                 :       23365 : &       timopt=timopt,tim_wfcopy=tim_wfcopy)
     619                 :             :      else
     620   [ +  -  +  -  :     4113243 :        blockvectorbr = blockvectorr
          -  +  -  -  -  
          -  -  -  +  +  
                   +  + ]
     621                 :             :        call wfcopy('D',vectsize*blocksize,gvnlxc,1,blockvectorvr,1,blocksize,iblock,'W',withbbloc=.false.,&
     622                 :        5550 : &       timopt=timopt,tim_wfcopy=tim_wfcopy)
     623                 :             :      end if
     624                 :             : 
     625                 :             :      call wfcopy('D',vectsize*blocksize,gwavef,1,blockvectorar,1,blocksize,iblock,'W',withbbloc=.false.,&
     626                 :       28915 : &     timopt=timopt,tim_wfcopy=tim_wfcopy)
     627                 :             : 
     628         [ -  + ]:       28915 :      if(use_linalg_gpu==1) then
     629                 :           0 :        call copy_on_gpu(blockvectorbr, blockvectorbr_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     630                 :             :        call gpu_xorthonormalize(blockvectorr_gpu,blockvectorbr_gpu,blocksize,mpi_enreg%comm_bandspinorfft,gram_gpu,vectsize,&
     631                 :           0 : &       x_cplx,timopt=timopt,tim_xortho=tim_xortho)
     632                 :           0 :        call copy_from_gpu(blockvectorr,  blockvectorr_gpu,  INT(cplx, c_size_t)*dp*vectsize*blocksize)
     633                 :           0 :        call gpu_xtrsm(cplx,'r','u','n','n',vectsize,blocksize,cone,gram_gpu,blocksize,blockvectorbr_gpu,vectsize)
     634                 :           0 :        call copy_from_gpu(blockvectorbr, blockvectorbr_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     635                 :           0 :        call copy_on_gpu(blockvectorar,   blockvectorar_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     636                 :           0 :        call gpu_xtrsm(cplx,'r','u','n','n',vectsize,blocksize,cone,gram_gpu,blocksize,blockvectorar_gpu,vectsize)
     637                 :           0 :        call copy_from_gpu(blockvectorar, blockvectorar_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     638         [ #  # ]:           0 :        if (gs_hamk%usepaw==0) then
     639                 :           0 :          call copy_on_gpu(blockvectorvr,   A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     640                 :           0 :          call gpu_xtrsm(cplx,'r','u','n','n',vectsize,blocksize,cone,gram_gpu,blocksize,A_gpu,vectsize)
     641                 :           0 :          call copy_from_gpu(blockvectorvr, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     642                 :             :        end if
     643                 :           0 :        call copy_from_gpu(gramrbr, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     644                 :             :      else
     645                 :             :        call abi_xorthonormalize(blockvectorr,blockvectorbr,blocksize,mpi_enreg%comm_bandspinorfft,gramrbr,vectsize,&
     646                 :       28915 : &       x_cplx,timopt=timopt,tim_xortho=tim_xortho)
     647                 :       28915 :        call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,gramrbr,blocksize,blockvectorbr,vectsize,x_cplx=x_cplx)
     648                 :       28915 :        call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,gramrbr,blocksize,blockvectorar,vectsize,x_cplx=x_cplx)
     649         [ +  + ]:       28915 :        if (gs_hamk%usepaw==0) then
     650                 :        5550 :          call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,gramrbr,blocksize,blockvectorvr,vectsize,x_cplx=x_cplx)
     651                 :             :        end if
     652                 :             :      end if
     653                 :             : 
     654         [ +  + ]:       28915 :      if(iterationnumber>1) then
     655         [ -  + ]:       23283 :        if(use_linalg_gpu==1) then
     656                 :           0 :          call copy_on_gpu(blockvectorp,    A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     657                 :           0 :          call copy_on_gpu(blockvectorbp,   C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     658                 :             :          call gpu_xorthonormalize(A_gpu,C_gpu,blocksize,mpi_enreg%comm_bandspinorfft,gram_gpu,vectsize,&
     659                 :           0 : &         x_cplx,timopt=timopt,tim_xortho=tim_xortho)
     660                 :           0 :          call copy_from_gpu(blockvectorp,  A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     661                 :           0 :          call gpu_xtrsm(cplx,'r','u','n','n',vectsize,blocksize,cone,gram_gpu,blocksize,C_gpu,vectsize)
     662                 :           0 :          call copy_from_gpu(blockvectorbp, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     663                 :           0 :          call copy_on_gpu(blockvectorap,   A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     664                 :           0 :          call gpu_xtrsm(cplx,'r','u','n','n',vectsize,blocksize,cone,gram_gpu,blocksize,A_gpu,vectsize)
     665                 :           0 :          call copy_from_gpu(blockvectorap, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     666         [ #  # ]:           0 :          if (gs_hamk%usepaw==0) then
     667                 :           0 :            call copy_on_gpu(blockvectorvp, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     668                 :           0 :            call gpu_xtrsm(cplx,'r','u','n','n',vectsize,blocksize,cone,gram_gpu,blocksize,A_gpu,vectsize)
     669                 :           0 :            call copy_from_gpu(blockvectorvp, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     670                 :             :          end if
     671                 :           0 :          call copy_from_gpu(grampbp, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     672                 :             :        else
     673                 :             : !        call orthonormalize(blockvectorp,blockvectorbp,blockvectorap)
     674                 :             :          call abi_xorthonormalize(blockvectorp,blockvectorbp,blocksize,mpi_enreg%comm_bandspinorfft,grampbp,vectsize,&
     675                 :       23283 : &         x_cplx,timopt=timopt,tim_xortho=tim_xortho)
     676                 :             : !        blockvectorap=matmul(blockvectorap,grampbp)
     677                 :       23283 :          call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,grampbp,blocksize,blockvectorbp,vectsize,x_cplx=x_cplx)
     678                 :       23283 :          call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,grampbp,blocksize,blockvectorap,vectsize,x_cplx=x_cplx)
     679         [ +  + ]:       23283 :          if (gs_hamk%usepaw==0) then
     680                 :        4722 :            call abi_xtrsm('r','u','n','n',vectsize,blocksize,cone,grampbp,blocksize,blockvectorvp,vectsize,x_cplx=x_cplx)
     681                 :             :          end if
     682                 :             :        end if
     683                 :             :      end if
     684                 :             : 
     685                 :       28915 :      activersize=blocksize
     686         [ +  + ]:       28915 :      if (iterationnumber==1) then
     687                 :       28915 :        activepsize=0
     688                 :             :        restart=1
     689                 :             :      else
     690                 :       23283 :        activepsize=blocksize
     691                 :       23283 :        restart=0
     692                 :             :      end if
     693                 :             : 
     694                 :             : !    gramxar=matmul((blockvectorax)^T,blockvectorr)
     695                 :             : !    gramrar=matmul((blockvectorar)^T,blockvectorr)
     696                 :             : !    gramxax=matmul((blockvectorax)^T,blockvectorx)
     697         [ -  + ]:       28915 :      if(use_linalg_gpu==1) then
     698                 :           0 :        call copy_on_gpu(blockvectorax, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     699                 :             : 
     700                 :             :        call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,cone,A_gpu,&
     701                 :           0 : &       vectsize,blockvectorr_gpu,vectsize,czero,gram_gpu,blocksize)
     702                 :           0 :        call copy_from_gpu(gramxar, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     703                 :             : 
     704                 :             :        call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorar_gpu,&
     705                 :           0 : &       vectsize,blockvectorr_gpu,vectsize,czero,gram_gpu,blocksize)
     706                 :           0 :        call copy_from_gpu(gramrar, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     707                 :             : 
     708                 :           0 :        call copy_on_gpu(blockvectorx, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     709                 :             :        call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,cone,A_gpu,&
     710                 :           0 : &       vectsize,C_gpu,vectsize,czero,gram_gpu,blocksize)
     711                 :             : 
     712                 :           0 :        call copy_from_gpu(gramxax, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     713                 :             :      else
     714                 :             :        call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorax,&
     715                 :       28915 : &       vectsize,blockvectorr,vectsize,czero,gramxar,blocksize,x_cplx=x_cplx)
     716                 :             :        call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorar,&
     717                 :       28915 : &       vectsize,blockvectorr,vectsize,czero,gramrar,blocksize,x_cplx=x_cplx)
     718                 :             :        call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorax,&
     719                 :       28915 : &       vectsize,blockvectorx,vectsize,czero,gramxax,blocksize,x_cplx=x_cplx)
     720                 :             :      end if
     721                 :             : 
     722   [ +  +  +  + ]:     2833751 :      transf3(:,:,1) = gramxar
     723   [ +  +  +  + ]:     2833751 :      transf3(:,:,2) = gramrar
     724   [ +  +  +  + ]:     2833751 :      transf3(:,:,3) = gramxax
     725         [ -  + ]:       28915 :      if(abs(dtset%timopt)==3) then
     726                 :           0 :        call timab(533,1,tsec)
     727                 :             :      end if
     728                 :       28915 :      call xmpi_sum(transf3,mpi_enreg%comm_bandspinorfft,ierr)
     729         [ -  + ]:       28915 :      if(abs(dtset%timopt)==3) then
     730                 :           0 :        call timab(533,2,tsec)
     731                 :             :      end if
     732                 :             : 
     733   [ +  -  +  -  :     2862666 :      gramxar = transf3(:,:,1)
          -  +  -  -  +  
                +  +  + ]
     734   [ +  -  +  -  :     2862666 :      gramrar = transf3(:,:,2)
          -  +  -  -  +  
                +  +  + ]
     735   [ +  -  +  -  :     2862666 :      gramxax = transf3(:,:,3)
          -  +  -  -  +  
                +  +  + ]
     736                 :             : 
     737                 :             : !    gramxbx=matmul((blockvectorbx)^T,blockvectorx)
     738                 :             : !    gramrbr=matmul((blockvectorbr)^T,blockvectorr)
     739                 :             : !    gramxbr=matmul((blockvectorbx)^T,blockvectorr)
     740                 :             : !    Note that the gramb matrix is more easier to construct than grama:
     741                 :             : !    i) <x|B|x>=<r|B|r>=<p|B|p>=(1;0)
     742                 :             : !    since the x, r and p blockvector are normalized
     743                 :             : !    ii) <r|B|x>=(0;0)
     744                 :             : !    since the x and r blockvector are orthogonalized
     745                 :             : !    iii) The <p|B|r> and <p|B|x> have to be computed.
     746   [ +  +  +  + ]:     2833751 :      gramxbx(:,:)=zero
     747   [ +  +  +  + ]:     2833751 :      gramrbr(:,:)=zero
     748   [ +  +  +  + ]:     2833751 :      gramxbr(:,:)=zero
     749         [ +  + ]:      192072 :      do iblocksize=1,blocksize
     750                 :      163157 :        gramxbx(cplx*(iblocksize-1)+1,iblocksize) = one
     751                 :      192072 :        gramrbr(cplx*(iblocksize-1)+1,iblocksize) = one
     752                 :             :      end do
     753                 :             : 
     754                 :             : !    ###########################################################################
     755                 :             : !    ################ PERFORM LOOP ON COND #####################################
     756                 :             : !    ###########################################################################
     757                 :             : 
     758                 :       28915 :      i1=0;i2=blocksize;i3=2*blocksize;i4=3*blocksize
     759         [ +  - ]:       28915 :      cond: do cond_try=1,2 !2 when restart
     760         [ +  + ]:       28915 :        if (restart==0) then
     761                 :             : 
     762                 :             : !        gramxap=matmul((blockvectorax)^T,blockvectorp)
     763                 :             : !        gramrap=matmul((blockvectorar)^T,blockvectorp)
     764                 :             : !        grampap=matmul((blockvectorap)^T,blockvectorp)
     765                 :             : !        gramxbp=matmul((blockvectorbx)^T,blockvectorp)
     766                 :             : !        gramrbp=matmul((blockvectorbr)^T,blockvectorp)
     767                 :             : !        grampbp=matmul((blockvectorbp)^T,blockvectorp)
     768         [ -  + ]:       23283 :          if(use_linalg_gpu==1) then
     769                 :           0 :            call copy_on_gpu(blockvectorp, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     770                 :           0 :            call copy_on_gpu(blockvectorax,A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     771                 :             :            call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,&
     772                 :           0 : &           cone,A_gpu,vectsize,C_gpu,vectsize,czero,gram_gpu,blocksize)
     773                 :           0 :            call copy_from_gpu(gramxap, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     774                 :             :            call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,&
     775                 :           0 : &           cone,blockvectorar_gpu,vectsize,C_gpu,vectsize,czero,gram_gpu,blocksize)
     776                 :           0 :            call copy_from_gpu(gramrap, gram_gpu,  INT(cplx, c_size_t)*dp*blocksize*blocksize)
     777                 :           0 :            call copy_on_gpu(blockvectorap, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     778                 :             :            call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,&
     779                 :           0 : &           cone,A_gpu,vectsize,C_gpu,vectsize,czero,gram_gpu,blocksize)
     780                 :           0 :            call copy_from_gpu(grampap, gram_gpu,  INT(cplx, c_size_t)*dp*blocksize*blocksize)
     781                 :           0 :            call copy_on_gpu(blockvectorbx, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     782                 :             :            call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,&
     783                 :           0 : &           cone,A_gpu,vectsize,C_gpu,vectsize,czero,gram_gpu,blocksize)
     784                 :           0 :            call copy_from_gpu(gramxbp, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     785                 :             :            call gpu_xgemm(cplx,cparam(cplx),'n',blocksize,blocksize,vectsize,&
     786                 :           0 : &           cone,blockvectorbr_gpu,vectsize,C_gpu,vectsize,czero,gram_gpu,blocksize)
     787                 :           0 :            call copy_from_gpu(gramrbp, gram_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     788                 :             :          else
     789                 :             :            call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorax,&
     790                 :       23283 : &           vectsize,blockvectorp,vectsize,czero,gramxap,blocksize,x_cplx=x_cplx)
     791                 :             :            call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorar,&
     792                 :       23283 : &           vectsize,blockvectorp,vectsize,czero,gramrap,blocksize,x_cplx=x_cplx)
     793                 :             :            call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorap,&
     794                 :       23283 : &           vectsize,blockvectorp,vectsize,czero,grampap,blocksize,x_cplx=x_cplx)
     795                 :             :            call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorbx,&
     796                 :       23283 : &           vectsize,blockvectorp,vectsize,czero,gramxbp,blocksize,x_cplx=x_cplx)
     797                 :             :            call abi_xgemm(cparam(cplx),'n',blocksize,blocksize,vectsize,cone,blockvectorbr,&
     798                 :       23283 : &           vectsize,blockvectorp,vectsize,czero,gramrbp,blocksize,x_cplx=x_cplx)
     799                 :             :          end if
     800                 :             : !        It's not necessary to compute the last one: <p|B|p>=(1;0) (see above)
     801   [ +  +  +  + ]:     2402221 :          transf5(:,:,1)=gramxap(:,:)
     802   [ +  +  +  + ]:     2402221 :          transf5(:,:,2)=gramrap(:,:)
     803   [ +  +  +  + ]:     2402221 :          transf5(:,:,3)=grampap(:,:)
     804   [ +  +  +  + ]:     2402221 :          transf5(:,:,4)=gramxbp(:,:)
     805   [ +  +  +  + ]:     2402221 :          transf5(:,:,5)=gramrbp(:,:)
     806         [ -  + ]:       23283 :          if(abs(dtset%timopt)==3) then
     807                 :           0 :            call timab(533,1,tsec)
     808                 :             :          end if
     809                 :       23283 :          call xmpi_sum(transf5,mpi_enreg%comm_bandspinorfft,ierr)
     810         [ -  + ]:       23283 :          if(abs(dtset%timopt)==3) then
     811                 :           0 :            call timab(533,2,tsec)
     812                 :             :          end if
     813   [ +  +  +  + ]:     2402221 :          gramxap(:,:)=transf5(:,:,1)
     814   [ +  +  +  + ]:     2402221 :          gramrap(:,:)=transf5(:,:,2)
     815   [ +  +  +  + ]:     2402221 :          grampap(:,:)=transf5(:,:,3)
     816   [ +  +  +  + ]:     2402221 :          gramxbp(:,:)=transf5(:,:,4)
     817   [ +  +  +  + ]:     2402221 :          gramrbp(:,:)=transf5(:,:,5)
     818   [ +  +  +  + ]:     2402221 :          grampbp(:,:)=zero
     819         [ +  + ]:      159752 :          do iblocksize=1,blocksize
     820                 :      159752 :            grampbp(cplx*(iblocksize-1)+1,iblocksize) = one
     821                 :             :          end do
     822                 :       23283 :          bigorder=i4
     823   [ -  +  +  -  :       93132 :          ABI_MALLOC(grama,(cplx*i4,i4))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     824   [ -  +  +  -  :       69849 :          ABI_MALLOC(gramb,(cplx*i4,i4))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     825   [ -  +  +  -  :       69849 :          ABI_MALLOC(eigen,(i4))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     826                 :             : !        ABI_MALLOC(coordx,(cplx*i4,blocksize))
     827   [ -  +  +  -  :       93132 :          ABI_MALLOC(coordx1,(cplx*blocksize,blocksize))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     828   [ -  +  +  -  :       69849 :          ABI_MALLOC(coordx2,(cplx*blocksize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     829   [ -  +  +  -  :       69849 :          ABI_MALLOC(coordx3,(cplx*blocksize,blocksize))
          -  +  +  -  -  
             +  -  +  -  
                      + ]
     830   [ +  +  +  +  :    41206539 :          grama(:,:)=zero;gramb(:,:)=zero
             +  +  +  + ]
     831   [ +  +  +  + ]:     2402221 :          grama(gramindex(i1+1):gramindex(i2)+cplx-1,i1+1:i2)=gramxax
     832   [ +  +  +  + ]:     2402221 :          grama(gramindex(i1+1):gramindex(i2)+cplx-1,i2+1:i3)=gramxar
     833   [ +  +  +  + ]:     2402221 :          grama(gramindex(i1+1):gramindex(i2)+cplx-1,i3+1:i4)=gramxap
     834         [ +  + ]:     2378938 :          grama(gramindex(i2+1):gramindex(i3)+cplx-1,i2+1:i3)=gramrar
     835   [ +  +  +  + ]:     2402221 :          grama(gramindex(i2+1):gramindex(i3)+cplx-1,i3+1:i4)=gramrap
     836   [ +  +  +  + ]:     2402221 :          grama(gramindex(i3+1):gramindex(i4)+cplx-1,i3+1:i4)=grampap
     837   [ +  +  +  + ]:     2402221 :          gramb(gramindex(i1+1):gramindex(i2)+cplx-1,i1+1:i2)=gramxbx
     838   [ +  +  +  + ]:     2402221 :          gramb(gramindex(i1+1):gramindex(i2)+cplx-1,i2+1:i3)=gramxbr
     839   [ +  +  +  + ]:     2402221 :          gramb(gramindex(i1+1):gramindex(i2)+cplx-1,i3+1:i4)=gramxbp
     840   [ +  +  +  + ]:     2402221 :          gramb(gramindex(i2+1):gramindex(i3)+cplx-1,i2+1:i3)=gramrbr
     841   [ +  +  +  + ]:     2402221 :          gramb(gramindex(i2+1):gramindex(i3)+cplx-1,i3+1:i4)=gramrbp
     842   [ +  +  +  + ]:     2425504 :          gramb(gramindex(i3+1):gramindex(i4)+cplx-1,i3+1:i4)=grampbp
     843                 :             :        else
     844                 :        5632 :          bigorder=i3
     845   [ -  +  +  -  :       22528 :          ABI_MALLOC(grama,(cplx*i3,i3))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     846   [ -  +  +  -  :       22528 :          ABI_MALLOC(gramb,(cplx*i3,i3))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     847   [ -  +  +  -  :       16896 :          ABI_MALLOC(eigen,(i3))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     848                 :             : !        ABI_MALLOC(coordx,(cplx*i3,blocksize))
     849   [ -  +  +  -  :       22528 :          ABI_MALLOC(coordx1,(cplx*blocksize,blocksize))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     850   [ -  +  +  -  :       22528 :          ABI_MALLOC(coordx2,(cplx*blocksize,blocksize))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     851   [ +  +  +  +  :     3306064 :          grama(:,:)=zero;gramb(:,:)=zero
             +  +  +  + ]
     852   [ +  +  +  + ]:      431530 :          grama(gramindex(i1+1):gramindex(i2)+cplx-1,i1+1:i2)=gramxax
     853   [ +  +  +  + ]:      431530 :          grama(gramindex(i1+1):gramindex(i2)+cplx-1,i2+1:i3)=gramxar
     854   [ +  +  +  + ]:      431530 :          grama(gramindex(i2+1):gramindex(i3)+cplx-1,i2+1:i3)=gramrar
     855   [ +  +  +  + ]:      431530 :          gramb(gramindex(i1+1):gramindex(i2)+cplx-1,i1+1:i2)=gramxbx
     856   [ +  +  +  + ]:      431530 :          gramb(gramindex(i1+1):gramindex(i2)+cplx-1,i2+1:i3)=gramxbr
     857   [ +  +  +  + ]:      431530 :          gramb(gramindex(i2+1):gramindex(i3)+cplx-1,i2+1:i3)=gramrbr
     858                 :             :        end if
     859                 :             : 
     860   [ -  +  +  -  :      115660 :        ABI_MALLOC(tmpgramb,(cplx*bigorder,bigorder))
          -  +  +  -  +  
          -  +  -  -  +  
                   -  + ]
     861   [ -  +  +  -  :       86745 :        ABI_MALLOC(tmpeigen,(bigorder))
          +  -  +  -  -  
                +  -  + ]
     862   [ +  -  +  -  :    22299674 :        tmpgramb=gramb
          -  +  -  -  +  
                +  +  + ]
     863                 :             : 
     864                 :             :        call abi_xheev('v','u',bigorder,tmpgramb,bigorder,tmpeigen,x_cplx=cplx,istwf_k=istwf_k, &
     865                 :       28915 : &       timopt=timopt,tim_xeigen=tim_xeigen,use_slk=dtset%use_slk,use_gpu_magma=use_lapack_gpu)
     866                 :             : 
     867                 :       28915 :        condestgramb=tmpeigen(bigorder)/tmpeigen(1)
     868         [ -  + ]:       28915 :        ABI_FREE(tmpgramb)
     869                 :       28915 :        ABI_FREE(tmpeigen)
     870                 :             : 
     871   [ +  -  -  + ]:       28915 :        if (condestgramb.gt.1d+5.or.condestgramb.lt.0.d0.or.info/=0) then
     872                 :           0 :          write(std_out,*)'condition number of the Gram matrix = ',condestgramb
     873         [ #  # ]:           0 :          if (cond_try==1.and.restart==0) then
     874         [ #  # ]:           0 :            ABI_FREE(grama)
     875         [ #  # ]:           0 :            ABI_FREE(gramb)
     876         [ #  # ]:           0 :            ABI_FREE(eigen)
     877                 :             : !          ABI_FREE(coordx)
     878         [ #  # ]:           0 :            ABI_FREE(coordx1)
     879         [ #  # ]:           0 :            ABI_FREE(coordx2)
     880         [ #  # ]:           0 :            if(bigorder==i4) then
     881         [ #  # ]:           0 :              ABI_FREE(coordx3)
     882                 :             :            end if
     883         [ #  # ]:           0 :            if (nrestart.gt.1) then
     884                 :           0 :              ABI_WARNING('the minimization is stopped for this block')
     885                 :           0 :              exit iter
     886                 :             :            else
     887                 :           0 :              restart=1
     888                 :           0 :              nrestart=nrestart+1
     889                 :           0 :              call wrtout(std_out,'Lobpcgwf: restart performed',"PERS")
     890                 :             :            end if
     891                 :             :          else
     892                 :           0 :            ABI_WARNING('Gramm matrix ill-conditionned: results may be unpredictable')
     893                 :             :          end if
     894                 :             :        else
     895                 :             :          exit cond
     896                 :             :        end if
     897                 :             :      end do cond
     898                 :             : 
     899                 :             : !    ###########################################################################
     900                 :             : !    ################ END LOOP ON COND #########################################
     901                 :             : !    ###########################################################################
     902                 :             : 
     903                 :             :      call abi_xhegv(1,'v','u',bigorder,grama,bigorder,gramb,bigorder,eigen,x_cplx=cplx,istwf_k=istwf_k, &
     904                 :       28915 :      timopt=timopt,tim_xeigen=tim_xeigen,use_slk=dtset%use_slk,use_gpu_magma=use_lapack_gpu)
     905                 :             : 
     906                 :       28915 :      deltae=-one
     907         [ +  + ]:      192072 :      do iblocksize=1,blocksize
     908         [ +  + ]:      479756 :        zvar(1:cplx)=lambda(cplx*(iblocksize-1)+1:cplx*(iblocksize-1)+cplx,iblocksize)
     909                 :      163157 :        deltae=max(deltae,abs(cmplx(zvar(1),zvar(2))-eigen(iblocksize)))
     910                 :             : #ifdef FC_CRAY
     911                 :             :        ! Weird numerical error occurs with Cray when using abi_xcopy
     912                 :             :        lambda(cplx*(iblocksize-1)+1,iblocksize) = eigen(iblocksize)
     913                 :             :        if (cplx==2) lambda(cplx*iblocksize,iblocksize) = zero
     914                 :             : #else
     915         [ +  + ]:      489471 :        zvar=(/eigen(iblocksize),zero/)
     916                 :      192072 :        call abi_xcopy(1,zvar,1,lambda(cplx*(iblocksize-1)+1,iblocksize),1,x_cplx=x_cplx)
     917                 :             : #endif
     918                 :             :      end do
     919                 :             : 
     920                 :             : !    DEBUG
     921                 :             : !    write(std_out,*)'eigen',eigen(1:blocksize)
     922                 :             : !    ENDDEBUG
     923                 :             : 
     924                 :             : !    coordx(1:bigorder*cplx,1:blocksize)=grama(1:bigorder*cplx,1:blocksize)
     925   [ +  +  +  + ]:     2833751 :      coordx1(:,:) =  grama(1+i1*cplx : i2*cplx,1:blocksize)
     926   [ +  +  +  + ]:     2833751 :      coordx2(:,:) =  grama(1+i2*cplx : i3*cplx,1:blocksize)
     927         [ +  + ]:       28915 :      if(bigorder==i4) then
     928   [ +  +  +  + ]:     2402221 :        coordx3(:,:) =  grama(1+i3*cplx : i4*cplx,1:blocksize)
     929                 :             :      end if
     930                 :             : 
     931                 :             : 
     932         [ -  + ]:       28915 :      if(use_linalg_gpu==1) then
     933                 :           0 :        call copy_on_gpu(coordx2, coordx2_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     934         [ #  # ]:           0 :        if(bigorder==i4) then
     935                 :           0 :          call copy_on_gpu(coordx3, coordx3_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
     936                 :             :        end if
     937                 :             :      end if
     938                 :             : 
     939         [ -  + ]:       28915 :      ABI_FREE(grama)
     940         [ -  + ]:       28915 :      ABI_FREE(gramb)
     941         [ -  + ]:       28915 :      ABI_FREE(eigen)
     942   [ +  +  +  - ]:       28915 :      if (restart==0 .and. iterationnumber >1) then
     943                 :             : 
     944                 :             : !      blockvectorp=matmul(blockvectorr,coordx(i2+1:i3,:))+&
     945                 :             : !      &               matmul(blockvectorp,coordx(i3+1:i4,:))
     946         [ -  + ]:       23283 :        if(use_linalg_gpu==1) then
     947                 :             : !        call copy_on_gpu(blockvectorr, blockvectorr_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     948                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,blockvectorr_gpu,&
     949                 :           0 : &         vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
     950                 :           0 :          call copy_on_gpu(blockvectorp, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     951                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,vectsize,&
     952                 :           0 : &         coordx3_gpu,blocksize,cone,C_gpu,vectsize)
     953                 :           0 :          call copy_from_gpu(blockvectorp, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     954                 :             :        else
     955                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorr,&
     956                 :       23283 : &         vectsize,coordx2,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
     957                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorp,&
     958                 :       23283 : &         vectsize,coordx3,blocksize,cone,blockvectordumm,vectsize,x_cplx=x_cplx)
     959   [ -  +  -  -  :    27286240 :          blockvectorp = blockvectordumm
             +  +  +  + ]
     960                 :             :        end if
     961                 :             : 
     962                 :             : !      blockvectorap=matmul(blockvectorar,coordx(i2+1:i3,:))+&
     963                 :             : !      &                matmul(blockvectorap,coordx(i3+1:i4,:))
     964         [ -  + ]:       23283 :        if(use_linalg_gpu==1) then
     965                 :             : !        call copy_on_gpu(blockvectorar,blockvectorar_gpu,cplx*dp*vectsize*blocksize)
     966                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,blockvectorar_gpu,&
     967                 :           0 : &         vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
     968                 :           0 :          call copy_on_gpu(blockvectorap, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     969                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,vectsize,&
     970                 :           0 : &         coordx3_gpu,blocksize,cone,C_gpu,vectsize)
     971                 :           0 :          call copy_from_gpu(blockvectorap, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     972                 :             :        else
     973                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorar,&
     974                 :       23283 : &         vectsize,coordx2,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
     975                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorap,&
     976                 :       23283 : &         vectsize,coordx3,blocksize,cone,blockvectordumm,vectsize,x_cplx=x_cplx)
     977   [ -  +  -  -  :    27286240 :          blockvectorap = blockvectordumm
             +  +  +  + ]
     978                 :             :        end if
     979                 :             : 
     980                 :             : 
     981                 :             : !      blockvectorvp=matmul(blockvectorvr,coordx(i2+1:i3,:))+&
     982                 :             : !      &                matmul(blockvectorvp,coordx(i3+1:i4,:))
     983         [ +  + ]:       23283 :        if (gs_hamk%usepaw==0) then
     984         [ -  + ]:        4722 :          if(use_linalg_gpu==1) then
     985                 :           0 :            call copy_on_gpu(blockvectorvr, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     986                 :             :            call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,&
     987                 :           0 : &           vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
     988                 :           0 :            call copy_on_gpu(blockvectorvp, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     989                 :             :            call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,&
     990                 :           0 : &           vectsize,coordx3_gpu,blocksize,cone,C_gpu,vectsize)
     991                 :           0 :            call copy_from_gpu(blockvectorvp, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
     992                 :             :          else
     993                 :             :            call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorvr,&
     994                 :        4722 : &           vectsize,coordx2,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
     995                 :             :            call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorvp,&
     996                 :        4722 : &           vectsize,coordx3,blocksize,cone,blockvectordumm,vectsize,x_cplx=x_cplx)
     997   [ +  -  +  -  :     3491147 :            blockvectorvp = blockvectordumm
          -  +  -  -  +  
                +  +  + ]
     998                 :             :          end if
     999                 :             :        end if
    1000                 :             : 
    1001                 :             : !      blockvectorbp=matmul(blockvectorbr,coordx(i2+1:i3,:))+&
    1002                 :             : !      &                matmul(blockvectorbp,coordx(i3+1:i4,:))
    1003         [ -  + ]:       23283 :        if(use_linalg_gpu==1) then
    1004                 :             : !        call copy_on_gpu(blockvectorbr,blockvectorbr_gpu,cplx*dp*vectsize*blocksize)
    1005                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,blockvectorbr_gpu,&
    1006                 :           0 : &         vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1007                 :           0 :          call copy_on_gpu(blockvectorbp, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1008                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,vectsize,&
    1009                 :           0 : &         coordx3_gpu,blocksize,cone,C_gpu,vectsize)
    1010                 :           0 :          call copy_from_gpu(blockvectorbp, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1011                 :             :        else
    1012                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorbr,&
    1013                 :       23283 : &         vectsize,coordx2,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
    1014                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorbp,&
    1015                 :       23283 : &         vectsize,coordx3,blocksize,cone,blockvectordumm,vectsize,x_cplx=x_cplx)
    1016   [ -  +  -  -  :    27286240 :          blockvectorbp = blockvectordumm
             +  +  +  + ]
    1017                 :             :        end if
    1018                 :             : 
    1019                 :             :      else
    1020                 :             : 
    1021                 :             : !      blockvectoSz =matmul(blockvectorr,coordx(i2+1:i3,:))
    1022         [ -  + ]:        5632 :        if(use_linalg_gpu==1) then
    1023                 :             : !        call copy_on_gpu(blockvectorr,blockvectorr_gpu,cplx*dp*vectsize*blocksize)
    1024                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,blockvectorr_gpu,&
    1025                 :           0 : &         vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1026                 :           0 :          call copy_from_gpu(blockvectorp, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1027                 :             :        else
    1028                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorr,&
    1029                 :        5632 : &         vectsize,coordx2,blocksize,czero,blockvectorp,vectsize,x_cplx=x_cplx)
    1030                 :             :        end if
    1031                 :             : 
    1032                 :             : !      blockvectorap=matmul(blockvectorar,coordx(i2+1:i3,:))
    1033         [ -  + ]:        5632 :        if(use_linalg_gpu==1) then
    1034                 :             : !        call copy_on_gpu(blockvectorar,blockvectorar_gpu,cplx*dp*vectsize*blocksize)
    1035                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,blockvectorar_gpu,&
    1036                 :           0 : &         vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1037                 :           0 :          call copy_from_gpu(blockvectorap, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1038                 :             :        else
    1039                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorar,&
    1040                 :        5632 : &         vectsize,coordx2,blocksize,czero,blockvectorap,vectsize,x_cplx=x_cplx)
    1041                 :             :        end if
    1042                 :             : !      blockvectorvp=matmul(blockvectorvr,coordx(i2+1:i3,:))
    1043         [ +  + ]:        5632 :        if (gs_hamk%usepaw==0) then
    1044         [ -  + ]:         828 :          if(use_linalg_gpu==1) then
    1045                 :           0 :            call copy_on_gpu(blockvectorvr, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1046                 :             :            call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,&
    1047                 :           0 : &           vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1048                 :           0 :            call copy_from_gpu(blockvectorvp, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1049                 :             :          else
    1050                 :             :            call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorvr,&
    1051                 :         828 : &           vectsize,coordx2,blocksize,czero,blockvectorvp,vectsize,x_cplx=x_cplx)
    1052                 :             :          end if
    1053                 :             :        end if
    1054                 :             : 
    1055                 :             : !      blockvectorbp=matmul(blockvectorbr,coordx(i2+1:i3,:))
    1056         [ -  + ]:        5632 :        if(use_linalg_gpu==1) then
    1057                 :             : !        call copy_on_gpu(blockvectorbr,blockvectorbr_gpu,cplx*dp*vectsize*blocksize)
    1058                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,blockvectorbr_gpu,&
    1059                 :           0 : &         vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1060                 :           0 :          call copy_from_gpu(blockvectorbp, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1061                 :             :        else
    1062                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorbr,&
    1063                 :        5632 : &         vectsize,coordx2,blocksize,czero,blockvectorbp,vectsize,x_cplx=x_cplx)
    1064                 :             :        end if
    1065                 :             :      end if
    1066                 :             : 
    1067         [ +  - ]:       28915 :      if(use_linalg_gpu==1) then
    1068                 :           0 :        call copy_on_gpu(coordx1, coordx2_gpu, INT(cplx, c_size_t)*dp*blocksize*blocksize)
    1069                 :             :      end if
    1070                 :             : 
    1071                 :             : !    blockvectorx = matmul(blockvectorx,coordx(i1+1:i2,:))+blockvectorp
    1072                 :             :      if(use_linalg_gpu==1) then
    1073                 :           0 :        call copy_on_gpu(blockvectorx, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1074                 :             :        call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,&
    1075                 :           0 : &       vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1076                 :           0 :        call copy_from_gpu(blockvectordumm, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1077                 :             :      else
    1078                 :             :        call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorx,&
    1079                 :       28915 : &       vectsize,coordx1,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
    1080                 :             :      end if
    1081   [ +  -  +  -  :    33616082 :      blockvectorx = blockvectordumm+blockvectorp
          -  +  -  -  +  
                +  +  + ]
    1082                 :             : 
    1083                 :             : !    blockvectorax= matmul(blockvectorax,coordx(i1+1:i2,:))+blockvectorap
    1084         [ -  + ]:       28915 :      if(use_linalg_gpu==1) then
    1085                 :           0 :        call copy_on_gpu(blockvectorax, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1086                 :             :        call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,&
    1087                 :           0 : &       vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1088                 :           0 :        call copy_from_gpu(blockvectordumm, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1089                 :             :      else
    1090                 :             :        call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorax,&
    1091                 :       28915 : &       vectsize,coordx1,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
    1092                 :             :      end if
    1093   [ +  -  +  -  :    33616082 :      blockvectorax = blockvectordumm+blockvectorap
          -  +  -  -  +  
                +  +  + ]
    1094                 :             : 
    1095                 :             : !    blockvectorvx= matmul(blockvectorvx,coordx(i1+1:i2,:))+blockvectorvp
    1096         [ +  + ]:       28915 :      if (gs_hamk%usepaw==0) then
    1097         [ -  + ]:        5550 :        if(use_linalg_gpu==1) then
    1098                 :           0 :          call copy_on_gpu(blockvectorvx, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1099                 :             :          call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,&
    1100                 :           0 : &         vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1101                 :           0 :          call copy_from_gpu(blockvectordumm, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1102                 :             :        else
    1103                 :             :          call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorvx,&
    1104                 :        5550 : &         vectsize,coordx1,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
    1105                 :             :        end if
    1106   [ +  -  +  -  :     4113243 :        blockvectorvx = blockvectordumm+blockvectorvp
          -  +  -  -  +  
                +  +  + ]
    1107                 :             :      end if
    1108                 :             : 
    1109                 :             : !    blockvectorbx= matmul(blockvectorbx,coordx(i1+1:i2,:))+blockvectorbp
    1110         [ -  + ]:       28915 :      if(use_linalg_gpu==1) then
    1111                 :           0 :        call copy_on_gpu(blockvectorbx, A_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1112                 :             :        call gpu_xgemm(cplx,'n','n',vectsize,blocksize,blocksize,cone,A_gpu,&
    1113                 :           0 : &       vectsize,coordx2_gpu,blocksize,czero,C_gpu,vectsize)
    1114                 :           0 :        call copy_from_gpu(blockvectordumm, C_gpu, INT(cplx, c_size_t)*dp*vectsize*blocksize)
    1115                 :             :      else
    1116                 :             :        call abi_xgemm('n','n',vectsize,blocksize,blocksize,cone,blockvectorbx,&
    1117                 :       28915 : &       vectsize,coordx1,blocksize,czero,blockvectordumm,vectsize,x_cplx=x_cplx)
    1118                 :             :      end if
    1119   [ +  -  +  -  :    33616082 :      blockvectorbx = blockvectordumm+blockvectorbp
          -  +  -  -  +  
                +  +  + ]
    1120                 :             : 
    1121                 :             : !    ABI_FREE(coordx)
    1122         [ -  + ]:       28915 :      ABI_FREE(coordx1)
    1123         [ -  + ]:       28915 :      ABI_FREE(coordx2)
    1124         [ +  + ]:       28915 :      if(bigorder==i4) then
    1125         [ -  + ]:       23283 :        ABI_FREE(coordx3)
    1126                 :             :      end if
    1127                 :             : 
    1128                 :             : !    Check convergence on energy and eventually exit
    1129         [ +  + ]:       28915 :      if (iterationnumber==1) then
    1130                 :        5632 :        deold=deltae
    1131         [ +  - ]:       23283 :      else if (iterationnumber>1) then
    1132   [ +  +  +  + ]:       23283 :        if ((abs(deltae)<0.005*abs(deold)).and.(iterationnumber/=maxiterations))then
    1133         [ -  + ]:        1077 :          if(prtvol>=10)then
    1134                 :             :            write(message, '(2(a,i4),1x,a,1p,e12.4,a,e12.4,a)' ) &
    1135                 :           0 : &           ' lobpcgwf: block',iblock,', line',iterationnumber,&
    1136                 :           0 : &           ', deltae=',deltae,' < 0.005*',deold,' =>skip lines !'
    1137                 :           0 :            call wrtout(std_out,message,'PERS')
    1138                 :             :          end if
    1139                 :             :          exit
    1140         [ +  + ]:       22206 :        else if (abs(deltae)>0.005*abs(deold)) then
    1141         [ -  + ]:       21669 :          if(prtvol>=10)then
    1142                 :             :            write(message, '(2(a,i4),1x,a,1p,e12.4,a,e12.4,a)' ) &
    1143                 :           0 : &           ' lobpcgwf: block',iblock,', line',iterationnumber,&
    1144                 :           0 : &           ', deltae=',deltae,' > 0.005*',deold,' =>keep on working !'
    1145                 :           0 :            call wrtout(std_out,message,'PERS')
    1146                 :             :          end if
    1147                 :             :        end if
    1148                 :             :      end if
    1149                 :             : 
    1150         [ -  + ]:      119142 :      if(abs(dtset%timopt)==4) then
    1151                 :           0 :        call timab(524,2,tsec)
    1152                 :             :      end if
    1153                 :             : 
    1154                 :             :    end do iter
    1155                 :             : 
    1156                 :             : !  ###########################################################################
    1157                 :             : !  ################## END LOOP ON NLINE ######################################
    1158                 :             : !  ###########################################################################
    1159                 :             : 
    1160         [ -  + ]:        5636 :    if(abs(dtset%timopt)==4) then
    1161                 :           0 :      call timab(525,1,tsec)
    1162                 :             :    end if
    1163                 :             : 
    1164         [ +  + ]:        5636 :    if (havetoprecon) then
    1165         [ +  - ]:        5628 :      if (gs_hamk%use_gbt == 0) then
    1166                 :             :        call xprecon(blockvectorbx,lambda,blocksize,&
    1167                 :             : &       iterationnumber,kinpw,mpi_enreg,npw_k,my_nspinor,&
    1168                 :        5628 : &       optekin,optpcon,pcon,blockvectorax,blockvectorr,vectsize,timopt=timopt,tim_xprecon=tim_xprecon)
    1169                 :             :      else
    1170                 :             :        call xprecon(blockvectorbx,lambda,blocksize,&
    1171                 :             : &       iterationnumber,kinpw,mpi_enreg,npw_k,my_nspinor,&
    1172                 :             : &       optekin,optpcon,pcon,blockvectorax,blockvectorr,vectsize,timopt=timopt,tim_xprecon=tim_xprecon,&
    1173   [ #  #  #  #  :           0 : &       kinpw_kp=kinpw_kp,pcon_kp=pcon_kp)
             #  #  #  # ]
    1174                 :             :      end if
    1175                 :             : 
    1176   [ +  -  -  +  :     6325490 :      residualnorms=sum(abs(blockvectorr)**2,dim=1)
          -  -  +  +  +  
                      + ]
    1177                 :             : 
    1178         [ -  + ]:        5628 :      if(abs(dtset%timopt)==3) then
    1179                 :           0 :        call timab(533,1,tsec)
    1180                 :             :      end if
    1181                 :        5628 :      call xmpi_sum(residualnorms,mpi_enreg%comm_bandspinorfft,ierr)
    1182         [ -  + ]:        5628 :      if(abs(dtset%timopt)==3) then
    1183                 :           0 :        call timab(533,2,tsec)
    1184                 :             :      end if
    1185                 :             : 
    1186         [ +  + ]:       32304 :      resid_k(bblocksize+1:bblocksize+blocksize)=residualnorms(1:blocksize)
    1187                 :             :    end if
    1188                 :             : 
    1189                 :             :    call wfcopy('I',vectsize*blocksize,blockvectorx,1,cg,1,blocksize,iblock,'C',withbbloc=.true.,&
    1190                 :        5636 : &   timopt=timopt,tim_wfcopy=tim_wfcopy)
    1191                 :             : 
    1192         [ +  + ]:        5636 :    if(gen_eigenpb) then
    1193                 :             :      call wfcopy('I',vectsize*blocksize,blockvectorbx,1,gsc,1,blocksize,iblock,'S',withbbloc=.true.,&
    1194                 :        4804 : &     timopt=timopt,tim_wfcopy=tim_wfcopy)
    1195                 :             :    end if
    1196                 :             : 
    1197                 :             : !  The Vnl+VFockACE part of the Hamiltonian is no more stored in the packed form such as it was the case for subvnlx(:).
    1198                 :             : !  Now, the full matrix is stored in totvnlx(:,:). This trick permits:
    1199                 :             : !  1) to avoid the reconstruction of the total matrix in vtowfk.F90 (double loop over bands)
    1200                 :             : !  2) to use two optimized matrix-matrix blas routine for general (in lobpcgccwf.F90) or hermitian (in vtowfk.F90)
    1201                 :             : !  operators, zgemm.f and zhemm.f respectively, rather than a triple loop in both cases.
    1202                 :        5636 :    iwavef=iblock*blocksize
    1203                 :        5636 :    isubh=1+2*bblocksize*(bblocksize+1)/2
    1204                 :             : 
    1205   [ -  +  +  -  :       22544 :    ABI_MALLOC(blockvectorz,(cplx*vectsize,iwavef))
          -  +  +  -  +  
          -  +  -  -  +  
                   -  + ]
    1206         [ +  + ]:        5636 :    if(bblocksize > 0 ) then
    1207   [ +  +  +  + ]:   153725699 :      blockvectorz(:,1:bblocksize) = blockvectory(:,1:bblocksize)
    1208                 :             :    end if
    1209   [ +  +  +  + ]:     6328558 :    blockvectorz(:,bblocksize+1:iwavef) = blockvectorx(:,1:blocksize)
    1210                 :             : 
    1211   [ -  +  +  -  :       22544 :    ABI_MALLOC(tsubham,(cplx*iwavef,blocksize))
          -  +  +  -  +  
          -  +  -  -  +  
                   -  + ]
    1212   [ +  +  +  + ]:     1700447 :    tsubham(:,:)=zero
    1213                 :             :    call abi_xgemm(cparam(cplx),'n',iwavef,blocksize,vectsize,cone,blockvectorz,vectsize,&
    1214                 :        5636 : &   blockvectorax,vectsize,czero,tsubham,iwavef,x_cplx=x_cplx)
    1215                 :             : 
    1216         [ +  + ]:        5636 :    if (gs_hamk%usepaw==0) then
    1217                 :             :      ! MG FIXME: Here gfortran4.9 allocates temporary array for C in abi_d2zgemm.
    1218                 :             :      call abi_xgemm(cparam(cplx),'n',blocksize,iwavef,vectsize,cone,blockvectorvx,vectsize,&
    1219   [ +  +  +  -  :        8560 : &     blockvectorz,vectsize,czero,totvnlx(cplx*bblocksize+1:cplx*iwavef,1:iwavef),blocksize,x_cplx=x_cplx)
          +  +  +  +  +  
             +  +  +  +  
                      + ]
    1220                 :             :    end if
    1221                 :             : 
    1222         [ +  + ]:       32336 :    do iblocksize=1,blocksize
    1223         [ +  + ]:      782714 :      do ii=1,bblocksize+iblocksize
    1224         [ +  + ]:      750378 :        if ( cplx == 1 ) then
    1225                 :        6297 :          subham(isubh)  = tsubham(ii,iblocksize)
    1226                 :        6297 :          subham(isubh+1)= zero
    1227                 :             :        else
    1228                 :      744081 :          subham(isubh)  = tsubham(2*ii-1,iblocksize)
    1229                 :      744081 :          subham(isubh+1)= tsubham(2*ii  ,iblocksize)
    1230                 :             :        end if
    1231                 :      777078 :        isubh=isubh+2
    1232                 :             :      end do
    1233                 :             :    end do
    1234         [ -  + ]:        5636 :    ABI_FREE(tsubham)
    1235         [ -  + ]:        5636 :    ABI_FREE(blockvectorz)
    1236                 :             : !  comm for subham and subvnlx are made in vtowfk
    1237                 :             : 
    1238         [ -  + ]:        5636 :    ABI_FREE(pcon)
    1239         [ -  + ]:        5636 :    ABI_SFREE(pcon_kp)
    1240         [ -  + ]:        5636 :    ABI_FREE(blockvectory)
    1241         [ -  + ]:        5636 :    ABI_FREE(blockvectorby)
    1242         [ -  + ]:        5636 :    ABI_FREE(gramyx)
    1243         [ -  + ]:        5636 :    ABI_FREE(blockvectorx)
    1244         [ -  + ]:        5636 :    ABI_FREE(blockvectorax)
    1245         [ -  + ]:        5636 :    ABI_FREE(blockvectorbx)
    1246         [ -  + ]:        5636 :    ABI_FREE(blockvectorr)
    1247         [ -  + ]:        5636 :    ABI_FREE(blockvectorar)
    1248         [ -  + ]:        5636 :    ABI_FREE(blockvectorbr)
    1249         [ -  + ]:        5636 :    ABI_FREE(blockvectorp)
    1250         [ -  + ]:        5636 :    ABI_FREE(blockvectorap)
    1251         [ -  + ]:        5636 :    ABI_FREE(blockvectorbp)
    1252         [ +  + ]:        5636 :    if (gs_hamk%usepaw==0) then
    1253         [ -  + ]:         832 :      ABI_FREE(blockvectorvx)
    1254         [ -  + ]:         832 :      ABI_FREE(blockvectorvp)
    1255         [ -  + ]:         832 :      ABI_FREE(blockvectorvr)
    1256                 :             :    end if
    1257         [ -  + ]:        5636 :    ABI_FREE(blockvectordumm)
    1258         [ -  + ]:        5636 :    ABI_FREE(gramxax)
    1259         [ -  + ]:        5636 :    ABI_FREE(gramxar)
    1260         [ -  + ]:        5636 :    ABI_FREE(gramxap)
    1261         [ -  + ]:        5636 :    ABI_FREE(gramrar)
    1262         [ -  + ]:        5636 :    ABI_FREE(gramrap)
    1263         [ -  + ]:        5636 :    ABI_FREE(grampap)
    1264         [ -  + ]:        5636 :    ABI_FREE(gramxbx)
    1265         [ -  + ]:        5636 :    ABI_FREE(gramxbr)
    1266         [ -  + ]:        5636 :    ABI_FREE(gramxbp)
    1267         [ -  + ]:        5636 :    ABI_FREE(gramrbr)
    1268         [ -  + ]:        5636 :    ABI_FREE(gramrbp)
    1269         [ -  + ]:        5636 :    ABI_FREE(grampbp)
    1270         [ -  + ]:        5636 :    ABI_FREE(transf3)
    1271         [ -  + ]:        5636 :    ABI_FREE(transf5)
    1272         [ -  + ]:        5636 :    ABI_FREE(lambda)
    1273         [ -  + ]:        5636 :    ABI_FREE(residualnorms)
    1274         [ -  + ]:       13522 :    if(use_linalg_gpu==1) then
    1275                 :           0 :      call dealloc_on_gpu(bblockvector_gpu)
    1276                 :           0 :      call dealloc_on_gpu(gram_gpu)
    1277                 :             :    end if
    1278                 :             : 
    1279                 :             :  end do  ! End big loop over bands inside blocks
    1280                 :             : 
    1281                 :             : #ifdef HAVE_OPENMP_OFFLOAD
    1282                 :             :  !$OMP TARGET EXIT DATA MAP(delete:cwavef,gwavef,gvnlxc,swavef) IF(dtset%gpu_option==ABI_GPU_OPENMP)
    1283                 :             : #endif
    1284         [ -  + ]:        2250 :  ABI_FREE(cwavef)
    1285         [ -  + ]:        2250 :  ABI_FREE(gwavef)
    1286         [ -  + ]:        2250 :  ABI_FREE(gvnlxc)
    1287         [ -  + ]:        2250 :  ABI_FREE(swavef)
    1288   [ +  -  +  +  :        6908 :  ABI_FREE(cprj_dum)
          -  +  -  +  -  
                      + ]
    1289                 :             : 
    1290         [ -  + ]:        2250 :  if(use_linalg_gpu==1) then
    1291                 :           0 :    call dealloc_on_gpu(blockvectorr_gpu)
    1292                 :           0 :    call dealloc_on_gpu(blockvectorar_gpu)
    1293                 :           0 :    call dealloc_on_gpu(blockvectorbr_gpu)
    1294                 :           0 :    call dealloc_on_gpu(A_gpu)
    1295                 :           0 :    call dealloc_on_gpu(C_gpu)
    1296                 :           0 :    call dealloc_on_gpu(coordx2_gpu)
    1297                 :           0 :    call dealloc_on_gpu(coordx3_gpu)
    1298                 :             :    !call gpu_linalg_shutdown()
    1299                 :             :  end if
    1300                 :             : 
    1301         [ -  + ]:        2250 :  if(abs(dtset%timopt)==4) then
    1302                 :           0 :    call timab(525,2,tsec)
    1303                 :             :  end if
    1304   [ -  +  -  +  :        2250 :  call timab(530,2,tsec)
          -  +  -  +  -  
          +  -  +  -  +  
          -  +  -  +  -  
          +  -  -  -  -  
          -  -  -  +  -  
          +  -  +  -  +  
          -  +  -  +  -  
          +  -  +  -  +  
          -  +  -  +  -  
                      + ]
    1305                 :             : 
    1306                 :             :  DBG_ENTER("COLL")
    1307                 :             : 
    1308                 :             :  contains
    1309                 :             : 
    1310                 :       28915 :    function gramindex(iblocksize)
    1311                 :             : 
    1312                 :             :    integer :: gramindex,iblocksize
    1313         [ +  + ]:      269780 :    gramindex=(iblocksize-1)*cplx+1
    1314                 :        2250 :  end function gramindex
    1315                 :             : 
    1316                 :             : end subroutine lobpcgwf
    1317                 :             : !!***
    1318                 :             : 
    1319                 :             : end module m_lobpcgwf_old
    1320                 :             : !!***
        

Generated by: LCOV version 2.3-1