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 : 24188 : 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 : 24188 : call timab(tim_init,1,tsec)
198 : 24188 : lobpcg%neigenpairs = neigenpairs
199 : 24188 : lobpcg%spacedim = spacedim
200 : 24188 : lobpcg%blockdim = blockdim
201 [ + + ]: 24188 : if (tolerance > 0.0) then
202 : 4568 : lobpcg%tolerance = tolerance
203 : : else
204 : 19620 : lobpcg%tolerance = 1.0e-20
205 : : end if
206 : 24188 : lobpcg%nline = nline
207 : 24188 : lobpcg%spacecom = spacecom
208 : 24188 : lobpcg%nblock = neigenpairs / blockdim
209 : 24188 : lobpcg%paral_kgb = paral_kgb
210 : 24188 : lobpcg%comm_rows = comm_rows
211 : 24188 : lobpcg%comm_cols = comm_cols
212 : 24188 : lobpcg%me_g0 = me_g0
213 : 24188 : lobpcg%me_g0_fft = me_g0_fft
214 : 24188 : lobpcg%gpu_option = gpu_option
215 : 24188 : lobpcg%gpu_thread_limit = 0
216 [ + - ]: 24188 : if(present(gpu_thread_limit)) lobpcg%gpu_thread_limit = gpu_thread_limit
217 : :
218 : 24188 : call lobpcg_allocateAll(lobpcg,space,me_g0)
219 : 24188 : call timab(tim_init,2,tsec)
220 : :
221 : 24188 : end subroutine lobpcg_init
222 : :
223 : :
224 : 24188 : 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 : 24188 : spacedim = lobpcg%spacedim
233 : 24188 : blockdim = lobpcg%blockdim
234 : :
235 : 24188 : call lobpcg_free(lobpcg) ! Make sure everything is not allocated and
236 : : ! pointer point to null()
237 : :
238 : 24188 : call xg_init(lobpcg%XWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
239 : 24188 : call xg_setBlock(lobpcg%XWP,lobpcg%X,spacedim,blockdim)
240 : 24188 : call xg_setBlock(lobpcg%XWP,lobpcg%W,spacedim,blockdim,fcol=blockdim+1)
241 : 24188 : call xg_setBlock(lobpcg%XWP,lobpcg%P,spacedim,blockdim,fcol=2*blockdim+1)
242 : 24188 : call xg_setBlock(lobpcg%XWP,lobpcg%XW,spacedim,2*blockdim)
243 : 24188 : call xg_setBlock(lobpcg%XWP,lobpcg%WP,spacedim,2*blockdim,fcol=blockdim+1)
244 : :
245 : 24188 : call xg_init(lobpcg%AXWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
246 : 24188 : call xg_setBlock(lobpcg%AXWP,lobpcg%AX,spacedim,blockdim)
247 : 24188 : call xg_setBlock(lobpcg%AXWP,lobpcg%AW,spacedim,blockdim,fcol=blockdim+1)
248 : 24188 : call xg_setBlock(lobpcg%AXWP,lobpcg%AP,spacedim,blockdim,fcol=2*blockdim+1)
249 : 24188 : call xg_setBlock(lobpcg%AXWP,lobpcg%AXW,spacedim,2*blockdim)
250 : 24188 : call xg_setBlock(lobpcg%AXWP,lobpcg%AWP,spacedim,2*blockdim,fcol=blockdim+1)
251 : :
252 : 24188 : call xg_init(lobpcg%BXWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
253 : 24188 : call xg_setBlock(lobpcg%BXWP,lobpcg%BX,spacedim,blockdim)
254 : 24188 : call xg_setBlock(lobpcg%BXWP,lobpcg%BW,spacedim,blockdim,fcol=blockdim+1)
255 : 24188 : call xg_setBlock(lobpcg%BXWP,lobpcg%BP,spacedim,blockdim,fcol=2*blockdim+1)
256 : 24188 : call xg_setBlock(lobpcg%BXWP,lobpcg%BXW,spacedim,2*blockdim)
257 : 24188 : call xg_setBlock(lobpcg%BXWP,lobpcg%BWP,spacedim,2*blockdim,fcol=blockdim+1)
258 : :
259 [ + + ]: 24188 : 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 : 6377 : lobpcg%AllBX0%self = lobpcg%BX
264 : 6377 : lobpcg%AllAX0%self = lobpcg%AX
265 : : end if
266 : :
267 : :
268 : 24188 : 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 : 72564 : subroutine lobpcg_run(lobpcg, X0, getAX_BX, pcond, eigen, occ, residu, prtvol, nspinor, isppol, ikpt, inonsc, istep, &
341 : : nbdbuf,pcond_nspinor)
342 : :
343 : : type(lobpcg_t) , intent(inout) :: lobpcg
344 : : type(xgBlock_t), intent(inout) :: X0 ! Full initial vectors
345 : : type(xgBlock_t), intent(inout) :: eigen ! Full initial eigen values
346 : : type(xgBlock_t), intent(inout) :: occ
347 : : type(xgBlock_t), intent(inout) :: residu
348 : : type(xgBlock_t), intent(in) :: pcond
349 : : integer , intent(in ) :: prtvol
350 : : integer , intent(in ) :: nspinor
351 : : integer , intent(in ) :: isppol,ikpt,inonsc,istep,nbdbuf
352 : : integer , intent(in ),optional :: pcond_nspinor
353 : :
354 : : type(xg_t) :: eigenvalues3N ! eigen values for Rayleight-Ritz
355 : : type(xg_t) :: residu_eff
356 : : type(xgBlock_t) :: eigenvaluesN ! eigen values for Rayleight-Ritz
357 : : type(xgBlock_t) :: eigenvalues2N ! eigen values for Rayleight-Ritz
358 : : logical :: skip,compute_residu,drop_null
359 : : integer :: blockdim, blockdim3, blockdim2
360 : : integer :: spacedim
361 : : integer :: iblock, nblock,nspinor_pcon
362 : : integer :: iline, nline
363 : : integer :: rows_tmp, cols_tmp, nband_eff, iband_min, iband_max
364 : : type(xgBlock_t) :: eigenBlock !
365 : : type(xgBlock_t) :: residuBlock,occBlock
366 : : double precision :: maxResidu, minResidu, dummy
367 : : double precision :: dlamch,tolerance
368 : : integer :: ierr = 0
369 : : integer :: nrestart
370 : : double precision :: tsec(2)
371 : : character(len=500) :: msg
372 : :
373 : : interface
374 : : subroutine getAX_BX(X,AX,BX)
375 : : use m_xg, only : xgBlock_t
376 : : type(xgBlock_t), intent(inout) :: X
377 : : type(xgBlock_t), intent(inout) :: AX
378 : : type(xgBlock_t), intent(inout) :: BX
379 : : end subroutine getAX_BX
380 : : end interface
381 : :
382 : : ! call timab(tim_run,1,tsec)
383 : :
384 : 24188 : lobpcg%prtvol = prtvol
385 : :
386 : 24188 : nspinor_pcon=nspinor
387 : 24188 : if (present(pcond_nspinor)) nspinor_pcon=pcond_nspinor
388 : :
389 : 24188 : tolerance=2*dlamch('E')
390 : :
391 : 24188 : blockdim = lobpcg%blockdim
392 : 24188 : blockdim2 = 2*blockdim
393 : 24188 : blockdim3 = 3*blockdim
394 : 24188 : spacedim = lobpcg%spacedim
395 : :
396 : 24188 : nblock = lobpcg%nblock
397 : 24188 : nline = lobpcg%nline
398 : :
399 [ + + ]: 24188 : if (nbdbuf>0) then
400 : 7392 : nband_eff = lobpcg%neigenpairs - nbdbuf
401 : : else
402 : 16796 : nband_eff = lobpcg%neigenpairs
403 : : end if
404 : :
405 : 24188 : call xgBlock_getSize(eigen,rows_tmp, cols_tmp)
406 [ - + - - ]: 24188 : if ( rows_tmp /= lobpcg%neigenpairs .and. cols_tmp /= 1 ) then
407 : 0 : ABI_ERROR("Error eigen size")
408 : : endif
409 : 24188 : call xgBlock_getSize(X0,rows_tmp, cols_tmp)
410 [ - + ]: 24188 : if ( rows_tmp /= lobpcg%spacedim ) then
411 : 0 : ABI_ERROR("Error X0 spacedim")
412 : : endif
413 [ - + ]: 24188 : if ( cols_tmp /= lobpcg%neigenpairs ) then
414 : 0 : ABI_ERROR("Error X0 npairs")
415 : : endif
416 : :
417 [ + + + + : 24188 : if (isppol==1.and.ikpt==1.and.inonsc==1.and.istep==1) then
+ + + + ]
418 : 692 : write(msg,'(a,es16.6)') ' lobpcg%tolerance(tolwfr_diago)=',lobpcg%tolerance
419 : 692 : call wrtout(std_out,msg,'COLL')
420 : : end if
421 : :
422 : 24188 : call xg_init(eigenvalues3N,SPACE_R,blockdim3,1, gpu_option=lobpcg%gpu_option)
423 : 24188 : call xg_setBlock(eigenvalues3N,eigenvaluesN,blockdim,1)
424 : 24188 : call xg_setBlock(eigenvalues3N,eigenvalues2N,blockdim2,1)
425 : :
426 : 24188 : call xgBlock_reshape(eigen,blockdim,nblock)
427 : 24188 : call xgBlock_reshape(residu,blockdim,nblock)
428 : 24188 : call xgBlock_reshape(occ,blockdim,nblock)
429 : :
430 : 24188 : lobpcg%AllX0 = X0
431 : :
432 : 24188 : call xg_init(residu_eff,SPACE_R,blockdim,1,gpu_option=ABI_GPU_DISABLED)
433 : :
434 [ + + ]: 24188 : if ( lobpcg%paral_kgb == 1 ) then
435 : 21672 : call timab(tim_transpose,1,tsec)
436 : : call xgTransposer_constructor(lobpcg%xgTransposerX,lobpcg%X,lobpcg%XColsRows,nspinor,&
437 : : STATE_LINALG,TRANS_ALL2ALL,lobpcg%comm_rows,lobpcg%comm_cols,0,0,lobpcg%me_g0_fft,&
438 : 21672 : gpu_option=lobpcg%gpu_option,gpu_thread_limit=lobpcg%gpu_thread_limit)
439 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerAX,lobpcg%xgTransposerX,&
440 : 21672 : lobpcg%AX,lobpcg%AXColsRows,STATE_LINALG)
441 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerBX,lobpcg%xgTransposerX,&
442 : 21672 : lobpcg%BX,lobpcg%BXColsRows,STATE_LINALG)
443 : :
444 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerW,lobpcg%xgTransposerX,&
445 : 21672 : lobpcg%W,lobpcg%WColsRows,STATE_LINALG)
446 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerAW,lobpcg%xgTransposerX,&
447 : 21672 : lobpcg%AW,lobpcg%AWColsRows,STATE_LINALG)
448 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerBW,lobpcg%xgTransposerX,&
449 : 21672 : lobpcg%BW,lobpcg%BWColsRows,STATE_LINALG)
450 : 21672 : call timab(tim_transpose,2,tsec)
451 : : else
452 : 2516 : call xgBlock_setBlock(lobpcg%X, lobpcg%XColsRows, spacedim, blockdim)
453 : 2516 : call xgBlock_setBlock(lobpcg%AX, lobpcg%AXColsRows, spacedim, blockdim)
454 : 2516 : call xgBlock_setBlock(lobpcg%BX, lobpcg%BXColsRows, spacedim, blockdim)
455 : 2516 : call xgBlock_setBlock(lobpcg%W, lobpcg%WColsRows, spacedim, blockdim)
456 : 2516 : call xgBlock_setBlock(lobpcg%AW, lobpcg%AWColsRows, spacedim, blockdim)
457 : 2516 : call xgBlock_setBlock(lobpcg%BW, lobpcg%BWColsRows, spacedim, blockdim)
458 : : end if
459 : :
460 : : !! Start big loop over blocks
461 [ + + ]: 164043 : do iblock = 1, nblock
462 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_BLOCK)
463 : 139855 : nrestart = 0
464 : :
465 : 139855 : call lobpcg_getX0(lobpcg,iblock)
466 : 139855 : call xgBlock_setBlock(residu,residuBlock,blockdim,1,fcol=iblock)
467 : 139855 : call xgBlock_setBlock(occ, occBlock, blockdim,1,fcol=iblock)
468 : :
469 [ + + ]: 139855 : if ( iblock > 1 ) then
470 : 115667 : call lobpcg_setPreviousX0_BX0(lobpcg,iblock)
471 : :
472 : : ! Orthogonalize current iblock X block With Respect To previous Blocks in B-basis
473 : 115667 : call lobpcg_orthoXwrtBlocks(lobpcg,lobpcg%X,iblock)
474 : : end if
475 : :
476 [ + + ]: 139855 : if (lobpcg%paral_kgb == 1) then
477 : 132500 : call timab(tim_transpose,1,tsec)
478 : 132500 : call xgTransposer_transpose(lobpcg%xgTransposerX,STATE_COLSROWS)
479 : 132500 : lobpcg%xgTransposerAX%state=STATE_COLSROWS
480 : 132500 : lobpcg%xgTransposerBX%state=STATE_COLSROWS
481 : 132500 : call timab(tim_transpose,2,tsec)
482 : : end if
483 : : ! Initialize some quantitites (AX and BX)
484 : 139855 : call timab(tim_ax_bx,1,tsec)
485 : 139855 : call getAX_BX(lobpcg%XColsRows,lobpcg%AXColsRows,lobpcg%BXColsRows)
486 : 139855 : call xgBlock_zero_im_g0(lobpcg%AXColsRows)
487 : 139855 : call xgBlock_zero_im_g0(lobpcg%BXColsRows)
488 : 139855 : call timab(tim_ax_bx,2,tsec)
489 [ + + ]: 139855 : if (lobpcg%paral_kgb == 1) then
490 : 132500 : call timab(tim_transpose,1,tsec)
491 : 132500 : call xgTransposer_transpose(lobpcg%xgTransposerX,STATE_LINALG)
492 : 132500 : call xgTransposer_transpose(lobpcg%xgTransposerAX,STATE_LINALG)
493 : 132500 : call xgTransposer_transpose(lobpcg%xgTransposerBX,STATE_LINALG)
494 : 132500 : call timab(tim_transpose,2,tsec)
495 : : end if
496 : :
497 : : ! B-orthonormalize X, BX and AX
498 : 139855 : call xg_Borthonormalize(lobpcg%X,lobpcg%BX,ierr,tim_Bortho_X,lobpcg%gpu_option,AX=lobpcg%AX)
499 : :
500 : : ! Do first RR on X to get the first eigen values
501 : 139855 : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvaluesN,ierr,lobpcg%prtvol,tim_RR_X,lobpcg%gpu_option)
502 : :
503 : 139855 : compute_residu = .true.
504 : :
505 [ + + ]: 697325 : do iline = 1, nline
506 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_LINE)
507 : :
508 [ - + ]: 568013 : if ( ierr /= 0 ) then
509 : : !ABI_COMMENT("Consider using more bands and nbdbuf if necessary.")
510 : 0 : ierr = 0
511 : : end if
512 : :
513 : : !write(*,*) " -> Iteration ", iline
514 : :
515 : : ! Compute AX-Lambda*BX
516 : 568013 : call lobpcg_getResidu(lobpcg,eigenvaluesN)
517 : :
518 : : ! Compute residu norm here !
519 : 568013 : call timab(tim_maxres,1,tsec)
520 : 568013 : call xgBlock_colwiseNorm2(lobpcg%W,residuBlock)
521 : 568013 : call timab(tim_maxres,2,tsec)
522 : :
523 : : ! Apply preconditioner
524 : 568013 : call timab(tim_pcond,1,tsec)
525 : 568013 : call xgBlock_apply_diag(lobpcg%W,pcond,nspinor_pcon)
526 : 568013 : call timab(tim_pcond,2,tsec)
527 : :
528 : 568013 : call timab(tim_nbdbuf,1,tsec)
529 [ + + ]: 568013 : if (nbdbuf>=0) then
530 : : ! There is a transfer from GPU to CPU in this copy
531 : 566645 : call xgBlock_copy(residuBlock,residu_eff%self)
532 : 566645 : iband_min = 1 + blockdim*(iblock-1)
533 : 566645 : iband_max = blockdim*iblock
534 [ + + ]: 566645 : if (iband_max<=nband_eff) then ! all bands of this block are below nband_eff
535 : 537356 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu)
536 [ + + ]: 29289 : else if (iband_min<=nband_eff) then ! some bands of this block are below nband_eff
537 : 21961 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu,row_bound=(nband_eff-iband_min+1))
538 : : else ! all bands of this block are above nband_eff
539 : 7328 : minResidu = 0.0
540 : 7328 : maxResidu = 0.0
541 : : end if
542 [ + - ]: 1368 : else if (nbdbuf==-101) then
543 : 1368 : call xgBlock_minmax(residuBlock,minResidu,dummy) ! Get minimum of true residuals
544 : : ! Compute effective residuals : res_eff = res * occ
545 : 1368 : call xgBlock_apply_diag(residuBlock,occBlock,1,Y=residu_eff%self)
546 : 1368 : call xgBlock_minmax(residu_eff%self,dummy,maxResidu) ! Get maximum of effective residuals
547 : : else
548 : 0 : ABI_ERROR('Bad value of nbdbuf')
549 : : end if
550 : 568013 : call timab(tim_nbdbuf,2,tsec)
551 [ + + ]: 568013 : if ( maxResidu < lobpcg%tolerance ) then
552 : : compute_residu = .false.
553 : : ABI_NVTX_END_RANGE()
554 : : exit
555 : : end if
556 : :
557 : : ! Orthonormalize with respect to previous blocks
558 [ + + ]: 557470 : if ( iblock > 1 ) then
559 : 445151 : call lobpcg_orthoXwrtBlocks(lobpcg,lobpcg%W,iblock)
560 : : end if
561 : :
562 [ + + ]: 557470 : if (lobpcg%paral_kgb == 1) then
563 : 526716 : call timab(tim_transpose,1,tsec)
564 : 526716 : call xgTransposer_transpose(lobpcg%xgTransposerW,STATE_COLSROWS)
565 : 526716 : lobpcg%xgTransposerAW%state=STATE_COLSROWS
566 : 526716 : lobpcg%xgTransposerBW%state=STATE_COLSROWS
567 : 526716 : call timab(tim_transpose,2,tsec)
568 : : end if
569 : : ! Apply A and B on W
570 : 557470 : call timab(tim_ax_bx,1,tsec)
571 : 557470 : call getAX_BX(lobpcg%WColsRows,lobpcg%AWColsRows,lobpcg%BWColsRows)
572 : 557470 : call xgBlock_zero_im_g0(lobpcg%AWColsRows)
573 : 557470 : call xgBlock_zero_im_g0(lobpcg%BWColsRows)
574 : 557470 : call timab(tim_ax_bx,2,tsec)
575 [ + + ]: 557470 : if (lobpcg%paral_kgb == 1) then
576 : 526716 : call timab(tim_transpose,1,tsec)
577 : 526716 : call xgTransposer_transpose(lobpcg%xgTransposerW,STATE_LINALG)
578 : 526716 : call xgTransposer_transpose(lobpcg%xgTransposerAW,STATE_LINALG)
579 : 526716 : call xgTransposer_transpose(lobpcg%xgTransposerBW,STATE_LINALG)
580 : 526716 : call timab(tim_transpose,2,tsec)
581 : : end if
582 : :
583 : : ! DO RR in the correct subspace
584 : : ! if residu starts to be too small, there is an accumulation error in
585 : : ! P with values such as 1e-29 that make the eigenvectors diverge
586 [ + + + + ]: 686782 : if ( iline == 1 .or. minResidu < 1e-27) then
587 : : ! Do RR on XW to get the eigen vectors
588 : : call xg_Borthonormalize(lobpcg%XW,lobpcg%BXW,ierr,tim_Bortho_XW,lobpcg%gpu_option,AX=lobpcg%AXW,&
589 : 134019 : & allow_regularize=.true.,regularized=drop_null) ! Do rotate AW
590 : 134019 : call xgBlock_zero(lobpcg%P)
591 : 134019 : call xgBlock_zero(lobpcg%AP)
592 : 134019 : call xgBlock_zero(lobpcg%BP)
593 [ - + ]: 134019 : if ( ierr /= 0 ) then
594 : 0 : ABI_COMMENT("B-orthonormalization (XW) did not work.")
595 : : end if
596 : : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvalues2N,ierr,lobpcg%prtvol,tim_RR_XW,lobpcg%gpu_option,&
597 : : & tolerance=tolerance,&
598 : : & XW=lobpcg%XW,AW=lobpcg%AW,BW=lobpcg%BW,P=lobpcg%P,AP=lobpcg%AP,BP=lobpcg%BP,WP=lobpcg%WP,&
599 : 134019 : & AWP=lobpcg%AWP,BWP=lobpcg%BWP,drop_null=drop_null)
600 [ - + ]: 134019 : if ( ierr /= 0 ) then
601 : 0 : ABI_WARNING("RayleighRitz (XW) did not work, but continue anyway.")
602 : : ABI_NVTX_END_RANGE()
603 : : exit
604 : : end if
605 : : else
606 : : ! B-orthonormalize P, BP
607 : 423451 : call xg_Borthonormalize(lobpcg%XWP%self,lobpcg%BXWP%self,ierr,tim_Bortho_XWP,lobpcg%gpu_option,AX=lobpcg%AXWP%self) ! Do rotate AW
608 : : ! Do RR on XWP to get the eigen vectors
609 [ + + ]: 423451 : if ( ierr == 0 ) then
610 : : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvalues3N%self,ierr,lobpcg%prtvol,tim_RR_XWP,lobpcg%gpu_option,&
611 : : & tolerance=tolerance,XW=lobpcg%XW,AW=lobpcg%AW,BW=lobpcg%BW,P=lobpcg%P,AP=lobpcg%AP,BP=lobpcg%BP,WP=lobpcg%WP,&
612 : 423432 : & AWP=lobpcg%AWP,BWP=lobpcg%BWP,XWP=lobpcg%XWP%self)
613 [ - + ]: 423432 : if ( ierr /= 0 ) then
614 : 0 : ABI_WARNING("RayleighRitz (XWP) did not work, but continue anyway.")
615 : : ABI_NVTX_END_RANGE()
616 : : exit
617 : : end if
618 : : else
619 : 19 : ABI_COMMENT("B-orthonormalization (XWP) did not work, try on XW.")
620 : : call xg_Borthonormalize(lobpcg%XW,lobpcg%BXW,ierr,tim_Bortho_XW,lobpcg%gpu_option,AX=lobpcg%AXW,&
621 : 19 : & allow_regularize=.true.,regularized=drop_null) ! Do rotate AW
622 [ - + ]: 19 : if ( ierr /= 0 ) then
623 : 0 : ABI_COMMENT("B-orthonormalization (XW) did not work.")
624 : : end if
625 : 19 : call xgBlock_zero(lobpcg%P)
626 : 19 : call xgBlock_zero(lobpcg%AP)
627 : 19 : call xgBlock_zero(lobpcg%BP)
628 : 19 : nrestart = nrestart + 1
629 : : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvalues2N,ierr,lobpcg%prtvol,tim_RR_XW,lobpcg%gpu_option,&
630 : : & tolerance=tolerance,&
631 : : & XW=lobpcg%XW,AW=lobpcg%AW,BW=lobpcg%BW,P=lobpcg%P,AP=lobpcg%AP,BP=lobpcg%BP,WP=lobpcg%WP,&
632 : 19 : & AWP=lobpcg%AWP,BWP=lobpcg%BWP,drop_null=drop_null)
633 [ - + ]: 19 : if ( ierr /= 0 ) then
634 : 0 : ABI_WARNING("RayleighRitz (XWP) did not work, but continue anyway.")
635 : : ABI_NVTX_END_RANGE()
636 : : exit
637 : : end if
638 : : end if
639 : : end if
640 : :
641 : : ABI_NVTX_END_RANGE()
642 : : end do
643 : :
644 [ + + ]: 139855 : if ( compute_residu ) then
645 : : ! Recompute AX-Lambda*BX for the last time
646 : 129312 : call lobpcg_getResidu(lobpcg,eigenvaluesN)
647 : : ! Recompute residu norm here !
648 : 129312 : call timab(tim_maxres,1,tsec)
649 : 129312 : call xgBlock_colwiseNorm2(lobpcg%W,residuBlock)
650 : 129312 : call timab(tim_maxres,2,tsec)
651 : : ! Apply preconditioner
652 : 129312 : call timab(tim_pcond,1,tsec)
653 : 129312 : call xgBlock_apply_diag(lobpcg%W,pcond,nspinor_pcon)
654 : 129312 : call timab(tim_pcond,2,tsec)
655 : :
656 : 129312 : call timab(tim_nbdbuf,1,tsec)
657 [ - + ]: 129312 : if(lobpcg%gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(residuBlock)
658 [ + + ]: 129312 : if (nbdbuf>=0) then
659 : 129157 : call xgBlock_copy(residuBlock,residu_eff%self)
660 : 129157 : iband_min = 1 + blockdim*(iblock-1)
661 : 129157 : iband_max = blockdim*iblock
662 [ + + ]: 129157 : if (iband_max<=nband_eff) then ! all bands of this block are below nband_eff
663 : 126356 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu)
664 [ + - ]: 2801 : else if (iband_min<=nband_eff) then ! some bands of this block are below nband_eff
665 : 2801 : call xgBlock_minmax(residu_eff%self,minResidu,maxResidu,row_bound=(nband_eff-iband_min+1))
666 : : else ! all bands of this block are above nband_eff
667 : 0 : minResidu = 0.0
668 : 0 : maxResidu = 0.0
669 : : end if
670 [ + - ]: 155 : else if (nbdbuf==-101) then
671 : 155 : call xgBlock_minmax(residuBlock,minResidu,dummy) ! Get minimum of true residuals
672 : : ! Compute effective residuals : res_eff = res * occ
673 : 155 : call xgBlock_apply_diag(residuBlock,occBlock,1,Y=residu_eff%self)
674 : 155 : call xgBlock_minmax(residu_eff%self,dummy,maxResidu) ! Get maximum of effective residuals
675 : : else
676 : 0 : ABI_ERROR('Bad value of nbdbuf')
677 : : end if
678 : 129312 : call timab(tim_nbdbuf,2,tsec)
679 : : end if
680 : :
681 [ + + + - ]: 139855 : if (prtvol==5.and.xmpi_comm_rank(lobpcg%spacecom)==0) then
682 : 1002 : write(msg,'(6(a,i4),2(a,es16.6))') 'lobpcg | istep=',istep,'| isppol=',isppol,'| ikpt=',ikpt,&
683 : 2004 : & '| inonsc=',inonsc,'| iblock=',iblock,'| nline_done=',iline-1,'| minRes=',minResidu,'| maxRes=',maxResidu
684 : 1002 : call wrtout(std_out,msg,'PERS')
685 : : end if
686 : :
687 : : ! Save eigenvalues
688 : 139855 : call timab(tim_copy,1,tsec)
689 : 139855 : call xgBlock_setBlock(eigen,eigenBlock,blockdim,1,fcol=iblock)
690 : 139855 : call xgBlock_copy(eigenvaluesN,eigenBlock)
691 : 139855 : call timab(tim_copy,2,tsec)
692 : :
693 : : ! Save new X in X0
694 : 139855 : call lobpcg_setX0(lobpcg,iblock)
695 : :
696 : : ! Copy previous BX into BX0 for previous block
697 [ + + ]: 164043 : if ( nblock > 1 ) then
698 : 133478 : call lobpcg_transferAX_BX(lobpcg,iblock)
699 : : end if
700 : :
701 : : ABI_NVTX_END_RANGE()
702 : : end do !! End iblock loop
703 : :
704 : 24188 : call xgBlock_reshape(eigen,blockdim*nblock,1)
705 : 24188 : call xgBlock_reshape(residu,blockdim*nblock,1)
706 : 24188 : call xgBlock_reshape(occ,blockdim*nblock,1)
707 : :
708 : 24188 : call xg_free(eigenvalues3N)
709 : 24188 : call xg_free(residu_eff)
710 : :
711 : 24188 : skip = .false.
712 [ - + ]: 24188 : if ( ierr /= 0 ) then
713 : 0 : ABI_COMMENT("Some errors happened, so H|Psi> and S|Psi> are computed before leaving")
714 [ # # ]: 0 : if ( lobpcg%paral_kgb == 1 ) then
715 : : call xgTransposer_constructor(lobpcg%xgTransposerAllX0,lobpcg%AllX0,lobpcg%AllX0ColsRows,nspinor,&
716 : : STATE_LINALG,TRANS_ALL2ALL,lobpcg%comm_rows,lobpcg%comm_cols,0,0,lobpcg%me_g0_fft,&
717 : 0 : gpu_option=lobpcg%gpu_option,gpu_thread_limit=lobpcg%gpu_thread_limit)
718 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerAllAX0,lobpcg%xgTransposerAllX0,&
719 : 0 : lobpcg%AllAX0%self,lobpcg%AllAX0ColsRows,STATE_LINALG)
720 : : call xgTransposer_copyConstructor(lobpcg%xgTransposerAllBX0,lobpcg%xgTransposerAllX0,&
721 : 0 : lobpcg%AllBX0%self,lobpcg%AllBX0ColsRows,STATE_LINALG)
722 : : else
723 : 0 : call xgBlock_setBlock(lobpcg%AllX0 , lobpcg%AllX0ColsRows , spacedim, lobpcg%neigenpairs)
724 : 0 : call xgBlock_setBlock(lobpcg%AllAX0%self, lobpcg%AllAX0ColsRows, spacedim, lobpcg%neigenpairs)
725 : 0 : call xgBlock_setBlock(lobpcg%AllBX0%self, lobpcg%AllBX0ColsRows, spacedim, lobpcg%neigenpairs)
726 : : end if
727 [ # # ]: 0 : if (lobpcg%paral_kgb == 1) then
728 : 0 : call timab(tim_transpose,1,tsec)
729 : 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllX0,STATE_COLSROWS)
730 : 0 : lobpcg%xgTransposerAllAX0%state=STATE_COLSROWS
731 : 0 : lobpcg%xgTransposerAllBX0%state=STATE_COLSROWS
732 : 0 : call timab(tim_transpose,2,tsec)
733 : : end if
734 : 0 : call timab(tim_ax_bx,1,tsec)
735 : 0 : call getAX_BX(lobpcg%AllX0ColsRows,lobpcg%AllAX0ColsRows,lobpcg%AllBX0ColsRows)
736 : 0 : call xgBlock_zero_im_g0(lobpcg%AllAX0ColsRows)
737 : 0 : call xgBlock_zero_im_g0(lobpcg%AllBX0ColsRows)
738 : 0 : call timab(tim_ax_bx,2,tsec)
739 [ # # ]: 0 : if (lobpcg%paral_kgb == 1) then
740 : 0 : call timab(tim_transpose,1,tsec)
741 : 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllX0,STATE_LINALG)
742 : 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllAX0,STATE_LINALG)
743 : 0 : call xgTransposer_transpose(lobpcg%xgTransposerAllBX0,STATE_LINALG)
744 : 0 : call timab(tim_transpose,2,tsec)
745 : : ! The transposers are constructed only with paral_kgb=1: free them only in that case
746 : 0 : call xgTransposer_free(lobpcg%xgTransposerAllX0)
747 : 0 : call xgTransposer_free(lobpcg%xgTransposerAllAX0)
748 : 0 : call xgTransposer_free(lobpcg%xgTransposerAllBX0)
749 : : end if
750 : : skip = .true.
751 : : end if
752 : :
753 : : if (.not.skip) then
754 [ + + ]: 24188 : if ( nblock > 1 ) then
755 : : call xg_Borthonormalize(X0,lobpcg%AllBX0%self,ierr,tim_Bortho_Xall,&
756 : 17811 : & lobpcg%gpu_option,AX=lobpcg%AllAX0%self) ! Do rotate AX
757 : : call xg_RayleighRitz(X0,lobpcg%AllAX0%self,lobpcg%AllBX0%self,eigen,ierr,lobpcg%prtvol,tim_RR_Xall,&
758 : 17811 : & lobpcg%gpu_option,tolerance=tolerance)
759 : : end if
760 : : end if
761 : :
762 [ + + ]: 24188 : if ( lobpcg%paral_kgb == 1 ) then
763 : 21672 : call xgTransposer_free(lobpcg%xgTransposerX)
764 : 21672 : call xgTransposer_free(lobpcg%xgTransposerAX)
765 : 21672 : call xgTransposer_free(lobpcg%xgTransposerBX)
766 : 21672 : call xgTransposer_free(lobpcg%xgTransposerW)
767 : 21672 : call xgTransposer_free(lobpcg%xgTransposerAW)
768 : 21672 : call xgTransposer_free(lobpcg%xgTransposerBW)
769 : : end if
770 : :
771 : : ! call timab(tim_run,2,tsec)
772 : :
773 [ + - ]: 24188 : end subroutine lobpcg_run
774 : :
775 : :
776 : 139855 : subroutine lobpcg_getX0(lobpcg,iblock)
777 : :
778 : : type(lobpcg_t), intent(inout) :: lobpcg
779 : : integer , intent(in ) :: iblock
780 : : integer :: blockdim
781 : : integer :: spacedim
782 : : double precision :: tsec(2)
783 : :
784 : 139855 : call timab(tim_copy,1,tsec)
785 : :
786 : 139855 : blockdim = lobpcg%blockdim
787 : 139855 : spacedim = lobpcg%spacedim
788 : :
789 : : !lobpcg%XWP(:,X+1:X+blockdim) = lobpcg%X0(:,(iblock-1)*blockdim+1:iblock*blockdim)
790 : 139855 : call xgBlock_setBlock(lobpcg%AllX0,lobpcg%X0,spacedim,blockdim,fcol=(iblock-1)*blockdim+1)
791 : 139855 : call xgBlock_copy(lobpcg%X0,lobpcg%X)
792 : :
793 : 139855 : call timab(tim_copy,2,tsec)
794 : :
795 : 139855 : end subroutine lobpcg_getX0
796 : :
797 : :
798 : 115667 : subroutine lobpcg_setPreviousX0_BX0(lobpcg,iblock)
799 : :
800 : : type(lobpcg_t) , intent(inout) :: lobpcg
801 : : integer , intent(in ) :: iblock
802 : :
803 [ - + ]: 115667 : if (iblock<2) then
804 : 0 : ABI_ERROR("iblock<2")
805 : : end if
806 : 115667 : call xg_setBlock(lobpcg%AllBX0,lobpcg%BX0,lobpcg%spacedim,(iblock-1)*lobpcg%blockdim)
807 : 115667 : call xgBlock_setBlock(lobpcg%AllX0,lobpcg%X0,lobpcg%spacedim,(iblock-1)*lobpcg%blockdim)
808 : 115667 : end subroutine lobpcg_setPreviousX0_BX0
809 : :
810 : :
811 : 560818 : subroutine lobpcg_orthoXwrtBlocks(lobpcg,var,iblock)
812 : :
813 : : type(lobpcg_t) , intent(inout) :: lobpcg
814 : : type(xgBlock_t), intent(inout) :: var
815 : : integer , intent(in ) :: iblock
816 : : integer :: previousBlock
817 : : integer :: blockdim
818 : : integer :: spacedim
819 : : integer :: space_buf
820 : : type(xg_t) :: buffer
821 : : double precision :: tsec(2)
822 : :
823 : 560818 : call timab(tim_ortho,1,tsec)
824 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_ORTHO_X_WRT)
825 : :
826 : 560818 : blockdim = lobpcg%blockdim
827 : 560818 : spacedim = lobpcg%spacedim
828 : 560818 : previousBlock = (iblock-1)*lobpcg%blockdim
829 : :
830 : 560818 : space_buf = space(var)
831 [ + + ]: 560818 : if (space(var)==SPACE_CR) then
832 : 22654 : space_buf = SPACE_R
833 : : end if
834 : 560818 : call xg_init(buffer,space_buf,previousBlock,blockdim,comm=lobpcg%spacecom,gpu_option=lobpcg%gpu_option)
835 : :
836 : : ! buffer = BX0^T*X
837 : 560818 : call xgBlock_gemm('t','n',1.0d0,lobpcg%BX0,var,0.d0,buffer%self,comm=lobpcg%spacecom)
838 : :
839 : : ! sum all process contribution
840 : : ! X = - X0*(BX0^T*X) + X
841 : 560818 : call xgBlock_gemm('n','n',-1.0d0,lobpcg%X0,buffer%self,1.0d0,var)
842 : :
843 : 560818 : call xg_free(buffer)
844 : :
845 : : ABI_NVTX_END_RANGE()
846 : 560818 : call timab(tim_ortho,2,tsec)
847 : :
848 : 560818 : end subroutine lobpcg_orthoXwrtBlocks
849 : :
850 : 697325 : subroutine lobpcg_getResidu(lobpcg,eigenvalues)
851 : :
852 : : type(lobpcg_t) , intent(inout) :: lobpcg
853 : : type(xgBlock_t), intent(in ) :: eigenvalues
854 : : double precision :: tsec(2)
855 : :
856 : 697325 : call timab(tim_maxres,1,tsec)
857 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_RESIDUE)
858 : : !lobpcg%XWP(1:spacedim,shiftW+iblock) = lobpcg%AXWP(:,shiftX+iblock) - lobpcg%BXWP(:,shiftX+iblock)*eigenvalues(iblock)
859 : 697325 : call xgBlock_colwiseCymax(lobpcg%W,eigenvalues,lobpcg%BX,lobpcg%AX)
860 : : ABI_NVTX_END_RANGE()
861 : 697325 : call timab(tim_maxres,2,tsec)
862 : 697325 : end subroutine lobpcg_getResidu
863 : :
864 : 139855 : subroutine lobpcg_setX0(lobpcg,iblock)
865 : :
866 : : type(lobpcg_t) , intent(inout) :: lobpcg
867 : : integer , intent(in ) :: iblock
868 : : type(xgBlock_t) :: Xtmp
869 : : integer :: blockdim
870 : : integer :: spacedim
871 : : double precision :: tsec(2)
872 : :
873 : 139855 : call timab(tim_copy,1,tsec)
874 : 139855 : blockdim = lobpcg%blockdim
875 : 139855 : spacedim = lobpcg%spacedim
876 : :
877 : : !X0(:,(iblock-1)*blockdim+1:iblock*blockdim) = lobpcg%XWP(:,lobpcg%X+1:lobpcg%X+blockdim)
878 : 139855 : call xgBlock_setBlock(lobpcg%AllX0,Xtmp,spacedim,blockdim,fcol=(iblock-1)*blockdim+1)
879 : 139855 : call xgBlock_copy(lobpcg%X,Xtmp)
880 : 139855 : call timab(tim_copy,2,tsec)
881 : :
882 : 139855 : end subroutine lobpcg_setX0
883 : :
884 : :
885 : 133478 : subroutine lobpcg_transferAX_BX(lobpcg,iblock)
886 : :
887 : : type(lobpcg_t), intent(inout) :: lobpcg
888 : : integer , intent(in ) :: iblock
889 : : type(xgBlock_t) :: CXtmp
890 : : integer :: firstcol
891 : : double precision :: tsec(2)
892 : :
893 : 133478 : call timab(tim_copy,1,tsec)
894 : :
895 [ - + ]: 133478 : if (iblock<1) then
896 : 0 : ABI_ERROR("iblock<1")
897 : : end if
898 : :
899 : : ! iblock goes from 1 to nblock-1 included
900 : 133478 : firstcol = (iblock-1)*lobpcg%blockdim+1 ! Start of each block
901 : :
902 : : ! BX
903 : 133478 : call xg_setBlock(lobpcg%AllBX0,CXtmp,lobpcg%spacedim,lobpcg%blockdim,fcol=firstcol)
904 : 133478 : call xgBlock_copy(lobpcg%BX,CXtmp)
905 : :
906 : : ! AX
907 : 133478 : call xg_setBlock(lobpcg%AllAX0,CXtmp,lobpcg%spacedim,lobpcg%blockdim,fcol=firstcol)
908 : 133478 : call xgBlock_copy(lobpcg%AX,CXtmp)
909 : :
910 : 133478 : call timab(tim_copy,2,tsec)
911 : :
912 : 133478 : end subroutine lobpcg_transferAX_BX
913 : :
914 : 48376 : subroutine lobpcg_free(lobpcg)
915 : :
916 : : type(lobpcg_t), intent(inout) :: lobpcg
917 : :
918 : 48376 : call xg_free(lobpcg%XWP)
919 : 48376 : call xg_free(lobpcg%AXWP)
920 : 48376 : call xg_free(lobpcg%BXWP)
921 : 48376 : call xg_free(lobpcg%AllAX0)
922 : 48376 : call xg_free(lobpcg%AllBX0)
923 : 48376 : end subroutine lobpcg_free
924 : :
925 [ # # # # : 0 : end module m_lobpcg2
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
926 : : !!***
|
