Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 :
6 : #include "abi_common.h"
7 :
8 : module m_lobpcg2_cprj
9 :
10 : use defs_basis
11 : use m_xg
12 : use m_xg_nonlop
13 : use m_xg_ortho_RR
14 : use m_xgTransposer
15 : use m_xgScalapack
16 : use defs_basis
17 : use m_abicore
18 : use m_errors
19 : use m_xomp
20 : #ifdef HAVE_OPENMP
21 : use omp_lib
22 : #endif
23 :
24 : use m_xmpi
25 : use m_time, only : timab
26 :
27 : implicit none
28 :
29 : private
30 :
31 : integer, parameter :: tim_init = 2040
32 : integer, parameter :: tim_free = 2041
33 : integer, parameter :: tim_cprj = 2043
34 : integer, parameter :: tim_ortho = 2044
35 : integer, parameter :: tim_transpose = 2039
36 : integer, parameter :: tim_maxres = 2046
37 : integer, parameter :: tim_ax_k = 2048
38 : integer, parameter :: tim_ax_v = 2049
39 : integer, parameter :: tim_ax_nl = 2050
40 : integer, parameter :: tim_pcond = 2047
41 :
42 : integer, parameter :: tim_Bortho_X = 2031
43 : integer, parameter :: tim_Bortho_XW = 2032
44 : integer, parameter :: tim_Bortho_XWP = 2033
45 : integer, parameter :: tim_Bortho_Xall = 2034
46 : integer, parameter :: tim_RR_X = 2035
47 : integer, parameter :: tim_RR_XW = 2036
48 : integer, parameter :: tim_RR_XWP = 2037
49 : integer, parameter :: tim_RR_Xall = 2038
50 :
51 : integer, parameter :: tim_copy = 2042
52 : integer, parameter :: tim_nbdbuf = 2045
53 :
54 : integer, parameter :: tim_enl = 2051
55 :
56 : type, public :: lobpcg_t
57 : logical :: is_nested ! For OpenMP nested region
58 : integer :: spacedim ! Space dimension for one vector (WFs)
59 : integer :: space ! Space type for one vector
60 : integer :: space_cprj ! Space type for one projector
61 : integer :: cprjdim ! Space dimension for one vector (cprj)
62 : integer :: neigenpairs ! Number of eigen values/vectors we want
63 : integer :: nblock ! Number of block in the space of dim spacedim
64 : integer :: blockdim ! Number of vectors in one block
65 : integer :: blockdim_cprj ! Number of vectors in one block (for cprj array)
66 : integer :: nline ! Number of line to perform
67 : integer :: spacecom ! Communicator for MPI
68 : integer :: paral_kgb ! paral_kgb formalism or not
69 : integer :: comm_rows ! communicator for rows
70 : integer :: comm_cols ! communicator for cols
71 : integer :: me_g0
72 : integer :: me_g0_fft
73 : double precision :: tolerance ! Tolerance on the residu to stop the minimization
74 : integer :: prtvol
75 : type(xg_nonlop_t) :: xg_nonlop ! contains data for nonlocal operations
76 : type(xgBlock_t) :: AllX0 ! Block of initial and final solution.
77 : type(xgBlock_t) :: AllcprjX0 ! Block of initial and final solution.
78 : type(xgBlock_t) :: X0
79 : type(xgBlock_t) :: cprjX0
80 : ! Dim is (cplx*spacedim,neigenpair)
81 : !double precision, allocatable :: eig(:) ! Actual eigen values
82 :
83 : ! Some denomination here
84 : ! X, W, P are the shifts where their block vectors starts.
85 : ! X is the current block we are solving
86 : ! W is the residu of the pencil A-lambda*B ie W=A*X-lambdaB*X
87 : ! P correspond to X_{iline} - X_{iline-1}
88 : ! if prefixed by A (AX, AW, AP) it means the result of A*X (A*W, A*P ! respectively)
89 : ! Idem for B (BX, BW, BP) which is the result of B*X (B*W, B*P)
90 : ! An extend to this is AXWP which is the vector [AX,AW,AP]
91 : ! and BXWP which is [BX,BW,BP] columnwise
92 : ! Now we only allocate the triple XWP matrix and each individual X, W and P
93 : ! Example : XWP=[x1,x2,x3,w1,w2,w3,p1,p2,p3] X points to the indice of x1-1, W points to w1-1 and
94 : ! P to p1-1
95 : type(xg_t) :: XWP
96 : type(xg_t) :: AXWP
97 : type(xg_t) :: cprjXWP
98 :
99 : type(xg_t) :: cprj_work
100 : type(xg_t) :: cprj_work2
101 :
102 : type(xgBlock_t) :: XColsRows
103 : type(xgBlock_t) :: AXColsRows
104 :
105 : type(xgTransposer_t) :: xgTransposerX
106 : type(xgTransposer_t) :: xgTransposerAX
107 :
108 : type(xgBlock_t) :: WColsRows
109 : type(xgBlock_t) :: AWColsRows
110 :
111 : type(xgTransposer_t) :: xgTransposerW
112 : type(xgTransposer_t) :: xgTransposerAW
113 :
114 : type(xgBlock_t) :: X ! Shift to apply to start reading the X values
115 : type(xgBlock_t) :: W ! Shift to apply to start reading the W values
116 : type(xgBlock_t) :: P ! Shift to apply to start reading the P values
117 : type(xgBlock_t) :: XW ! Shift to apply to start reading the P values
118 : type(xgBlock_t) :: WP ! Shift to apply to start reading the P values
119 :
120 : type(xgBlock_t) :: AX ! Shift to apply to start reading the X values
121 : type(xgBlock_t) :: AW ! Shift to apply to start reading the W values
122 : type(xgBlock_t) :: AP ! Shift to apply to start reading the P values
123 : type(xgBlock_t) :: AXW ! Shift to apply to start reading the P values
124 : type(xgBlock_t) :: AWP ! Shift to apply to start reading the P values
125 :
126 : type(xgBlock_t) :: cprjX ! Shift to apply to start reading the X values
127 : type(xgBlock_t) :: cprjW ! Shift to apply to start reading the W values
128 : type(xgBlock_t) :: cprjP ! Shift to apply to start reading the P values
129 : type(xgBlock_t) :: cprjXW ! Shift to apply to start reading the P values
130 : type(xgBlock_t) :: cprjWP ! Shift to apply to start reading the P values
131 :
132 : type(xg_t) :: AllAX0
133 :
134 : type(xgBlock_t) :: AllX0ColsRows
135 : type(xgBlock_t) :: AllAX0ColsRows
136 :
137 : type(xgTransposer_t) :: xgTransposerAllX0
138 : type(xgTransposer_t) :: xgTransposerAllAX0
139 :
140 : end type lobpcg_t
141 :
142 : public :: lobpcg_init
143 : public :: lobpcg_memInfo
144 : public :: lobpcg_run_cprj
145 : public :: lobpcg_free
146 :
147 : contains
148 :
149 :
150 4036 : subroutine lobpcg_init(lobpcg, bandpp, neigenpairs, spacedim, cprjdim, blockdim, tolerance, nline, &
151 : & space, space_cprj, spacecom, paral_kgb, xg_nonlop, comm_rows, comm_cols, me_g0, me_g0_fft)
152 :
153 : type(lobpcg_t) , intent(inout) :: lobpcg
154 : type(xg_nonlop_t), intent(in ) :: xg_nonlop
155 : integer , intent(in ) :: neigenpairs
156 : integer , intent(in ) :: spacedim
157 : integer , intent(in ) :: cprjdim
158 : integer , intent(in ) :: blockdim
159 : integer , intent(in ) :: bandpp
160 : double precision , intent(in ) :: tolerance
161 : integer , intent(in ) :: nline
162 : integer , intent(in ) :: comm_rows,comm_cols
163 : integer , intent(in ) :: space
164 : integer , intent(in ) :: space_cprj
165 : integer , intent(in ) :: spacecom
166 : integer , intent(in ) :: paral_kgb
167 : integer , intent(in ) :: me_g0
168 : integer , intent(in ) :: me_g0_fft
169 : double precision :: tsec(2)
170 :
171 4036 : call timab(tim_init,1,tsec)
172 4036 : lobpcg%neigenpairs = neigenpairs
173 4036 : lobpcg%spacedim = spacedim
174 4036 : lobpcg%space = space
175 4036 : lobpcg%space_cprj = space_cprj
176 4036 : lobpcg%cprjdim = cprjdim
177 4036 : lobpcg%nblock = neigenpairs / blockdim
178 4036 : lobpcg%blockdim = blockdim
179 4036 : lobpcg%blockdim_cprj = bandpp*xg_nonlop%nspinor
180 4036 : if (tolerance > 0.0) then
181 0 : lobpcg%tolerance = tolerance
182 : else
183 4036 : lobpcg%tolerance = 1e-20
184 : end if
185 4036 : lobpcg%nline = nline
186 4036 : lobpcg%spacecom = spacecom
187 4036 : lobpcg%paral_kgb = paral_kgb
188 4036 : lobpcg%comm_rows = comm_rows
189 4036 : lobpcg%comm_cols = comm_cols
190 4036 : lobpcg%xg_nonlop = xg_nonlop
191 4036 : lobpcg%me_g0 = me_g0
192 4036 : lobpcg%me_g0_fft = me_g0_fft
193 :
194 4036 : call lobpcg_allocateAll(lobpcg)
195 4036 : call timab(tim_init,2,tsec)
196 :
197 4036 : end subroutine lobpcg_init
198 :
199 :
200 4036 : subroutine lobpcg_allocateAll(lobpcg)
201 :
202 : type(lobpcg_t) , intent(inout) :: lobpcg
203 : integer :: spacedim
204 : integer :: cprjdim
205 : integer :: blockdim
206 : integer :: blockdim_cprj
207 :
208 4036 : spacedim = lobpcg%spacedim
209 4036 : blockdim = lobpcg%blockdim
210 :
211 4036 : cprjdim = lobpcg%cprjdim
212 4036 : blockdim_cprj = lobpcg%blockdim_cprj
213 :
214 4036 : call lobpcg_free(lobpcg) ! Make sure everything is not allocated and
215 : ! pointer point to null()
216 :
217 4036 : call xg_init(lobpcg%XWP,lobpcg%space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=lobpcg%me_g0)
218 4036 : call xg_setBlock(lobpcg%XWP,lobpcg%X,spacedim,blockdim)
219 4036 : call xg_setBlock(lobpcg%XWP,lobpcg%W,spacedim,blockdim,fcol=blockdim+1)
220 4036 : call xg_setBlock(lobpcg%XWP,lobpcg%P,spacedim,blockdim,fcol=2*blockdim+1)
221 4036 : call xg_setBlock(lobpcg%XWP,lobpcg%XW,spacedim,2*blockdim)
222 4036 : call xg_setBlock(lobpcg%XWP,lobpcg%WP,spacedim,2*blockdim,fcol=blockdim+1)
223 :
224 4036 : call xg_init(lobpcg%AXWP,lobpcg%space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=lobpcg%me_g0)
225 4036 : call xg_setBlock(lobpcg%AXWP,lobpcg%AX,spacedim,blockdim)
226 4036 : call xg_setBlock(lobpcg%AXWP,lobpcg%AW,spacedim,blockdim,fcol=blockdim+1)
227 4036 : call xg_setBlock(lobpcg%AXWP,lobpcg%AP,spacedim,blockdim,fcol=2*blockdim+1)
228 4036 : call xg_setBlock(lobpcg%AXWP,lobpcg%AXW,spacedim,2*blockdim)
229 4036 : call xg_setBlock(lobpcg%AXWP,lobpcg%AWP,spacedim,2*blockdim,fcol=blockdim+1)
230 :
231 4036 : call xg_init(lobpcg%cprjXWP,lobpcg%space_cprj,cprjdim,3*blockdim_cprj,lobpcg%spacecom)
232 4036 : call xg_setBlock(lobpcg%cprjXWP,lobpcg%cprjX, cprjdim, blockdim_cprj)
233 4036 : call xg_setBlock(lobpcg%cprjXWP,lobpcg%cprjW, cprjdim, blockdim_cprj,fcol= blockdim_cprj+1)
234 4036 : call xg_setBlock(lobpcg%cprjXWP,lobpcg%cprjP, cprjdim, blockdim_cprj,fcol= 2*blockdim_cprj+1)
235 4036 : call xg_setBlock(lobpcg%cprjXWP,lobpcg%cprjXW,cprjdim,2*blockdim_cprj)
236 4036 : call xg_setBlock(lobpcg%cprjXWP,lobpcg%cprjWP,cprjdim,2*blockdim_cprj,fcol= blockdim_cprj+1)
237 :
238 4036 : call xg_init(lobpcg%cprj_work ,lobpcg%space_cprj,cprjdim,blockdim_cprj,lobpcg%spacecom)
239 4036 : call xg_init(lobpcg%cprj_work2,lobpcg%space_cprj,cprjdim,blockdim_cprj,lobpcg%spacecom)
240 :
241 4036 : if ( lobpcg%nblock /= 1 ) then
242 4004 : call xg_init(lobpcg%AllAX0,lobpcg%space,spacedim,lobpcg%neigenpairs,lobpcg%spacecom,me_g0=lobpcg%me_g0)
243 : else
244 32 : lobpcg%AllAX0%self = lobpcg%AX
245 : end if
246 :
247 4036 : end subroutine lobpcg_allocateAll
248 :
249 :
250 0 : function lobpcg_memInfo(neigenpairs, spacedim, blockdim, space) result(arraymem)
251 :
252 : integer , intent(in ) :: neigenpairs
253 : integer , intent(in ) :: spacedim
254 : integer , intent(in ) :: blockdim
255 : integer , intent(in ) :: space
256 : double precision :: memXWP
257 : double precision :: memAXWP
258 : double precision :: memBXWP
259 : double precision :: memAllBX0
260 : double precision :: memAllAX0
261 : double precision :: memeigenvalues3N
262 : double precision :: membufferOrtho
263 : double precision :: membufferBOrtho
264 : double precision :: memsubA
265 : double precision :: memsubB
266 : double precision :: memsubBtmp
267 : double precision :: memvec
268 : double precision :: memRR
269 : double precision :: maxmemTmp
270 : double precision :: cplx
271 : integer :: nblock
272 : double precision :: arraymem(2)
273 :
274 0 : cplx = 1 ; if ( space == SPACE_C ) cplx = 2
275 0 : nblock = neigenpairs/blockdim
276 :
277 : ! Permanent in lobpcg
278 0 : memXWP = cplx* kind(1.d0) * spacedim * 3*blockdim
279 0 : memAXWP = cplx* kind(1.d0) * spacedim * 3*blockdim
280 0 : memBXWP = cplx* kind(1.d0) * spacedim * 3*blockdim
281 0 : if ( nblock > 1 ) then
282 0 : memAllAX0 = cplx * kind(1.d0) * spacedim * 3*blockdim
283 0 : memAllBX0 = cplx * kind(1.d0) * spacedim * 3*blockdim
284 0 : membufferOrtho = cplx * kind(1.d0) * blockdim * (nblock-1) * blockdim
285 : else
286 : memAllAX0 = 0
287 : memAllBX0 = 0
288 : membufferOrtho = 0
289 : endif
290 0 : memeigenvalues3N = kind(1.d0) * 3*blockdim
291 :
292 : ! Temporary arrays
293 0 : membufferBOrtho = cplx * kind(1.d0) * 2*blockdim * 2*blockdim ! For the moment being, only Bortho with X or WP at the same time
294 0 : memsubA = cplx * kind(1.d0) * 3*blockdim * 3*blockdim
295 0 : memsubB = cplx * kind(1.d0) * 3*blockdim * 3*blockdim
296 0 : memsubBtmp = cplx * kind(1.d0) * spacedim * blockdim
297 0 : memvec = cplx * kind(1.d0) * 3*blockdim * blockdim
298 0 : memRR = max(memsubA+memsubB+memvec,memsubBtmp+memvec)
299 :
300 0 : maxmemTmp = max( membufferBOrtho,memRR,membufferOrtho )
301 :
302 0 : arraymem(1) = memXWP+memAXWP+memBXWP+memAllAX0+memAllBX0+memeigenvalues3N
303 0 : arraymem(2) = maxmemTmp
304 :
305 0 : end function lobpcg_memInfo
306 :
307 :
308 12108 : subroutine lobpcg_run_cprj(lobpcg, X0, cprjX0, getAX, kin, pcond, eigen, occ, residu, enl, prtvol, nspinor, &
309 : isppol, ikpt, inonsc, istep, nbdbuf)
310 :
311 : type(lobpcg_t) , intent(inout) :: lobpcg
312 : type(xgBlock_t), intent(inout) :: X0 ! Full initial vectors
313 : type(xgBlock_t), intent(inout) :: cprjX0 ! Full initial cprj
314 : type(xgBlock_t), intent(inout) :: eigen ! Full initial eigen values
315 : type(xgBlock_t), intent(inout) :: occ
316 : type(xgBlock_t), intent(inout) :: residu
317 : type(xgBlock_t), intent(inout) :: enl
318 : type(xgBlock_t), intent(in ) :: kin
319 : type(xgBlock_t), intent(in ) :: pcond
320 : integer , intent(in ) :: prtvol
321 : integer , intent(in ) :: nspinor
322 : integer , intent(in ) :: isppol,ikpt,inonsc,istep,nbdbuf
323 :
324 : type(xg_t) :: eigenvalues3N ! eigen values for Rayleight-Ritz
325 : type(xg_t) :: residu_eff
326 : type(xgBlock_t) :: eigenvaluesN ! eigen values for Rayleight-Ritz
327 : type(xgBlock_t) :: eigenvalues2N ! eigen values for Rayleight-Ritz
328 : logical :: skip,compute_residu
329 : integer :: blockdim, blockdim3, blockdim2, blockdim_cprj
330 : integer :: spacedim
331 : integer :: iblock, nblock
332 : integer :: iline, nline
333 : integer :: rows_tmp, cols_tmp, nband_eff, iband_min, iband_max
334 : integer,parameter :: gpu_option=ABI_GPU_DISABLED
335 : type(xgBlock_t) :: eigenBlock !
336 : type(xgBlock_t) :: residuBlock,occBlock
337 : type(xg_t):: cprj_work_all
338 : double precision :: maxResidu, minResidu, dummy
339 : double precision :: dlamch,tolerance
340 : integer :: ierr = 0
341 : integer :: nrestart
342 : double precision :: tsec(2)
343 : character(len=500) :: msg
344 :
345 : interface
346 : subroutine getAX(X,AX)
347 : use m_xg, only : xgBlock_t
348 : type(xgBlock_t), intent(inout) :: X
349 : type(xgBlock_t), intent(inout) :: AX
350 : end subroutine getAX
351 : end interface
352 :
353 4036 : lobpcg%prtvol = prtvol
354 :
355 8072 : tolerance=2*dlamch('E')
356 :
357 4036 : blockdim = lobpcg%blockdim
358 4036 : blockdim2 = 2*blockdim
359 4036 : blockdim3 = 3*blockdim
360 4036 : blockdim_cprj = lobpcg%blockdim_cprj
361 4036 : spacedim = lobpcg%spacedim
362 :
363 4036 : nblock = lobpcg%nblock
364 4036 : nline = lobpcg%nline
365 :
366 4036 : if (nbdbuf>0) then
367 4036 : nband_eff = lobpcg%neigenpairs - nbdbuf
368 : else
369 0 : nband_eff = lobpcg%neigenpairs
370 : end if
371 :
372 4036 : call xgBlock_getSize(eigen,rows_tmp, cols_tmp)
373 4036 : if ( rows_tmp /= lobpcg%neigenpairs .and. cols_tmp /= 1 ) then
374 0 : ABI_ERROR("Error eigen size")
375 : endif
376 4036 : call xgBlock_getSize(X0,rows_tmp, cols_tmp)
377 4036 : if ( rows_tmp /= lobpcg%spacedim ) then
378 0 : ABI_ERROR("Error X0 spacedim")
379 : endif
380 4036 : if ( cols_tmp /= lobpcg%neigenpairs ) then
381 0 : ABI_ERROR("Error X0 npairs")
382 : endif
383 :
384 4036 : if (isppol==1.and.ikpt==1.and.inonsc==1.and.istep==1) then
385 42 : write(msg,'(a,es16.6)') ' lobpcg%tolerance(tolwfr_diago)=',lobpcg%tolerance
386 42 : call wrtout(std_out,msg,'COLL')
387 : end if
388 :
389 4036 : call xg_init(eigenvalues3N,SPACE_R,blockdim3,1)
390 4036 : call xg_setBlock(eigenvalues3N,eigenvaluesN,blockdim,1)
391 4036 : call xg_setBlock(eigenvalues3N,eigenvalues2N,blockdim2,1)
392 :
393 4036 : call xgBlock_reshape(eigen,blockdim,nblock)
394 4036 : call xgBlock_reshape(residu,blockdim,nblock)
395 4036 : call xgBlock_reshape(occ,blockdim,nblock)
396 :
397 4036 : lobpcg%AllX0 = X0
398 4036 : lobpcg%AllcprjX0 = cprjX0
399 :
400 4036 : call xg_init(residu_eff,SPACE_R,blockdim,1,gpu_option=ABI_GPU_DISABLED)
401 :
402 4036 : if ( lobpcg%paral_kgb == 1 ) then
403 2684 : call timab(tim_transpose,1,tsec)
404 : call xgTransposer_constructor(lobpcg%xgTransposerX,lobpcg%X,lobpcg%XColsRows,nspinor,&
405 2684 : STATE_LINALG,TRANS_ALL2ALL,lobpcg%comm_rows,lobpcg%comm_cols,0,0,lobpcg%me_g0_fft)
406 : call xgTransposer_copyConstructor(lobpcg%xgTransposerAX,lobpcg%xgTransposerX,&
407 2684 : lobpcg%AX,lobpcg%AXColsRows,STATE_LINALG)
408 : call xgTransposer_copyConstructor(lobpcg%xgTransposerW,lobpcg%xgTransposerX,&
409 2684 : lobpcg%W,lobpcg%WColsRows,STATE_LINALG)
410 : call xgTransposer_copyConstructor(lobpcg%xgTransposerAW,lobpcg%xgTransposerX,&
411 2684 : lobpcg%AW,lobpcg%AWColsRows,STATE_LINALG)
412 2684 : call timab(tim_transpose,2,tsec)
413 : else
414 1352 : call xgBlock_setBlock(lobpcg%X, lobpcg%XColsRows, spacedim, blockdim)
415 1352 : call xgBlock_setBlock(lobpcg%AX, lobpcg%AXColsRows, spacedim, blockdim)
416 1352 : call xgBlock_setBlock(lobpcg%W, lobpcg%WColsRows, spacedim, blockdim)
417 1352 : call xgBlock_setBlock(lobpcg%AW, lobpcg%AWColsRows, spacedim, blockdim)
418 : end if
419 :
420 : !! Start big loop over blocks
421 20628 : do iblock = 1, nblock
422 16592 : nrestart = 0
423 :
424 16592 : call lobpcg_getX0(lobpcg,iblock)
425 16592 : call xgBlock_setBlock(residu,residuBlock,blockdim,1,fcol=iblock)
426 16592 : call xgBlock_setBlock(occ, occBlock, blockdim,1,fcol=iblock)
427 :
428 16592 : call timab(tim_cprj,1,tsec)
429 16592 : call xg_nonlop_getcprj(lobpcg%xg_nonlop,lobpcg%X,lobpcg%cprjX,lobpcg%cprj_work%self)
430 16592 : call timab(tim_cprj,2,tsec)
431 :
432 16592 : if ( iblock > 1 ) then
433 12556 : call lobpcg_setPreviousX0(lobpcg,iblock)
434 : ! Orthogonalize current iblock X block With Respect To previous Blocks in B-basis
435 12556 : call lobpcg_orthoXwrtBlocks(lobpcg,lobpcg%X,lobpcg%cprjX,iblock,lobpcg%cprj_work%self)
436 : end if
437 :
438 16592 : if (lobpcg%paral_kgb == 1) then
439 11184 : call timab(tim_transpose,1,tsec)
440 11184 : call xgTransposer_transpose(lobpcg%xgTransposerX,STATE_COLSROWS)
441 11184 : lobpcg%xgTransposerAX%state=STATE_COLSROWS
442 11184 : call timab(tim_transpose,2,tsec)
443 : end if
444 :
445 : ! Initialize some quantitites (AX)
446 16592 : call timab(tim_ax_v,1,tsec)
447 16592 : call getAX(lobpcg%XColsRows,lobpcg%AXColsRows)
448 16592 : call xgBlock_zero_im_g0(lobpcg%AXColsRows)
449 16592 : call timab(tim_ax_v,2,tsec)
450 16592 : if (lobpcg%paral_kgb == 1) then
451 11184 : call timab(tim_transpose,1,tsec)
452 11184 : lobpcg%xgTransposerX%state=STATE_LINALG
453 11184 : call xgTransposer_transpose(lobpcg%xgTransposerAX,STATE_LINALG)
454 11184 : call timab(tim_transpose,2,tsec)
455 : end if
456 :
457 16592 : call timab(tim_ax_k,1,tsec)
458 16592 : call xgBlock_add_diag(lobpcg%X,kin,nspinor,lobpcg%AX)
459 16592 : call timab(tim_ax_k,2,tsec)
460 :
461 : ! B-orthonormalize X, BX and AX
462 16592 : call xg_Borthonormalize_cprj(lobpcg%xg_nonlop,lobpcg%X,lobpcg%cprjX,ierr,tim_Bortho_X,gpu_option,AX=lobpcg%AX)
463 :
464 : ! Do first RR on X to get the first eigen values
465 : call xg_RayleighRitz_cprj(lobpcg%xg_nonlop,lobpcg%X,lobpcg%cprjX,lobpcg%AX,eigenvaluesN,ierr,&
466 16592 : & lobpcg%prtvol,tim_RR_X,gpu_option,add_Anl=.True.)
467 :
468 16592 : compute_residu = .true.
469 :
470 75096 : do iline = 1, nline
471 :
472 67717 : if ( ierr /= 0 ) then
473 : !ABI_COMMENT("Consider using more bands and nbdbuf if necessary.")
474 0 : ierr = 0
475 : end if
476 :
477 : ! Compute AX-Lambda*BX
478 67717 : call timab(tim_copy,1,tsec)
479 67717 : call xgBlock_copy(lobpcg%AX,lobpcg%W)
480 67717 : call timab(tim_copy,2,tsec)
481 : ! Add the non-local part: (A_nl - Lambda*B)X
482 67717 : call timab(tim_ax_nl,1,tsec)
483 67717 : if (lobpcg%xg_nonlop%paw) then
484 : call xg_nonlop_getHmeSX(lobpcg%xg_nonlop,lobpcg%X,lobpcg%cprjX,lobpcg%W,eigenvaluesN,&
485 55920 : lobpcg%cprj_work%self,lobpcg%cprj_work2%self)
486 : else
487 : call xg_nonlop_getHX(lobpcg%xg_nonlop,lobpcg%W,lobpcg%cprjX,&
488 11797 : lobpcg%cprj_work%self,lobpcg%cprj_work2%self)
489 11797 : call xgBlock_yxmax(lobpcg%W,eigenvaluesN,lobpcg%X)
490 : end if
491 67717 : call timab(tim_ax_nl,2,tsec)
492 :
493 : ! Compute residu norm here !
494 67717 : call timab(tim_maxres,1,tsec)
495 67717 : call xgBlock_colwiseNorm2(lobpcg%W,residuBlock)
496 67717 : call timab(tim_maxres,2,tsec)
497 :
498 : ! Apply preconditioner
499 67717 : call timab(tim_pcond,1,tsec)
500 67717 : call xgBlock_apply_diag(lobpcg%W,pcond,nspinor)
501 67717 : call timab(tim_pcond,2,tsec)
502 :
503 67717 : call timab(tim_nbdbuf,1,tsec)
504 67717 : if (nbdbuf>=0) then
505 : ! There is a transfer from GPU to CPU in this copy
506 67717 : call xgBlock_copy(residuBlock,residu_eff%self)
507 67717 : iband_min = 1 + blockdim*(iblock-1)
508 67717 : iband_max = blockdim*iblock
509 67717 : if (iband_max<=nband_eff) then ! all bands of this block are below nband_eff
510 57269 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu)
511 10448 : else if (iband_min<=nband_eff) then ! some bands of this block are below nband_eff
512 3120 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu,row_bound=(nband_eff-iband_min+1))
513 : else ! all bands of this block are above nband_eff
514 7328 : minResidu = 0.0
515 7328 : maxResidu = 0.0
516 : end if
517 0 : else if (nbdbuf==-101) then
518 0 : call xgBlock_minmax(residuBlock,minResidu,dummy) ! Get minimum of true residuals
519 : ! Compute effective residuals : res_eff = res * occ
520 0 : call xgBlock_apply_diag(residuBlock,occBlock,1,Y=residu_eff%self)
521 0 : call xgBlock_minmax(residu_eff%self,dummy,maxResidu) ! Get maximum of effective residuals
522 : else
523 0 : ABI_ERROR('Bad value of nbdbuf')
524 : end if
525 67717 : call timab(tim_nbdbuf,2,tsec)
526 67717 : if ( maxResidu < lobpcg%tolerance ) then
527 : compute_residu = .false.
528 : exit
529 : end if
530 :
531 58504 : call timab(tim_cprj,1,tsec)
532 58504 : call xg_nonlop_getcprj(lobpcg%xg_nonlop,lobpcg%W,lobpcg%cprjW,lobpcg%cprj_work%self)
533 58504 : call timab(tim_cprj,2,tsec)
534 :
535 : ! Orthonormalize with respect to previous blocks
536 58504 : if ( iblock > 1 ) then
537 32398 : call lobpcg_orthoXwrtBlocks(lobpcg,lobpcg%W,lobpcg%cprjW,iblock,lobpcg%cprj_work%self)
538 : end if
539 :
540 58504 : if (lobpcg%paral_kgb == 1) then
541 37406 : call timab(tim_transpose,1,tsec)
542 37406 : call xgTransposer_transpose(lobpcg%xgTransposerW,STATE_COLSROWS)
543 37406 : lobpcg%xgTransposerAW%state=STATE_COLSROWS
544 37406 : call timab(tim_transpose,2,tsec)
545 : end if
546 : ! Apply A and B on W
547 58504 : call timab(tim_ax_v,1,tsec)
548 58504 : call getAX(lobpcg%WColsRows,lobpcg%AWColsRows)
549 58504 : call xgBlock_zero_im_g0(lobpcg%AWColsRows)
550 58504 : call timab(tim_ax_v,2,tsec)
551 58504 : if (lobpcg%paral_kgb == 1) then
552 37406 : call timab(tim_transpose,1,tsec)
553 37406 : lobpcg%xgTransposerW%state=STATE_LINALG
554 37406 : call xgTransposer_transpose(lobpcg%xgTransposerAW,STATE_LINALG)
555 37406 : call timab(tim_transpose,2,tsec)
556 : end if
557 :
558 58504 : call timab(tim_ax_k,1,tsec)
559 58504 : call xgBlock_add_diag(lobpcg%W,kin,nspinor,lobpcg%AW)
560 58504 : call timab(tim_ax_k,2,tsec)
561 :
562 : ! DO RR in the correct subspace
563 : ! if residu starts to be too small, there is an accumulation error in
564 : ! P with values such as 1e-29 that make the eigenvectors diverge
565 65883 : if ( iline == 1 .or. minResidu < 1e-27) then
566 : ! Do RR on XW to get the eigen vectors
567 : call xg_Borthonormalize_cprj(lobpcg%xg_nonlop,lobpcg%XW,lobpcg%cprjXW,&
568 8938 : & ierr,tim_Bortho_XW,gpu_option,AX=lobpcg%AXW,blockdim_cprj=blockdim_cprj)
569 8938 : call xgBlock_zero(lobpcg%P)
570 8938 : call xgBlock_zero(lobpcg%AP)
571 8938 : call xgBlock_zero(lobpcg%cprjP)
572 8938 : if ( ierr /= 0 ) then
573 0 : ABI_COMMENT("B-orthonormalization (XW) did not work.")
574 : end if
575 : call xg_RayleighRitz_cprj(lobpcg%xg_nonlop,lobpcg%X,lobpcg%cprjX,lobpcg%AX,eigenvaluesN,ierr,&
576 : & lobpcg%prtvol,tim_RR_XW,gpu_option,tolerance=tolerance,&
577 : & XW=lobpcg%XW,W=lobpcg%W,cprjXW=lobpcg%cprjXW,cprjW=lobpcg%cprjW,AW=lobpcg%AW,&
578 : & P=lobpcg%P,cprjP=lobpcg%cprjP,AP=lobpcg%AP,WP=lobpcg%WP,cprjWP=lobpcg%cprjWP,&
579 8938 : & AWP=lobpcg%AWP,blockdim_cprj=blockdim_cprj,add_Anl=.True.)
580 8938 : if ( ierr /= 0 ) then
581 0 : ABI_WARNING("RayleighRitz (XW) did not work, but continue anyway.")
582 : exit
583 : end if
584 : else
585 : ! B-orthonormalize P, BP
586 : call xg_Borthonormalize_cprj(lobpcg%xg_nonlop,lobpcg%XWP%self,lobpcg%cprjXWP%self,&
587 49566 : & ierr,tim_Bortho_XWP,gpu_option,AX=lobpcg%AXWP%self,blockdim_cprj=blockdim_cprj)
588 : ! Do RR on XWP to get the eigen vectors
589 49566 : if ( ierr == 0 ) then
590 : call xg_RayleighRitz_cprj(lobpcg%xg_nonlop,lobpcg%X,lobpcg%cprjX,lobpcg%AX,eigenvaluesN,ierr,&
591 : & lobpcg%prtvol,tim_RR_XWP,gpu_option,tolerance=tolerance,&
592 : & XW=lobpcg%XW,W=lobpcg%W,cprjXW=lobpcg%cprjXW,cprjW=lobpcg%cprjW,AW=lobpcg%AW,&
593 : & P=lobpcg%P,cprjP=lobpcg%cprjP,AP=lobpcg%AP,WP=lobpcg%WP,cprjWP=lobpcg%cprjWP,&
594 49566 : & AWP=lobpcg%AWP,XWP=lobpcg%XWP%self,cprjXWP=lobpcg%cprjXWP%self,blockdim_cprj=blockdim_cprj,add_Anl=.True.)
595 49566 : if ( ierr /= 0 ) then
596 0 : ABI_WARNING("RayleighRitz (XWP) did not work, but continue anyway.")
597 : exit
598 : end if
599 : else
600 0 : ABI_COMMENT("B-orthonormalization (XWP) did not work, try on XW.")
601 : call xg_Borthonormalize_cprj(lobpcg%xg_nonlop,lobpcg%XW,lobpcg%cprjXW,&
602 0 : & ierr,tim_Bortho_XW,gpu_option,AX=lobpcg%AXW,blockdim_cprj=blockdim_cprj)
603 0 : if ( ierr /= 0 ) then
604 0 : ABI_COMMENT("B-orthonormalization (XW) did not work.")
605 : end if
606 0 : call xgBlock_zero(lobpcg%P)
607 0 : call xgBlock_zero(lobpcg%AP)
608 0 : call xgBlock_zero(lobpcg%cprjP)
609 0 : nrestart = nrestart + 1
610 : call xg_RayleighRitz_cprj(lobpcg%xg_nonlop,lobpcg%X,lobpcg%cprjX,lobpcg%AX,eigenvaluesN,ierr,&
611 : & lobpcg%prtvol,tim_RR_XW,gpu_option,tolerance=tolerance,&
612 : & XW=lobpcg%XW,W=lobpcg%W,cprjXW=lobpcg%cprjXW,cprjW=lobpcg%cprjW,AW=lobpcg%AW,&
613 : & P=lobpcg%P,cprjP=lobpcg%cprjP,AP=lobpcg%AP,WP=lobpcg%WP,cprjWP=lobpcg%cprjWP,&
614 0 : & AWP=lobpcg%AWP,blockdim_cprj=blockdim_cprj,add_Anl=.True.)
615 0 : if ( ierr /= 0 ) then
616 0 : ABI_WARNING("RayleighRitz (XW) did not work, but continue anyway.")
617 : exit
618 : end if
619 : end if
620 : end if
621 :
622 : end do
623 :
624 16592 : if ( compute_residu ) then
625 : ! Recompute AX-Lambda*BX for the last time
626 7379 : call timab(tim_copy,1,tsec)
627 7379 : call xgBlock_copy(lobpcg%AX,lobpcg%W)
628 7379 : call timab(tim_copy,2,tsec)
629 7379 : call timab(tim_ax_nl,1,tsec)
630 7379 : if (lobpcg%xg_nonlop%paw) then
631 : call xg_nonlop_getHmeSX(lobpcg%xg_nonlop,lobpcg%X,lobpcg%cprjX,lobpcg%W,eigenvaluesN,&
632 6305 : lobpcg%cprj_work%self,lobpcg%cprj_work2%self)
633 : else
634 : call xg_nonlop_getHX(lobpcg%xg_nonlop,lobpcg%W,lobpcg%cprjX,&
635 1074 : lobpcg%cprj_work%self,lobpcg%cprj_work2%self)
636 1074 : call xgBlock_yxmax(lobpcg%W,eigenvaluesN,lobpcg%X)
637 : end if
638 7379 : call timab(tim_ax_nl,2,tsec)
639 : ! Recompute residu norm here !
640 7379 : call timab(tim_maxres,1,tsec)
641 7379 : call xgBlock_colwiseNorm2(lobpcg%W,residuBlock)
642 7379 : call timab(tim_maxres,2,tsec)
643 : ! Apply preconditioner
644 7379 : call timab(tim_pcond,1,tsec)
645 7379 : call xgBlock_apply_diag(lobpcg%W,pcond,nspinor)
646 7379 : call timab(tim_pcond,2,tsec)
647 :
648 7379 : call timab(tim_nbdbuf,1,tsec)
649 7379 : if (nbdbuf>=0) then
650 7379 : call xgBlock_copy(residuBlock,residu_eff%self)
651 7379 : iband_min = 1 + blockdim*(iblock-1)
652 7379 : iband_max = blockdim*iblock
653 7379 : if (iband_max<=nband_eff) then ! all bands of this block are below nband_eff
654 6871 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu)
655 508 : else if (iband_min<=nband_eff) then ! some bands of this block are below nband_eff
656 508 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu,row_bound=(nband_eff-iband_min+1))
657 : else ! all bands of this block are above nband_eff
658 0 : minResidu = 0.0
659 0 : maxResidu = 0.0
660 : end if
661 0 : else if (nbdbuf==-101) then
662 0 : call xgBlock_minmax(residuBlock,minResidu,dummy) ! Get minimum of true residuals
663 : ! Compute effective residuals : res_eff = res * occ
664 0 : call xgBlock_apply_diag(residuBlock,occBlock,1,Y=residu_eff%self)
665 0 : call xgBlock_minmax(residu_eff%self,dummy,maxResidu) ! Get maximum of effective residuals
666 : else
667 0 : ABI_ERROR('Bad value of nbdbuf')
668 : end if
669 7379 : call timab(tim_nbdbuf,2,tsec)
670 : end if
671 :
672 16592 : if (prtvol==5.and.xmpi_comm_rank(lobpcg%spacecom)==0) then
673 0 : write(msg,'(6(a,i4),2(a,es16.6))') 'lobpcg | istep=',istep,'| isppol=',isppol,'| ikpt=',ikpt,&
674 0 : & '| inonsc=',inonsc,'| iblock=',iblock,'| nline_done=',iline-1,'| minRes=',minResidu,'| maxRes=',maxResidu
675 0 : call wrtout(std_out,msg,'PERS')
676 : end if
677 :
678 : ! Save eigenvalues
679 16592 : call timab(tim_copy,1,tsec)
680 16592 : call xgBlock_setBlock(eigen,eigenBlock,blockdim,1,fcol=iblock)
681 16592 : call xgBlock_copy(eigenvaluesN,eigenBlock)
682 16592 : call timab(tim_copy,2,tsec)
683 :
684 : ! Save new X in X0
685 16592 : call lobpcg_setX0(lobpcg,iblock)
686 :
687 : ! Copy previous BX into BX0 for previous block
688 20628 : if ( nblock > 1 ) then
689 16560 : call lobpcg_transferAX_BX(lobpcg,iblock)
690 : end if
691 :
692 : end do !! End iblock loop
693 :
694 4036 : call xgBlock_reshape(eigen,blockdim*nblock,1)
695 4036 : call xgBlock_reshape(residu,blockdim*nblock,1)
696 4036 : call xgBlock_reshape(occ,blockdim*nblock,1)
697 :
698 4036 : call xg_free(eigenvalues3N)
699 4036 : call xg_free(residu_eff)
700 :
701 4036 : skip = .false.
702 4036 : if ( ierr /= 0 ) then
703 0 : ABI_COMMENT("Some errors happened, so H|Psi> and S|Psi> are computed before leaving")
704 0 : if ( lobpcg%paral_kgb == 1 ) then
705 : call xgTransposer_constructor(lobpcg%xgTransposerAllX0,lobpcg%AllX0,lobpcg%AllX0ColsRows,nspinor,&
706 0 : STATE_LINALG,TRANS_ALL2ALL,lobpcg%comm_rows,lobpcg%comm_cols,0,0,lobpcg%me_g0_fft)
707 : call xgTransposer_copyConstructor(lobpcg%xgTransposerAllAX0,lobpcg%xgTransposerAllX0,&
708 0 : lobpcg%AllAX0%self,lobpcg%AllAX0ColsRows,STATE_LINALG)
709 : else
710 0 : call xgBlock_setBlock(lobpcg%AllX0 , lobpcg%AllX0ColsRows , spacedim, lobpcg%neigenpairs)
711 0 : call xgBlock_setBlock(lobpcg%AllAX0%self, lobpcg%AllAX0ColsRows, spacedim, lobpcg%neigenpairs)
712 : end if
713 0 : if (lobpcg%paral_kgb == 1) then
714 0 : call timab(tim_transpose,1,tsec)
715 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllX0,STATE_COLSROWS)
716 0 : lobpcg%xgTransposerAllAX0%state=STATE_COLSROWS
717 0 : call timab(tim_transpose,2,tsec)
718 : end if
719 0 : call timab(tim_ax_v,1,tsec)
720 0 : call getAX(lobpcg%AllX0ColsRows,lobpcg%AllAX0ColsRows)
721 0 : call xgBlock_zero_im_g0(lobpcg%AllAX0ColsRows)
722 0 : call timab(tim_ax_v,2,tsec)
723 0 : if (lobpcg%paral_kgb == 1) then
724 0 : call timab(tim_transpose,1,tsec)
725 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllX0,STATE_LINALG)
726 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllAX0,STATE_LINALG)
727 0 : call timab(tim_transpose,2,tsec)
728 : end if
729 :
730 0 : call timab(tim_ax_k,1,tsec)
731 0 : call xgBlock_add_diag(lobpcg%AllX0,kin,nspinor,lobpcg%AllAX0%self)
732 0 : call timab(tim_ax_k,2,tsec)
733 :
734 0 : call xgTransposer_free(lobpcg%xgTransposerAllX0)
735 0 : call xgTransposer_free(lobpcg%xgTransposerAllAX0)
736 : skip = .true.
737 : end if
738 :
739 : if (.not.skip) then
740 :
741 4036 : if ( nblock > 1 ) then
742 4004 : call timab(tim_cprj,1,tsec)
743 4004 : call xg_init(cprj_work_all,space(cprjX0),rows(cprjX0),cols(cprjX0),comm(cprjX0))
744 : ! cprj computation could be avoided, with MPI comm instead
745 4004 : call xg_nonlop_getcprj(lobpcg%xg_nonlop,X0,cprjX0,cprj_work_all%self)
746 4004 : call xg_free(cprj_work_all)
747 4004 : call timab(tim_cprj,2,tsec)
748 : call xg_Borthonormalize_cprj(lobpcg%xg_nonlop,X0,cprjX0,&
749 4004 : & ierr,tim_Bortho_Xall,gpu_option,AX=lobpcg%AllAX0%self)
750 : call xg_RayleighRitz_cprj(lobpcg%xg_nonlop,X0,cprjX0,lobpcg%AllAX0%self,eigen,ierr,&
751 4004 : & lobpcg%prtvol,tim_RR_Xall,gpu_option,add_Anl=.True.)
752 : end if
753 :
754 : end if
755 :
756 4036 : if ( lobpcg%paral_kgb == 1 ) then
757 2684 : call xgTransposer_free(lobpcg%xgTransposerX)
758 2684 : call xgTransposer_free(lobpcg%xgTransposerAX)
759 2684 : call xgTransposer_free(lobpcg%xgTransposerW)
760 2684 : call xgTransposer_free(lobpcg%xgTransposerAW)
761 : end if
762 :
763 4036 : if (.not.lobpcg%xg_nonlop%paw) then
764 816 : call timab(tim_enl,1,tsec)
765 816 : call xg_init(cprj_work_all,space(cprjX0),rows(cprjX0),cols(cprjX0),comm(cprjX0))
766 816 : call xg_nonlop_colwiseXHX(lobpcg%xg_nonlop,cprjX0,cprj_work_all%self,enl)
767 816 : call xg_free(cprj_work_all)
768 816 : call timab(tim_enl,2,tsec)
769 : end if
770 :
771 4036 : end subroutine lobpcg_run_cprj
772 :
773 16592 : subroutine lobpcg_getX0(lobpcg,iblock)
774 :
775 : type(lobpcg_t) , intent(inout) :: lobpcg
776 : integer , intent(in ) :: iblock
777 : type(xgBlock_t) :: Xtmp,cprjXtmp
778 : integer :: blockdim
779 : integer :: spacedim
780 : integer :: cprjdim
781 : integer :: blockdim_cprj
782 : double precision :: tsec(2)
783 :
784 16592 : call timab(tim_copy,1,tsec)
785 :
786 16592 : blockdim = lobpcg%blockdim
787 16592 : spacedim = lobpcg%spacedim
788 :
789 16592 : cprjdim = lobpcg%cprjdim
790 16592 : blockdim_cprj = lobpcg%blockdim_cprj
791 :
792 16592 : call xgBlock_setBlock(lobpcg%AllX0,Xtmp,spacedim,blockdim,fcol=(iblock-1)*blockdim+1)
793 16592 : call xgBlock_copy(Xtmp,lobpcg%X)
794 :
795 16592 : call xgBlock_setBlock(lobpcg%AllcprjX0,cprjXtmp,cprjdim,blockdim_cprj,fcol=(iblock-1)*blockdim_cprj+1)
796 16592 : call xgBlock_copy(cprjXtmp,lobpcg%cprjX)
797 :
798 16592 : call timab(tim_copy,2,tsec)
799 :
800 16592 : end subroutine lobpcg_getX0
801 :
802 :
803 12556 : subroutine lobpcg_setPreviousX0(lobpcg,iblock)
804 :
805 : type(lobpcg_t) , intent(inout) :: lobpcg
806 : integer , intent(in ) :: iblock
807 :
808 12556 : if (iblock<2) then
809 0 : ABI_ERROR("iblock<2")
810 : end if
811 12556 : call xgBlock_setBlock(lobpcg%AllX0 ,lobpcg%X0 ,lobpcg%spacedim,(iblock-1)*lobpcg%blockdim)
812 12556 : call xgBlock_setBlock(lobpcg%AllcprjX0,lobpcg%cprjX0,lobpcg%cprjdim ,(iblock-1)*lobpcg%blockdim_cprj)
813 12556 : end subroutine lobpcg_setPreviousX0
814 :
815 44954 : subroutine lobpcg_orthoXwrtBlocks(lobpcg,var,cprjvar,iblock,cprj_work)
816 :
817 : type(lobpcg_t) , intent(inout) :: lobpcg
818 : type(xgBlock_t), intent(inout) :: var
819 : type(xgBlock_t), intent(inout) :: cprjvar
820 : type(xgBlock_t), intent(inout) :: cprj_work
821 : integer , intent(in ) :: iblock
822 : integer :: previousBlock
823 : integer :: blockdim
824 : integer :: spacedim
825 : integer :: space_buf
826 : integer :: nspinor,cprjdim
827 : type(xg_t) :: buffer
828 : double precision :: tsec(2)
829 : type(xgBlock_t) :: cprjX0_spinor,cprj_work_spinor
830 :
831 44954 : call timab(tim_ortho,1,tsec)
832 :
833 44954 : if (iblock<2) then
834 0 : ABI_ERROR("iblock<2")
835 : end if
836 :
837 44954 : if (cols(cprjvar)/=cols(cprj_work)) then
838 0 : ABI_ERROR("cprjvar and cprj_work should have same number of columns")
839 : end if
840 44954 : blockdim = lobpcg%blockdim
841 44954 : spacedim = lobpcg%spacedim
842 44954 : previousBlock = (iblock-1)*lobpcg%blockdim
843 :
844 44954 : cprjdim = lobpcg%xg_nonlop%cprjdim
845 44954 : nspinor = lobpcg%xg_nonlop%nspinor
846 :
847 44954 : space_buf = space(var)
848 44954 : if (space(var)==SPACE_CR) then
849 20858 : space_buf = SPACE_R
850 : end if
851 44954 : call xg_init(buffer,space_buf,previousBlock,blockdim,lobpcg%spacecom)
852 :
853 : ! buffer = X0^T*X
854 44954 : call xgBlock_gemm('t','n',1.0d0,lobpcg%X0,var,0.d0,buffer%self,comm=lobpcg%spacecom)
855 :
856 : ! Add the nonlocal part if paw
857 44954 : if (lobpcg%xg_nonlop%paw) then
858 37654 : call xg_nonlop_getXSY(lobpcg%xg_nonlop,lobpcg%cprjX0,cprjvar,cprj_work,buffer%self)
859 : end if
860 :
861 : ! sum all process contribution
862 : ! X = - X0*(BX0^T*X) + X
863 44954 : call xgBlock_gemm('n','n',-1.0d0,lobpcg%X0,buffer%self,1.0d0,var)
864 :
865 44954 : call xgBlock_zero(cprj_work)
866 44954 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,nspinor,COLS2ROWS)
867 44954 : call xgBlock_reshape_spinor(lobpcg%cprjX0,cprjX0_spinor,nspinor,COLS2ROWS)
868 : call xgBlock_gemm_mpi_cyclic_permutation(cprjX0_spinor,buffer%self,cprj_work_spinor,&
869 44954 : & lobpcg%xg_nonlop%me_band,blocksize=lobpcg%blockdim_cprj/nspinor,comm=lobpcg%xg_nonlop%comm_band)
870 44954 : call xgBlock_saxpy(cprjvar,-1.0d0,cprj_work)
871 :
872 44954 : call xg_free(buffer)
873 :
874 44954 : call timab(tim_ortho,2,tsec)
875 :
876 44954 : end subroutine lobpcg_orthoXwrtBlocks
877 :
878 16592 : subroutine lobpcg_setX0(lobpcg,iblock)
879 :
880 : type(lobpcg_t) , intent(inout) :: lobpcg
881 : integer , intent(in ) :: iblock
882 : type(xgBlock_t) :: Xtmp,cprjXtmp
883 : integer :: blockdim
884 : integer :: spacedim
885 : integer :: cprjdim
886 : integer :: blockdim_cprj
887 : double precision :: tsec(2)
888 :
889 16592 : call timab(tim_copy,1,tsec)
890 16592 : blockdim = lobpcg%blockdim
891 16592 : spacedim = lobpcg%spacedim
892 :
893 16592 : cprjdim = lobpcg%cprjdim
894 16592 : blockdim_cprj = lobpcg%blockdim_cprj
895 :
896 : !X0(:,(iblock-1)*blockdim+1:iblock*blockdim) = lobpcg%XWP(:,lobpcg%X+1:lobpcg%X+blockdim)
897 16592 : call xgBlock_setBlock(lobpcg%AllX0,Xtmp,spacedim,blockdim,fcol=(iblock-1)*blockdim+1)
898 16592 : call xgBlock_copy(lobpcg%X,Xtmp)
899 :
900 16592 : call xgBlock_setBlock(lobpcg%AllcprjX0,cprjXtmp,cprjdim,blockdim_cprj,fcol=(iblock-1)*blockdim_cprj+1)
901 16592 : call xgBlock_copy(lobpcg%cprjX,cprjXtmp)
902 16592 : call timab(tim_copy,2,tsec)
903 :
904 16592 : end subroutine lobpcg_setX0
905 :
906 16560 : subroutine lobpcg_transferAX_BX(lobpcg,iblock)
907 :
908 : type(lobpcg_t), intent(inout) :: lobpcg
909 : integer , intent(in ) :: iblock
910 : type(xgBlock_t) :: CXtmp
911 : integer :: firstcol
912 : double precision :: tsec(2)
913 :
914 16560 : call timab(tim_copy,1,tsec)
915 :
916 16560 : if (iblock<1) then
917 0 : ABI_ERROR("iblock<1")
918 : end if
919 :
920 : ! iblock goes from 1 to nblock-1 included
921 16560 : firstcol = (iblock-1)*lobpcg%blockdim+1 ! Start of each block
922 :
923 : ! AX
924 16560 : call xg_setBlock(lobpcg%AllAX0,CXtmp,lobpcg%spacedim,lobpcg%blockdim,fcol=firstcol)
925 16560 : call xgBlock_copy(lobpcg%AX,CXtmp)
926 :
927 16560 : call timab(tim_copy,2,tsec)
928 :
929 16560 : end subroutine lobpcg_transferAX_BX
930 :
931 8072 : subroutine lobpcg_free(lobpcg)
932 :
933 : type(lobpcg_t), intent(inout) :: lobpcg
934 :
935 8072 : call xg_free(lobpcg%XWP)
936 8072 : call xg_free(lobpcg%AXWP)
937 8072 : call xg_free(lobpcg%cprjXWP)
938 8072 : call xg_free(lobpcg%cprj_work)
939 8072 : call xg_free(lobpcg%cprj_work2)
940 8072 : call xg_free(lobpcg%AllAX0)
941 :
942 8072 : end subroutine lobpcg_free
943 :
944 0 : end module m_lobpcg2_cprj
945 : !!***
|
