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