Branch data Line data Source code
1 : : !!****f* ABINIT/m_lobpcg2
2 : : !! NAME
3 : : !! m_lobpcg2
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the types and routines used to apply the
7 : : !! LOBPCG method (second version introduced by J. Bieder), using the xg_tools.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2015-2026 ABINIT group (J. Bieder, L. Baguet)
11 : : !! This file is distributed under the terms of the
12 : : !! gnu general public license, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !! for the initials of contributors, see ~abinit/doc/developers/contributors.txt .
15 : : !!
16 : : !! SOURCE
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : :
22 : : #include "abi_common.h"
23 : :
24 : : ! nvtx related macro definition
25 : : #include "nvtx_macros.h"
26 : :
27 : : module m_lobpcg2
28 : :
29 : : use m_xg
30 : : use m_xgTransposer
31 : : use m_xg_ortho_RR
32 : : use defs_basis
33 : : use m_abicore
34 : : use m_errors
35 : : use m_xomp
36 : : #ifdef HAVE_OPENMP
37 : : use omp_lib
38 : : #endif
39 : : use m_xmpi
40 : : use, intrinsic :: iso_c_binding, only: c_size_t
41 : :
42 : : #if defined(HAVE_GPU_MARKERS)
43 : : use m_nvtx_data
44 : : #endif
45 : :
46 : : use m_time, only : timab
47 : :
48 : : implicit none
49 : :
50 : : private
51 : :
52 : : integer, parameter :: VAR_X = 1000
53 : : integer, parameter :: VAR_W = 1001
54 : : integer, parameter :: VAR_P = 1002
55 : : integer, parameter :: VAR_XW = 1010
56 : : integer, parameter :: VAR_WP = 1011
57 : : integer, parameter :: VAR_XWP = 1100
58 : :
59 : : integer, parameter :: tim_init = 1651
60 : : integer, parameter :: tim_free = 1652
61 : : integer, parameter :: tim_copy = 1653
62 : : integer, parameter :: tim_getAX_BX = 1654
63 : : integer, parameter :: tim_ortho = 1655
64 : : integer, parameter :: tim_nbdbuf = 1656
65 : : ! integer, parameter :: tim_RR = 1657
66 : : integer, parameter :: tim_maxres = 1658
67 : : integer, parameter :: tim_ax_bx = 1659
68 : : integer, parameter :: tim_pcond = 1660
69 : : ! integer, parameter :: tim_hegv = 1661
70 : :
71 : : integer, parameter :: tim_Bortho_X = 1641
72 : : integer, parameter :: tim_Bortho_XW = 1642
73 : : integer, parameter :: tim_Bortho_XWP = 1643
74 : : integer, parameter :: tim_Bortho_Xall = 1644
75 : : integer, parameter :: tim_RR_X = 1645
76 : : integer, parameter :: tim_RR_XW = 1646
77 : : integer, parameter :: tim_RR_XWP = 1647
78 : : integer, parameter :: tim_RR_Xall = 1648
79 : :
80 : : integer, parameter :: tim_transpose = 1649
81 : :
82 : : type, public :: lobpcg_t
83 : : logical :: is_nested ! For OpenMP nested region
84 : : integer :: spacedim ! Space dimension for one vector
85 : : integer :: neigenpairs ! Number of eigen values/vectors we want
86 : : integer :: nblock ! Number of block in the space of dim spacedim
87 : : integer :: blockdim ! Number of vectors in one block
88 : : integer :: nline ! Number of line to perform
89 : : integer :: spacecom ! Communicator for MPI
90 : : integer :: paral_kgb ! paral_kgb formalism or not
91 : : integer :: comm_rows ! communicator for rows
92 : : integer :: comm_cols ! communicator for cols
93 : : integer :: me_g0 ! =1 if the processor have G=0 (linalg representation)
94 : : integer :: me_g0_fft ! =1 if the processor have G=0 (fft_representation)
95 : : integer :: gpu_option ! Which GPU version is used (0=none)
96 : : integer :: gpu_thread_limit ! When GPU is enabled, how many CPU threads to use in sensitive areas
97 : : double precision :: tolerance ! Tolerance on the residu to stop the minimization
98 : : integer :: prtvol
99 : : type(xgBlock_t) :: AllX0 ! Block of initial and final solution.
100 : : type(xgBlock_t) :: X0 ! Block of initial and final solution.
101 : : ! Dim is (cplx*spacedim,neigenpair)
102 : : !double precision, allocatable :: eig(:) ! Actual eigen values
103 : :
104 : : ! Some denomination here
105 : : ! X, W, P are the shifts where their block vectors starts.
106 : : ! X is the current block we are solving
107 : : ! W is the residu of the pencil A-lambda*B ie W=A*X-lambdaB*X
108 : : ! P correspond to X_{iline} - X_{iline-1}
109 : : ! if prefixed by A (AX, AW, AP) it means the result of A*X (A*W, A*P ! respectively)
110 : : ! Idem for B (BX, BW, BP) which is the result of B*X (B*W, B*P)
111 : : ! An extend to this is AXWP which is the vector [AX,AW,AP]
112 : : ! and BXWP which is [BX,BW,BP] columnwise
113 : : ! Now we only allocate the triple XWP matrix and each individual X, W and P
114 : : ! 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
115 : : ! P to p1-1
116 : : type(xg_t) :: XWP
117 : : type(xg_t) :: AXWP
118 : : type(xg_t) :: BXWP
119 : :
120 : : type(xgBlock_t) :: XColsRows
121 : : type(xgBlock_t) :: AXColsRows
122 : : type(xgBlock_t) :: BXColsRows
123 : :
124 : : type(xgTransposer_t) :: xgTransposerX
125 : : type(xgTransposer_t) :: xgTransposerAX
126 : : type(xgTransposer_t) :: xgTransposerBX
127 : :
128 : : type(xgBlock_t) :: WColsRows
129 : : type(xgBlock_t) :: AWColsRows
130 : : type(xgBlock_t) :: BWColsRows
131 : :
132 : : type(xgTransposer_t) :: xgTransposerW
133 : : type(xgTransposer_t) :: xgTransposerAW
134 : : type(xgTransposer_t) :: xgTransposerBW
135 : :
136 : : type(xgBlock_t) :: X ! Shift to apply to start reading the X values
137 : : type(xgBlock_t) :: W ! Shift to apply to start reading the W values
138 : : type(xgBlock_t) :: P ! Shift to apply to start reading the P values
139 : : type(xgBlock_t) :: XW ! Shift to apply to start reading the P values
140 : : type(xgBlock_t) :: WP ! Shift to apply to start reading the P values
141 : :
142 : : type(xgBlock_t) :: AX ! Shift to apply to start reading the X values
143 : : type(xgBlock_t) :: AW ! Shift to apply to start reading the W values
144 : : type(xgBlock_t) :: AP ! Shift to apply to start reading the P values
145 : : type(xgBlock_t) :: AXW ! Shift to apply to start reading the P values
146 : : type(xgBlock_t) :: AWP ! Shift to apply to start reading the P values
147 : :
148 : : type(xgBlock_t) :: BX ! Shift to apply to start reading the X values
149 : : type(xgBlock_t) :: BW ! Shift to apply to start reading the W values
150 : : type(xgBlock_t) :: BP ! Shift to apply to start reading the P values
151 : : type(xgBlock_t) :: BXW ! Shift to apply to start reading the P values
152 : : type(xgBlock_t) :: BWP ! Shift to apply to start reading the P values
153 : :
154 : : type(xg_t) :: AllBX0
155 : : type(xg_t) :: AllAX0
156 : : type(xgBlock_t) :: BX0
157 : :
158 : : type(xgBlock_t) :: AllX0ColsRows
159 : : type(xgBlock_t) :: AllAX0ColsRows
160 : : type(xgBlock_t) :: AllBX0ColsRows
161 : :
162 : : type(xgTransposer_t) :: xgTransposerAllX0
163 : : type(xgTransposer_t) :: xgTransposerAllAX0
164 : : type(xgTransposer_t) :: xgTransposerAllBX0
165 : :
166 : : ! Variable for work for lapack
167 : : end type lobpcg_t
168 : :
169 : : public :: lobpcg_init
170 : : public :: lobpcg_memInfo
171 : : public :: lobpcg_run
172 : : public :: lobpcg_free
173 : :
174 : : contains
175 : :
176 : :
177 : 24171 : subroutine lobpcg_init(lobpcg, neigenpairs, spacedim, blockdim, tolerance, nline, &
178 : : & space, spacecom, paral_kgb, comm_rows, comm_cols, me_g0, me_g0_fft, gpu_option, &
179 : : & gpu_thread_limit)
180 : :
181 : : type(lobpcg_t) , intent(inout) :: lobpcg
182 : : integer , intent(in ) :: neigenpairs
183 : : integer , intent(in ) :: spacedim
184 : : integer , intent(in ) :: blockdim
185 : : double precision, intent(in ) :: tolerance
186 : : integer , intent(in ) :: nline
187 : : integer , intent(in ) :: comm_rows,comm_cols
188 : : integer , intent(in ) :: space
189 : : integer , intent(in ) :: spacecom
190 : : integer , intent(in ) :: paral_kgb
191 : : integer , intent(in ) :: me_g0
192 : : integer , intent(in ) :: me_g0_fft
193 : : integer , intent(in ) :: gpu_option
194 : : integer,optional, intent(in ) :: gpu_thread_limit
195 : : double precision :: tsec(2)
196 : :
197 : 24171 : call timab(tim_init,1,tsec)
198 : 24171 : lobpcg%neigenpairs = neigenpairs
199 : 24171 : lobpcg%spacedim = spacedim
200 : 24171 : lobpcg%blockdim = blockdim
201 [ + + ]: 24171 : if (tolerance > 0.0) then
202 : 4568 : lobpcg%tolerance = tolerance
203 : : else
204 : 19603 : lobpcg%tolerance = 1.0e-20
205 : : end if
206 : 24171 : lobpcg%nline = nline
207 : 24171 : lobpcg%spacecom = spacecom
208 : 24171 : lobpcg%nblock = neigenpairs / blockdim
209 : 24171 : lobpcg%paral_kgb = paral_kgb
210 : 24171 : lobpcg%comm_rows = comm_rows
211 : 24171 : lobpcg%comm_cols = comm_cols
212 : 24171 : lobpcg%me_g0 = me_g0
213 : 24171 : lobpcg%me_g0_fft = me_g0_fft
214 : 24171 : lobpcg%gpu_option = gpu_option
215 : 24171 : lobpcg%gpu_thread_limit = 0
216 [ + - ]: 24171 : if(present(gpu_thread_limit)) lobpcg%gpu_thread_limit = gpu_thread_limit
217 : :
218 : 24171 : call lobpcg_allocateAll(lobpcg,space,me_g0)
219 : 24171 : call timab(tim_init,2,tsec)
220 : :
221 : 24171 : end subroutine lobpcg_init
222 : :
223 : :
224 : 24171 : subroutine lobpcg_allocateAll(lobpcg,space,me_g0)
225 : :
226 : : type(lobpcg_t) , intent(inout) :: lobpcg
227 : : integer , intent(in ) :: space
228 : : integer , intent(in ) :: me_g0
229 : : integer :: spacedim
230 : : integer :: blockdim
231 : :
232 : 24171 : spacedim = lobpcg%spacedim
233 : 24171 : blockdim = lobpcg%blockdim
234 : :
235 : 24171 : call lobpcg_free(lobpcg) ! Make sure everything is not allocated and
236 : : ! pointer point to null()
237 : :
238 : 24171 : call xg_init(lobpcg%XWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
239 : 24171 : call xg_setBlock(lobpcg%XWP,lobpcg%X,spacedim,blockdim)
240 : 24171 : call xg_setBlock(lobpcg%XWP,lobpcg%W,spacedim,blockdim,fcol=blockdim+1)
241 : 24171 : call xg_setBlock(lobpcg%XWP,lobpcg%P,spacedim,blockdim,fcol=2*blockdim+1)
242 : 24171 : call xg_setBlock(lobpcg%XWP,lobpcg%XW,spacedim,2*blockdim)
243 : 24171 : call xg_setBlock(lobpcg%XWP,lobpcg%WP,spacedim,2*blockdim,fcol=blockdim+1)
244 : :
245 : 24171 : call xg_init(lobpcg%AXWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
246 : 24171 : call xg_setBlock(lobpcg%AXWP,lobpcg%AX,spacedim,blockdim)
247 : 24171 : call xg_setBlock(lobpcg%AXWP,lobpcg%AW,spacedim,blockdim,fcol=blockdim+1)
248 : 24171 : call xg_setBlock(lobpcg%AXWP,lobpcg%AP,spacedim,blockdim,fcol=2*blockdim+1)
249 : 24171 : call xg_setBlock(lobpcg%AXWP,lobpcg%AXW,spacedim,2*blockdim)
250 : 24171 : call xg_setBlock(lobpcg%AXWP,lobpcg%AWP,spacedim,2*blockdim,fcol=blockdim+1)
251 : :
252 : 24171 : call xg_init(lobpcg%BXWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
253 : 24171 : call xg_setBlock(lobpcg%BXWP,lobpcg%BX,spacedim,blockdim)
254 : 24171 : call xg_setBlock(lobpcg%BXWP,lobpcg%BW,spacedim,blockdim,fcol=blockdim+1)
255 : 24171 : call xg_setBlock(lobpcg%BXWP,lobpcg%BP,spacedim,blockdim,fcol=2*blockdim+1)
256 : 24171 : call xg_setBlock(lobpcg%BXWP,lobpcg%BXW,spacedim,2*blockdim)
257 : 24171 : call xg_setBlock(lobpcg%BXWP,lobpcg%BWP,spacedim,2*blockdim,fcol=blockdim+1)
258 : :
259 [ + + ]: 24171 : if ( lobpcg%nblock /= 1 ) then
260 : 17811 : call xg_init(lobpcg%AllBX0,space,spacedim,lobpcg%neigenpairs,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
261 : 17811 : call xg_init(lobpcg%AllAX0,space,spacedim,lobpcg%neigenpairs,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
262 : : else
263 : 6360 : lobpcg%AllBX0%self = lobpcg%BX
264 : 6360 : lobpcg%AllAX0%self = lobpcg%AX
265 : : end if
266 : :
267 : :
268 : 24171 : end subroutine lobpcg_allocateAll
269 : :
270 : :
271 [ # # ]: 0 : function lobpcg_memInfo(neigenpairs, spacedim, space, paral_kgb, blockdim) result(arraymem)
272 : :
273 : : integer , intent(in ) :: neigenpairs
274 : : integer , intent(in ) :: spacedim
275 : : integer , intent(in ) :: blockdim
276 : : integer , intent(in ) :: space
277 : : integer , intent(in ) :: paral_kgb
278 : : integer(kind=c_size_t) :: memXWP
279 : : integer(kind=c_size_t) :: memAXWP
280 : : integer(kind=c_size_t) :: memBXWP
281 : : integer(kind=c_size_t) :: mem_xgTransposer
282 : : integer(kind=c_size_t) :: memAllBX0
283 : : integer(kind=c_size_t) :: memAllAX0
284 : : integer(kind=c_size_t) :: memeigenvalues3N
285 : : integer(kind=c_size_t) :: membufferOrtho
286 : : integer(kind=c_size_t) :: membufferBOrtho
287 : : integer(kind=c_size_t) :: memsubA
288 : : integer(kind=c_size_t) :: memsubB
289 : : integer(kind=c_size_t) :: memsubBtmp
290 : : integer(kind=c_size_t) :: memvec
291 : : integer(kind=c_size_t) :: memRR
292 : : integer(kind=c_size_t) :: maxmemTmp
293 : : integer(kind=c_size_t) :: cplx
294 : : integer :: nblock
295 : : integer(kind=c_size_t) :: arraymem(2)
296 : :
297 [ # # ]: 0 : cplx = 1 ; if ( space == SPACE_C ) cplx = 2
298 : 0 : nblock = neigenpairs/blockdim
299 : :
300 : : ! Permanent in lobpcg
301 : 0 : memXWP = int(cplx,c_size_t)* kind(1.d0) * spacedim * 3*blockdim
302 : 0 : memAXWP = int(cplx,c_size_t)* kind(1.d0) * spacedim * 3*blockdim
303 : 0 : memBXWP = int(cplx,c_size_t)* kind(1.d0) * spacedim * 3*blockdim
304 [ # # ]: 0 : if(paral_kgb > 1) then
305 : : ! xgtransposer *_ColRows buffers for X, AX, BX, W, AW, BW + internal send/recvbuf
306 : 0 : mem_xgTransposer = int(cplx,c_size_t)* kind(1.d0) * spacedim * 7*blockdim
307 : : else
308 : : mem_xgTransposer = 0
309 : : end if
310 [ # # ]: 0 : if ( nblock > 1 ) then
311 : 0 : memAllAX0 = int(cplx,c_size_t) * kind(1.d0) * spacedim * 3*blockdim
312 : 0 : memAllBX0 = int(cplx,c_size_t) * kind(1.d0) * spacedim * 3*blockdim
313 : 0 : membufferOrtho = int(cplx,c_size_t) * kind(1.d0) * blockdim * (nblock-1) * blockdim
314 : 0 : mem_xgTransposer = mem_xgTransposer + int(cplx,c_size_t) * kind(1.d0) * spacedim * 3*blockdim
315 : : else
316 : : memAllAX0 = 0
317 : : memAllBX0 = 0
318 : : membufferOrtho = 0
319 : : endif
320 : 0 : memeigenvalues3N = kind(1.d0) * 3*blockdim
321 : :
322 : : ! Temporary arrays
323 : :
324 : : ! For the moment being, only Bortho with X or WP at the same time
325 : 0 : membufferBOrtho = int(cplx,c_size_t) * kind(1.d0) * 2*blockdim * 2*blockdim
326 : 0 : memsubA = int(cplx,c_size_t) * kind(1.d0) * 3*blockdim * 3*blockdim
327 : 0 : memsubB = int(cplx,c_size_t) * kind(1.d0) * 3*blockdim * 3*blockdim
328 : 0 : memsubBtmp = int(cplx,c_size_t) * kind(1.d0) * spacedim * blockdim
329 : 0 : memvec = int(cplx,c_size_t) * kind(1.d0) * 3*blockdim * blockdim
330 : 0 : memRR = max(memsubA+memsubB+memvec,memsubBtmp+memvec)
331 : :
332 : 0 : maxmemTmp = max( membufferBOrtho,memRR,membufferOrtho )
333 : :
334 : 0 : arraymem(1) = memXWP+memAXWP+memBXWP+memAllAX0+memAllBX0+memeigenvalues3N+mem_xgTransposer
335 : 0 : arraymem(2) = maxmemTmp
336 : :
337 : 0 : end function lobpcg_memInfo
338 : :
339 : :
340 : 72513 : subroutine lobpcg_run(lobpcg, X0, getAX_BX, pcond, eigen, occ, residu, prtvol, nspinor, isppol, ikpt, inonsc, istep, nbdbuf)
341 : :
342 : : type(lobpcg_t) , intent(inout) :: lobpcg
343 : : type(xgBlock_t), intent(inout) :: X0 ! Full initial vectors
344 : : type(xgBlock_t), intent(inout) :: eigen ! Full initial eigen values
345 : : type(xgBlock_t), intent(inout) :: occ
346 : : type(xgBlock_t), intent(inout) :: residu
347 : : type(xgBlock_t), intent(in) :: pcond
348 : : integer , intent(in ) :: prtvol
349 : : integer , intent(in ) :: nspinor
350 : : integer , intent(in ) :: isppol,ikpt,inonsc,istep,nbdbuf
351 : :
352 : : type(xg_t) :: eigenvalues3N ! eigen values for Rayleight-Ritz
353 : : type(xg_t) :: residu_eff
354 : : type(xgBlock_t) :: eigenvaluesN ! eigen values for Rayleight-Ritz
355 : : type(xgBlock_t) :: eigenvalues2N ! eigen values for Rayleight-Ritz
356 : : logical :: skip,compute_residu
357 : : integer :: blockdim, blockdim3, blockdim2
358 : : integer :: spacedim
359 : : integer :: iblock, nblock
360 : : integer :: iline, nline
361 : : integer :: rows_tmp, cols_tmp, nband_eff, iband_min, iband_max
362 : : type(xgBlock_t) :: eigenBlock !
363 : : type(xgBlock_t) :: residuBlock,occBlock
364 : : double precision :: maxResidu, minResidu, dummy
365 : : double precision :: dlamch,tolerance
366 : : integer :: ierr = 0
367 : : integer :: nrestart
368 : : double precision :: tsec(2)
369 : : character(len=500) :: msg
370 : :
371 : : interface
372 : : subroutine getAX_BX(X,AX,BX)
373 : : use m_xg, only : xgBlock_t
374 : : type(xgBlock_t), intent(inout) :: X
375 : : type(xgBlock_t), intent(inout) :: AX
376 : : type(xgBlock_t), intent(inout) :: BX
377 : : end subroutine getAX_BX
378 : : end interface
379 : :
380 : : ! call timab(tim_run,1,tsec)
381 : :
382 : 24171 : lobpcg%prtvol = prtvol
383 : :
384 : 48342 : tolerance=2*dlamch('E')
385 : :
386 : 24171 : blockdim = lobpcg%blockdim
387 : 24171 : blockdim2 = 2*blockdim
388 : 24171 : blockdim3 = 3*blockdim
389 : 24171 : spacedim = lobpcg%spacedim
390 : :
391 : 24171 : nblock = lobpcg%nblock
392 : 24171 : nline = lobpcg%nline
393 : :
394 [ + + ]: 24171 : if (nbdbuf>0) then
395 : 7392 : nband_eff = lobpcg%neigenpairs - nbdbuf
396 : : else
397 : 16779 : nband_eff = lobpcg%neigenpairs
398 : : end if
399 : :
400 : 24171 : call xgBlock_getSize(eigen,rows_tmp, cols_tmp)
401 [ - + - - ]: 24171 : if ( rows_tmp /= lobpcg%neigenpairs .and. cols_tmp /= 1 ) then
402 : 0 : ABI_ERROR("Error eigen size")
403 : : endif
404 : 24171 : call xgBlock_getSize(X0,rows_tmp, cols_tmp)
405 [ - + ]: 24171 : if ( rows_tmp /= lobpcg%spacedim ) then
406 : 0 : ABI_ERROR("Error X0 spacedim")
407 : : endif
408 [ - + ]: 24171 : if ( cols_tmp /= lobpcg%neigenpairs ) then
409 : 0 : ABI_ERROR("Error X0 npairs")
410 : : endif
411 : :
412 [ + + + + : 24171 : if (isppol==1.and.ikpt==1.and.inonsc==1.and.istep==1) then
+ + + + ]
413 : 691 : write(msg,'(a,es16.6)') ' lobpcg%tolerance(tolwfr_diago)=',lobpcg%tolerance
414 : 691 : call wrtout(std_out,msg,'COLL')
415 : : end if
416 : :
417 : 24171 : call xg_init(eigenvalues3N,SPACE_R,blockdim3,1, gpu_option=lobpcg%gpu_option)
418 : 24171 : call xg_setBlock(eigenvalues3N,eigenvaluesN,blockdim,1)
419 : 24171 : call xg_setBlock(eigenvalues3N,eigenvalues2N,blockdim2,1)
420 : :
421 : 24171 : call xgBlock_reshape(eigen,blockdim,nblock)
422 : 24171 : call xgBlock_reshape(residu,blockdim,nblock)
423 : 24171 : call xgBlock_reshape(occ,blockdim,nblock)
424 : :
425 : 24171 : lobpcg%AllX0 = X0
426 : :
427 : 24171 : call xg_init(residu_eff,SPACE_R,blockdim,1,gpu_option=ABI_GPU_DISABLED)
428 : :
429 [ + + ]: 24171 : if ( lobpcg%paral_kgb == 1 ) then
430 : 21672 : call timab(tim_transpose,1,tsec)
431 : : call xgTransposer_constructor(lobpcg%xgTransposerX,lobpcg%X,lobpcg%XColsRows,nspinor,&
432 : : STATE_LINALG,TRANS_ALL2ALL,lobpcg%comm_rows,lobpcg%comm_cols,0,0,lobpcg%me_g0_fft,&
433 : 21672 : gpu_option=lobpcg%gpu_option,gpu_thread_limit=lobpcg%gpu_thread_limit)
434 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerAX,lobpcg%xgTransposerX,&
435 : 21672 : lobpcg%AX,lobpcg%AXColsRows,STATE_LINALG)
436 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerBX,lobpcg%xgTransposerX,&
437 : 21672 : lobpcg%BX,lobpcg%BXColsRows,STATE_LINALG)
438 : :
439 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerW,lobpcg%xgTransposerX,&
440 : 21672 : lobpcg%W,lobpcg%WColsRows,STATE_LINALG)
441 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerAW,lobpcg%xgTransposerX,&
442 : 21672 : lobpcg%AW,lobpcg%AWColsRows,STATE_LINALG)
443 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerBW,lobpcg%xgTransposerX,&
444 : 21672 : lobpcg%BW,lobpcg%BWColsRows,STATE_LINALG)
445 : 21672 : call timab(tim_transpose,2,tsec)
446 : : else
447 : 2499 : call xgBlock_setBlock(lobpcg%X, lobpcg%XColsRows, spacedim, blockdim)
448 : 2499 : call xgBlock_setBlock(lobpcg%AX, lobpcg%AXColsRows, spacedim, blockdim)
449 : 2499 : call xgBlock_setBlock(lobpcg%BX, lobpcg%BXColsRows, spacedim, blockdim)
450 : 2499 : call xgBlock_setBlock(lobpcg%W, lobpcg%WColsRows, spacedim, blockdim)
451 : 2499 : call xgBlock_setBlock(lobpcg%AW, lobpcg%AWColsRows, spacedim, blockdim)
452 : 2499 : call xgBlock_setBlock(lobpcg%BW, lobpcg%BWColsRows, spacedim, blockdim)
453 : : end if
454 : :
455 : : !! Start big loop over blocks
456 [ + + ]: 164009 : do iblock = 1, nblock
457 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_BLOCK)
458 : 139838 : nrestart = 0
459 : :
460 : 139838 : call lobpcg_getX0(lobpcg,iblock)
461 : 139838 : call xgBlock_setBlock(residu,residuBlock,blockdim,1,fcol=iblock)
462 : 139838 : call xgBlock_setBlock(occ, occBlock, blockdim,1,fcol=iblock)
463 : :
464 [ + + ]: 139838 : if ( iblock > 1 ) then
465 : 115667 : call lobpcg_setPreviousX0_BX0(lobpcg,iblock)
466 : :
467 : : ! Orthogonalize current iblock X block With Respect To previous Blocks in B-basis
468 : 115667 : call lobpcg_orthoXwrtBlocks(lobpcg,lobpcg%X,iblock)
469 : : end if
470 : :
471 [ + + ]: 139838 : if (lobpcg%paral_kgb == 1) then
472 : 132500 : call timab(tim_transpose,1,tsec)
473 : 132500 : call xgTransposer_transpose(lobpcg%xgTransposerX,STATE_COLSROWS)
474 : 132500 : lobpcg%xgTransposerAX%state=STATE_COLSROWS
475 : 132500 : lobpcg%xgTransposerBX%state=STATE_COLSROWS
476 : 132500 : call timab(tim_transpose,2,tsec)
477 : : end if
478 : : ! Initialize some quantitites (AX and BX)
479 : 139838 : call timab(tim_ax_bx,1,tsec)
480 : 139838 : call getAX_BX(lobpcg%XColsRows,lobpcg%AXColsRows,lobpcg%BXColsRows)
481 : 139838 : call xgBlock_zero_im_g0(lobpcg%AXColsRows)
482 : 139838 : call xgBlock_zero_im_g0(lobpcg%BXColsRows)
483 : 139838 : call timab(tim_ax_bx,2,tsec)
484 [ + + ]: 139838 : if (lobpcg%paral_kgb == 1) then
485 : 132500 : call timab(tim_transpose,1,tsec)
486 : 132500 : call xgTransposer_transpose(lobpcg%xgTransposerX,STATE_LINALG)
487 : 132500 : call xgTransposer_transpose(lobpcg%xgTransposerAX,STATE_LINALG)
488 : 132500 : call xgTransposer_transpose(lobpcg%xgTransposerBX,STATE_LINALG)
489 : 132500 : call timab(tim_transpose,2,tsec)
490 : : end if
491 : :
492 : : ! B-orthonormalize X, BX and AX
493 : 139838 : call xg_Borthonormalize(lobpcg%X,lobpcg%BX,ierr,tim_Bortho_X,lobpcg%gpu_option,AX=lobpcg%AX)
494 : :
495 : : ! Do first RR on X to get the first eigen values
496 : 139838 : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvaluesN,ierr,lobpcg%prtvol,tim_RR_X,lobpcg%gpu_option)
497 : :
498 : 139838 : compute_residu = .true.
499 : :
500 [ + + ]: 697252 : do iline = 1, nline
501 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_LINE)
502 : :
503 [ - + ]: 567950 : if ( ierr /= 0 ) then
504 : : !ABI_COMMENT("Consider using more bands and nbdbuf if necessary.")
505 : 0 : ierr = 0
506 : : end if
507 : :
508 : : !write(*,*) " -> Iteration ", iline
509 : :
510 : : ! Compute AX-Lambda*BX
511 : 567950 : call lobpcg_getResidu(lobpcg,eigenvaluesN)
512 : :
513 : : ! Compute residu norm here !
514 : 567950 : call timab(tim_maxres,1,tsec)
515 : 567950 : call xgBlock_colwiseNorm2(lobpcg%W,residuBlock)
516 : 567950 : call timab(tim_maxres,2,tsec)
517 : :
518 : : ! Apply preconditioner
519 : 567950 : call timab(tim_pcond,1,tsec)
520 : 567950 : call xgBlock_apply_diag(lobpcg%W,pcond,nspinor)
521 : 567950 : call timab(tim_pcond,2,tsec)
522 : :
523 : 567950 : call timab(tim_nbdbuf,1,tsec)
524 [ + + ]: 567950 : if (nbdbuf>=0) then
525 : : ! There is a transfer from GPU to CPU in this copy
526 : 566582 : call xgBlock_copy(residuBlock,residu_eff%self)
527 : 566582 : iband_min = 1 + blockdim*(iblock-1)
528 : 566582 : iband_max = blockdim*iblock
529 [ + + ]: 566582 : if (iband_max<=nband_eff) then ! all bands of this block are below nband_eff
530 : 537293 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu)
531 [ + + ]: 29289 : else if (iband_min<=nband_eff) then ! some bands of this block are below nband_eff
532 : 21961 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu,row_bound=(nband_eff-iband_min+1))
533 : : else ! all bands of this block are above nband_eff
534 : 7328 : minResidu = 0.0
535 : 7328 : maxResidu = 0.0
536 : : end if
537 [ + - ]: 1368 : else if (nbdbuf==-101) then
538 : 1368 : call xgBlock_minmax(residuBlock,minResidu,dummy) ! Get minimum of true residuals
539 : : ! Compute effective residuals : res_eff = res * occ
540 : 1368 : call xgBlock_apply_diag(residuBlock,occBlock,1,Y=residu_eff%self)
541 : 1368 : call xgBlock_minmax(residu_eff%self,dummy,maxResidu) ! Get maximum of effective residuals
542 : : else
543 : 0 : ABI_ERROR('Bad value of nbdbuf')
544 : : end if
545 : 567950 : call timab(tim_nbdbuf,2,tsec)
546 [ + + ]: 567950 : if ( maxResidu < lobpcg%tolerance ) then
547 : : compute_residu = .false.
548 : : ABI_NVTX_END_RANGE()
549 : : exit
550 : : end if
551 : :
552 : : ! Orthonormalize with respect to previous blocks
553 [ + + ]: 557414 : if ( iblock > 1 ) then
554 : 445151 : call lobpcg_orthoXwrtBlocks(lobpcg,lobpcg%W,iblock)
555 : : end if
556 : :
557 [ + + ]: 557414 : if (lobpcg%paral_kgb == 1) then
558 : 526716 : call timab(tim_transpose,1,tsec)
559 : 526716 : call xgTransposer_transpose(lobpcg%xgTransposerW,STATE_COLSROWS)
560 : 526716 : lobpcg%xgTransposerAW%state=STATE_COLSROWS
561 : 526716 : lobpcg%xgTransposerBW%state=STATE_COLSROWS
562 : 526716 : call timab(tim_transpose,2,tsec)
563 : : end if
564 : : ! Apply A and B on W
565 : 557414 : call timab(tim_ax_bx,1,tsec)
566 : 557414 : call getAX_BX(lobpcg%WColsRows,lobpcg%AWColsRows,lobpcg%BWColsRows)
567 : 557414 : call xgBlock_zero_im_g0(lobpcg%AWColsRows)
568 : 557414 : call xgBlock_zero_im_g0(lobpcg%BWColsRows)
569 : 557414 : call timab(tim_ax_bx,2,tsec)
570 [ + + ]: 557414 : if (lobpcg%paral_kgb == 1) then
571 : 526716 : call timab(tim_transpose,1,tsec)
572 : 526716 : call xgTransposer_transpose(lobpcg%xgTransposerW,STATE_LINALG)
573 : 526716 : call xgTransposer_transpose(lobpcg%xgTransposerAW,STATE_LINALG)
574 : 526716 : call xgTransposer_transpose(lobpcg%xgTransposerBW,STATE_LINALG)
575 : 526716 : call timab(tim_transpose,2,tsec)
576 : : end if
577 : :
578 : : ! DO RR in the correct subspace
579 : : ! if residu starts to be too small, there is an accumulation error in
580 : : ! P with values such as 1e-29 that make the eigenvectors diverge
581 [ + + + + ]: 686716 : if ( iline == 1 .or. minResidu < 1e-27) then
582 : : ! Do RR on XW to get the eigen vectors
583 : 134002 : call xg_Borthonormalize(lobpcg%XW,lobpcg%BXW,ierr,tim_Bortho_XW,lobpcg%gpu_option,AX=lobpcg%AXW) ! Do rotate AW
584 : 134002 : call xgBlock_zero(lobpcg%P)
585 : 134002 : call xgBlock_zero(lobpcg%AP)
586 : 134002 : call xgBlock_zero(lobpcg%BP)
587 [ - + ]: 134002 : if ( ierr /= 0 ) then
588 : 0 : ABI_COMMENT("B-orthonormalization (XW) did not work.")
589 : : end if
590 : : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvalues2N,ierr,lobpcg%prtvol,tim_RR_XW,lobpcg%gpu_option,&
591 : : & tolerance=tolerance,&
592 : : & XW=lobpcg%XW,AW=lobpcg%AW,BW=lobpcg%BW,P=lobpcg%P,AP=lobpcg%AP,BP=lobpcg%BP,WP=lobpcg%WP,&
593 : 134002 : & AWP=lobpcg%AWP,BWP=lobpcg%BWP)
594 [ - + ]: 134002 : if ( ierr /= 0 ) then
595 : 0 : ABI_WARNING("RayleighRitz (XW) did not work, but continue anyway.")
596 : : ABI_NVTX_END_RANGE()
597 : : exit
598 : : end if
599 : : else
600 : : ! B-orthonormalize P, BP
601 : 423412 : call xg_Borthonormalize(lobpcg%XWP%self,lobpcg%BXWP%self,ierr,tim_Bortho_XWP,lobpcg%gpu_option,AX=lobpcg%AXWP%self) ! Do rotate AW
602 : : ! Do RR on XWP to get the eigen vectors
603 [ + + ]: 423412 : if ( ierr == 0 ) then
604 : : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvalues3N%self,ierr,lobpcg%prtvol,tim_RR_XWP,lobpcg%gpu_option,&
605 : : & tolerance=tolerance,XW=lobpcg%XW,AW=lobpcg%AW,BW=lobpcg%BW,P=lobpcg%P,AP=lobpcg%AP,BP=lobpcg%BP,WP=lobpcg%WP,&
606 : 423404 : & AWP=lobpcg%AWP,BWP=lobpcg%BWP,XWP=lobpcg%XWP%self)
607 [ - + ]: 423404 : if ( ierr /= 0 ) then
608 : 0 : ABI_WARNING("RayleighRitz (XWP) did not work, but continue anyway.")
609 : : ABI_NVTX_END_RANGE()
610 : : exit
611 : : end if
612 : : else
613 : 8 : ABI_COMMENT("B-orthonormalization (XWP) did not work, try on XW.")
614 : 8 : call xg_Borthonormalize(lobpcg%XW,lobpcg%BXW,ierr,tim_Bortho_XW,lobpcg%gpu_option,AX=lobpcg%AXW) ! Do rotate AW
615 [ - + ]: 8 : if ( ierr /= 0 ) then
616 : 0 : ABI_COMMENT("B-orthonormalization (XW) did not work.")
617 : : end if
618 : 8 : call xgBlock_zero(lobpcg%P)
619 : 8 : call xgBlock_zero(lobpcg%AP)
620 : 8 : call xgBlock_zero(lobpcg%BP)
621 : 8 : nrestart = nrestart + 1
622 : : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvalues2N,ierr,lobpcg%prtvol,tim_RR_XW,lobpcg%gpu_option,&
623 : : & tolerance=tolerance,&
624 : : & XW=lobpcg%XW,AW=lobpcg%AW,BW=lobpcg%BW,P=lobpcg%P,AP=lobpcg%AP,BP=lobpcg%BP,WP=lobpcg%WP,&
625 : 8 : & AWP=lobpcg%AWP,BWP=lobpcg%BWP)
626 [ - + ]: 8 : if ( ierr /= 0 ) then
627 : 0 : ABI_WARNING("RayleighRitz (XWP) did not work, but continue anyway.")
628 : : ABI_NVTX_END_RANGE()
629 : : exit
630 : : end if
631 : : end if
632 : : end if
633 : :
634 : : ABI_NVTX_END_RANGE()
635 : : end do
636 : :
637 [ + + ]: 139838 : if ( compute_residu ) then
638 : : ! Recompute AX-Lambda*BX for the last time
639 : 129302 : call lobpcg_getResidu(lobpcg,eigenvaluesN)
640 : : ! Recompute residu norm here !
641 : 129302 : call timab(tim_maxres,1,tsec)
642 : 129302 : call xgBlock_colwiseNorm2(lobpcg%W,residuBlock)
643 : 129302 : call timab(tim_maxres,2,tsec)
644 : : ! Apply preconditioner
645 : 129302 : call timab(tim_pcond,1,tsec)
646 : 129302 : call xgBlock_apply_diag(lobpcg%W,pcond,nspinor)
647 : 129302 : call timab(tim_pcond,2,tsec)
648 : :
649 : 129302 : call timab(tim_nbdbuf,1,tsec)
650 [ - + ]: 129302 : if(lobpcg%gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(residuBlock)
651 [ + + ]: 129302 : if (nbdbuf>=0) then
652 : 129147 : call xgBlock_copy(residuBlock,residu_eff%self)
653 : 129147 : iband_min = 1 + blockdim*(iblock-1)
654 : 129147 : iband_max = blockdim*iblock
655 [ + + ]: 129147 : if (iband_max<=nband_eff) then ! all bands of this block are below nband_eff
656 : 126346 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu)
657 [ + - ]: 2801 : else if (iband_min<=nband_eff) then ! some bands of this block are below nband_eff
658 : 2801 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu,row_bound=(nband_eff-iband_min+1))
659 : : else ! all bands of this block are above nband_eff
660 : 0 : minResidu = 0.0
661 : 0 : maxResidu = 0.0
662 : : end if
663 [ + - ]: 155 : else if (nbdbuf==-101) then
664 : 155 : call xgBlock_minmax(residuBlock,minResidu,dummy) ! Get minimum of true residuals
665 : : ! Compute effective residuals : res_eff = res * occ
666 : 155 : call xgBlock_apply_diag(residuBlock,occBlock,1,Y=residu_eff%self)
667 : 155 : call xgBlock_minmax(residu_eff%self,dummy,maxResidu) ! Get maximum of effective residuals
668 : : else
669 : 0 : ABI_ERROR('Bad value of nbdbuf')
670 : : end if
671 : 129302 : call timab(tim_nbdbuf,2,tsec)
672 : : end if
673 : :
674 [ + + + - ]: 139838 : if (prtvol==5.and.xmpi_comm_rank(lobpcg%spacecom)==0) then
675 : 1002 : write(msg,'(6(a,i4),2(a,es16.6))') 'lobpcg | istep=',istep,'| isppol=',isppol,'| ikpt=',ikpt,&
676 : 2004 : & '| inonsc=',inonsc,'| iblock=',iblock,'| nline_done=',iline-1,'| minRes=',minResidu,'| maxRes=',maxResidu
677 : 1002 : call wrtout(std_out,msg,'PERS')
678 : : end if
679 : :
680 : : ! Save eigenvalues
681 : 139838 : call timab(tim_copy,1,tsec)
682 : 139838 : call xgBlock_setBlock(eigen,eigenBlock,blockdim,1,fcol=iblock)
683 : 139838 : call xgBlock_copy(eigenvaluesN,eigenBlock)
684 : 139838 : call timab(tim_copy,2,tsec)
685 : :
686 : : ! Save new X in X0
687 : 139838 : call lobpcg_setX0(lobpcg,iblock)
688 : :
689 : : ! Copy previous BX into BX0 for previous block
690 [ + + ]: 164009 : if ( nblock > 1 ) then
691 : 133478 : call lobpcg_transferAX_BX(lobpcg,iblock)
692 : : end if
693 : :
694 : : ABI_NVTX_END_RANGE()
695 : : end do !! End iblock loop
696 : :
697 : 24171 : call xgBlock_reshape(eigen,blockdim*nblock,1)
698 : 24171 : call xgBlock_reshape(residu,blockdim*nblock,1)
699 : 24171 : call xgBlock_reshape(occ,blockdim*nblock,1)
700 : :
701 : 24171 : call xg_free(eigenvalues3N)
702 : 24171 : call xg_free(residu_eff)
703 : :
704 : 24171 : skip = .false.
705 [ - + ]: 24171 : if ( ierr /= 0 ) then
706 : 0 : ABI_COMMENT("Some errors happened, so H|Psi> and S|Psi> are computed before leaving")
707 [ # # ]: 0 : if ( lobpcg%paral_kgb == 1 ) then
708 : : call xgTransposer_constructor(lobpcg%xgTransposerAllX0,lobpcg%AllX0,lobpcg%AllX0ColsRows,nspinor,&
709 : : STATE_LINALG,TRANS_ALL2ALL,lobpcg%comm_rows,lobpcg%comm_cols,0,0,lobpcg%me_g0_fft,&
710 : 0 : gpu_option=lobpcg%gpu_option,gpu_thread_limit=lobpcg%gpu_thread_limit)
711 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerAllAX0,lobpcg%xgTransposerAllX0,&
712 : 0 : lobpcg%AllAX0%self,lobpcg%AllAX0ColsRows,STATE_LINALG)
713 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerAllBX0,lobpcg%xgTransposerAllX0,&
714 : 0 : lobpcg%AllBX0%self,lobpcg%AllBX0ColsRows,STATE_LINALG)
715 : : else
716 : 0 : call xgBlock_setBlock(lobpcg%AllX0 , lobpcg%AllX0ColsRows , spacedim, lobpcg%neigenpairs)
717 : 0 : call xgBlock_setBlock(lobpcg%AllAX0%self, lobpcg%AllAX0ColsRows, spacedim, lobpcg%neigenpairs)
718 : 0 : call xgBlock_setBlock(lobpcg%AllBX0%self, lobpcg%AllBX0ColsRows, spacedim, lobpcg%neigenpairs)
719 : : end if
720 [ # # ]: 0 : if (lobpcg%paral_kgb == 1) then
721 : 0 : call timab(tim_transpose,1,tsec)
722 : 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllX0,STATE_COLSROWS)
723 : 0 : lobpcg%xgTransposerAllAX0%state=STATE_COLSROWS
724 : 0 : lobpcg%xgTransposerAllBX0%state=STATE_COLSROWS
725 : 0 : call timab(tim_transpose,2,tsec)
726 : : end if
727 : 0 : call timab(tim_ax_bx,1,tsec)
728 : 0 : call getAX_BX(lobpcg%AllX0ColsRows,lobpcg%AllAX0ColsRows,lobpcg%AllBX0ColsRows)
729 : 0 : call xgBlock_zero_im_g0(lobpcg%AllAX0ColsRows)
730 : 0 : call xgBlock_zero_im_g0(lobpcg%AllBX0ColsRows)
731 : 0 : call timab(tim_ax_bx,2,tsec)
732 [ # # ]: 0 : if (lobpcg%paral_kgb == 1) then
733 : 0 : call timab(tim_transpose,1,tsec)
734 : 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllX0,STATE_LINALG)
735 : 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllAX0,STATE_LINALG)
736 : 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllBX0,STATE_LINALG)
737 : 0 : call timab(tim_transpose,2,tsec)
738 : : end if
739 : 0 : call xgTransposer_free(lobpcg%xgTransposerAllX0)
740 : 0 : call xgTransposer_free(lobpcg%xgTransposerAllAX0)
741 : 0 : call xgTransposer_free(lobpcg%xgTransposerAllBX0)
742 : : skip = .true.
743 : : end if
744 : :
745 : : if (.not.skip) then
746 [ + + ]: 24171 : if ( nblock > 1 ) then
747 : : call xg_Borthonormalize(X0,lobpcg%AllBX0%self,ierr,tim_Bortho_Xall,&
748 : 17811 : & lobpcg%gpu_option,AX=lobpcg%AllAX0%self) ! Do rotate AX
749 : : call xg_RayleighRitz(X0,lobpcg%AllAX0%self,lobpcg%AllBX0%self,eigen,ierr,lobpcg%prtvol,tim_RR_Xall,&
750 : 17811 : & lobpcg%gpu_option,tolerance=tolerance)
751 : : end if
752 : : end if
753 : :
754 [ + + ]: 24171 : if ( lobpcg%paral_kgb == 1 ) then
755 : 21672 : call xgTransposer_free(lobpcg%xgTransposerX)
756 : 21672 : call xgTransposer_free(lobpcg%xgTransposerAX)
757 : 21672 : call xgTransposer_free(lobpcg%xgTransposerBX)
758 : 21672 : call xgTransposer_free(lobpcg%xgTransposerW)
759 : 21672 : call xgTransposer_free(lobpcg%xgTransposerAW)
760 : 21672 : call xgTransposer_free(lobpcg%xgTransposerBW)
761 : : end if
762 : :
763 : : ! call timab(tim_run,2,tsec)
764 : :
765 : 24171 : end subroutine lobpcg_run
766 : :
767 : :
768 : 139838 : subroutine lobpcg_getX0(lobpcg,iblock)
769 : :
770 : : type(lobpcg_t), intent(inout) :: lobpcg
771 : : integer , intent(in ) :: iblock
772 : : integer :: blockdim
773 : : integer :: spacedim
774 : : double precision :: tsec(2)
775 : :
776 : 139838 : call timab(tim_copy,1,tsec)
777 : :
778 : 139838 : blockdim = lobpcg%blockdim
779 : 139838 : spacedim = lobpcg%spacedim
780 : :
781 : : !lobpcg%XWP(:,X+1:X+blockdim) = lobpcg%X0(:,(iblock-1)*blockdim+1:iblock*blockdim)
782 : 139838 : call xgBlock_setBlock(lobpcg%AllX0,lobpcg%X0,spacedim,blockdim,fcol=(iblock-1)*blockdim+1)
783 : 139838 : call xgBlock_copy(lobpcg%X0,lobpcg%X)
784 : :
785 : 139838 : call timab(tim_copy,2,tsec)
786 : :
787 : 139838 : end subroutine lobpcg_getX0
788 : :
789 : :
790 : 115667 : subroutine lobpcg_setPreviousX0_BX0(lobpcg,iblock)
791 : :
792 : : type(lobpcg_t) , intent(inout) :: lobpcg
793 : : integer , intent(in ) :: iblock
794 : :
795 [ - + ]: 115667 : if (iblock<2) then
796 : 0 : ABI_ERROR("iblock<2")
797 : : end if
798 : 115667 : call xg_setBlock(lobpcg%AllBX0,lobpcg%BX0,lobpcg%spacedim,(iblock-1)*lobpcg%blockdim)
799 : 115667 : call xgBlock_setBlock(lobpcg%AllX0,lobpcg%X0,lobpcg%spacedim,(iblock-1)*lobpcg%blockdim)
800 : 115667 : end subroutine lobpcg_setPreviousX0_BX0
801 : :
802 : :
803 : 560818 : subroutine lobpcg_orthoXwrtBlocks(lobpcg,var,iblock)
804 : :
805 : : type(lobpcg_t) , intent(inout) :: lobpcg
806 : : type(xgBlock_t), intent(inout) :: var
807 : : integer , intent(in ) :: iblock
808 : : integer :: previousBlock
809 : : integer :: blockdim
810 : : integer :: spacedim
811 : : integer :: space_buf
812 : : type(xg_t) :: buffer
813 : : double precision :: tsec(2)
814 : :
815 : 560818 : call timab(tim_ortho,1,tsec)
816 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_ORTHO_X_WRT)
817 : :
818 : 560818 : blockdim = lobpcg%blockdim
819 : 560818 : spacedim = lobpcg%spacedim
820 : 560818 : previousBlock = (iblock-1)*lobpcg%blockdim
821 : :
822 : 560818 : space_buf = space(var)
823 [ + + ]: 560818 : if (space(var)==SPACE_CR) then
824 : 22654 : space_buf = SPACE_R
825 : : end if
826 : 560818 : call xg_init(buffer,space_buf,previousBlock,blockdim,comm=lobpcg%spacecom,gpu_option=lobpcg%gpu_option)
827 : :
828 : : ! buffer = BX0^T*X
829 : 560818 : call xgBlock_gemm('t','n',1.0d0,lobpcg%BX0,var,0.d0,buffer%self,comm=lobpcg%spacecom)
830 : :
831 : : ! sum all process contribution
832 : : ! X = - X0*(BX0^T*X) + X
833 : 560818 : call xgBlock_gemm('n','n',-1.0d0,lobpcg%X0,buffer%self,1.0d0,var)
834 : :
835 : 560818 : call xg_free(buffer)
836 : :
837 : : ABI_NVTX_END_RANGE()
838 : 560818 : call timab(tim_ortho,2,tsec)
839 : :
840 : 560818 : end subroutine lobpcg_orthoXwrtBlocks
841 : :
842 : 697252 : subroutine lobpcg_getResidu(lobpcg,eigenvalues)
843 : :
844 : : type(lobpcg_t) , intent(inout) :: lobpcg
845 : : type(xgBlock_t), intent(in ) :: eigenvalues
846 : : double precision :: tsec(2)
847 : :
848 : 697252 : call timab(tim_maxres,1,tsec)
849 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_RESIDUE)
850 : : !lobpcg%XWP(1:spacedim,shiftW+iblock) = lobpcg%AXWP(:,shiftX+iblock) - lobpcg%BXWP(:,shiftX+iblock)*eigenvalues(iblock)
851 : 697252 : call xgBlock_colwiseCymax(lobpcg%W,eigenvalues,lobpcg%BX,lobpcg%AX)
852 : : ABI_NVTX_END_RANGE()
853 : 697252 : call timab(tim_maxres,2,tsec)
854 : 697252 : end subroutine lobpcg_getResidu
855 : :
856 : 139838 : subroutine lobpcg_setX0(lobpcg,iblock)
857 : :
858 : : type(lobpcg_t) , intent(inout) :: lobpcg
859 : : integer , intent(in ) :: iblock
860 : : type(xgBlock_t) :: Xtmp
861 : : integer :: blockdim
862 : : integer :: spacedim
863 : : double precision :: tsec(2)
864 : :
865 : 139838 : call timab(tim_copy,1,tsec)
866 : 139838 : blockdim = lobpcg%blockdim
867 : 139838 : spacedim = lobpcg%spacedim
868 : :
869 : : !X0(:,(iblock-1)*blockdim+1:iblock*blockdim) = lobpcg%XWP(:,lobpcg%X+1:lobpcg%X+blockdim)
870 : 139838 : call xgBlock_setBlock(lobpcg%AllX0,Xtmp,spacedim,blockdim,fcol=(iblock-1)*blockdim+1)
871 : 139838 : call xgBlock_copy(lobpcg%X,Xtmp)
872 : 139838 : call timab(tim_copy,2,tsec)
873 : :
874 : 139838 : end subroutine lobpcg_setX0
875 : :
876 : :
877 : 133478 : subroutine lobpcg_transferAX_BX(lobpcg,iblock)
878 : :
879 : : type(lobpcg_t), intent(inout) :: lobpcg
880 : : integer , intent(in ) :: iblock
881 : : type(xgBlock_t) :: CXtmp
882 : : integer :: firstcol
883 : : double precision :: tsec(2)
884 : :
885 : 133478 : call timab(tim_copy,1,tsec)
886 : :
887 [ - + ]: 133478 : if (iblock<1) then
888 : 0 : ABI_ERROR("iblock<1")
889 : : end if
890 : :
891 : : ! iblock goes from 1 to nblock-1 included
892 : 133478 : firstcol = (iblock-1)*lobpcg%blockdim+1 ! Start of each block
893 : :
894 : : ! BX
895 : 133478 : call xg_setBlock(lobpcg%AllBX0,CXtmp,lobpcg%spacedim,lobpcg%blockdim,fcol=firstcol)
896 : 133478 : call xgBlock_copy(lobpcg%BX,CXtmp)
897 : :
898 : : ! AX
899 : 133478 : call xg_setBlock(lobpcg%AllAX0,CXtmp,lobpcg%spacedim,lobpcg%blockdim,fcol=firstcol)
900 : 133478 : call xgBlock_copy(lobpcg%AX,CXtmp)
901 : :
902 : 133478 : call timab(tim_copy,2,tsec)
903 : :
904 : 133478 : end subroutine lobpcg_transferAX_BX
905 : :
906 : 48342 : subroutine lobpcg_free(lobpcg)
907 : :
908 : : type(lobpcg_t), intent(inout) :: lobpcg
909 : :
910 : 48342 : call xg_free(lobpcg%XWP)
911 : 48342 : call xg_free(lobpcg%AXWP)
912 : 48342 : call xg_free(lobpcg%BXWP)
913 : 48342 : call xg_free(lobpcg%AllAX0)
914 : 48342 : call xg_free(lobpcg%AllBX0)
915 : 48342 : end subroutine lobpcg_free
916 : :
917 [ # # # # : 0 : end module m_lobpcg2
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
918 : : !!***
|
