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 : : !! 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 : : !!***
|
