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 24063 : 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 24063 : call timab(tim_init,1,tsec)
198 24063 : lobpcg%neigenpairs = neigenpairs
199 24063 : lobpcg%spacedim = spacedim
200 24063 : lobpcg%blockdim = blockdim
201 24063 : if (tolerance > 0.0) then
202 4568 : lobpcg%tolerance = tolerance
203 : else
204 19495 : lobpcg%tolerance = 1.0e-20
205 : end if
206 24063 : lobpcg%nline = nline
207 24063 : lobpcg%spacecom = spacecom
208 24063 : lobpcg%nblock = neigenpairs / blockdim
209 24063 : lobpcg%paral_kgb = paral_kgb
210 24063 : lobpcg%comm_rows = comm_rows
211 24063 : lobpcg%comm_cols = comm_cols
212 24063 : lobpcg%me_g0 = me_g0
213 24063 : lobpcg%me_g0_fft = me_g0_fft
214 24063 : lobpcg%gpu_option = gpu_option
215 24063 : lobpcg%gpu_thread_limit = 0
216 24063 : if(present(gpu_thread_limit)) lobpcg%gpu_thread_limit = gpu_thread_limit
217 :
218 24063 : call lobpcg_allocateAll(lobpcg,space,me_g0)
219 24063 : call timab(tim_init,2,tsec)
220 :
221 24063 : end subroutine lobpcg_init
222 :
223 :
224 24063 : 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 24063 : spacedim = lobpcg%spacedim
233 24063 : blockdim = lobpcg%blockdim
234 :
235 24063 : call lobpcg_free(lobpcg) ! Make sure everything is not allocated and
236 : ! pointer point to null()
237 :
238 24063 : call xg_init(lobpcg%XWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
239 24063 : call xg_setBlock(lobpcg%XWP,lobpcg%X,spacedim,blockdim)
240 24063 : call xg_setBlock(lobpcg%XWP,lobpcg%W,spacedim,blockdim,fcol=blockdim+1)
241 24063 : call xg_setBlock(lobpcg%XWP,lobpcg%P,spacedim,blockdim,fcol=2*blockdim+1)
242 24063 : call xg_setBlock(lobpcg%XWP,lobpcg%XW,spacedim,2*blockdim)
243 24063 : call xg_setBlock(lobpcg%XWP,lobpcg%WP,spacedim,2*blockdim,fcol=blockdim+1)
244 :
245 24063 : call xg_init(lobpcg%AXWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
246 24063 : call xg_setBlock(lobpcg%AXWP,lobpcg%AX,spacedim,blockdim)
247 24063 : call xg_setBlock(lobpcg%AXWP,lobpcg%AW,spacedim,blockdim,fcol=blockdim+1)
248 24063 : call xg_setBlock(lobpcg%AXWP,lobpcg%AP,spacedim,blockdim,fcol=2*blockdim+1)
249 24063 : call xg_setBlock(lobpcg%AXWP,lobpcg%AXW,spacedim,2*blockdim)
250 24063 : call xg_setBlock(lobpcg%AXWP,lobpcg%AWP,spacedim,2*blockdim,fcol=blockdim+1)
251 :
252 24063 : call xg_init(lobpcg%BXWP,space,spacedim,3*blockdim,lobpcg%spacecom,me_g0=me_g0,gpu_option=lobpcg%gpu_option)
253 24063 : call xg_setBlock(lobpcg%BXWP,lobpcg%BX,spacedim,blockdim)
254 24063 : call xg_setBlock(lobpcg%BXWP,lobpcg%BW,spacedim,blockdim,fcol=blockdim+1)
255 24063 : call xg_setBlock(lobpcg%BXWP,lobpcg%BP,spacedim,blockdim,fcol=2*blockdim+1)
256 24063 : call xg_setBlock(lobpcg%BXWP,lobpcg%BXW,spacedim,2*blockdim)
257 24063 : call xg_setBlock(lobpcg%BXWP,lobpcg%BWP,spacedim,2*blockdim,fcol=blockdim+1)
258 :
259 24063 : 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 6252 : lobpcg%AllBX0%self = lobpcg%BX
264 6252 : lobpcg%AllAX0%self = lobpcg%AX
265 : end if
266 :
267 :
268 24063 : 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 72189 : 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 24063 : lobpcg%prtvol = prtvol
383 :
384 48126 : tolerance=2*dlamch('E')
385 :
386 24063 : blockdim = lobpcg%blockdim
387 24063 : blockdim2 = 2*blockdim
388 24063 : blockdim3 = 3*blockdim
389 24063 : spacedim = lobpcg%spacedim
390 :
391 24063 : nblock = lobpcg%nblock
392 24063 : nline = lobpcg%nline
393 :
394 24063 : if (nbdbuf>0) then
395 7392 : nband_eff = lobpcg%neigenpairs - nbdbuf
396 : else
397 16671 : nband_eff = lobpcg%neigenpairs
398 : end if
399 :
400 24063 : call xgBlock_getSize(eigen,rows_tmp, cols_tmp)
401 24063 : if ( rows_tmp /= lobpcg%neigenpairs .and. cols_tmp /= 1 ) then
402 0 : ABI_ERROR("Error eigen size")
403 : endif
404 24063 : call xgBlock_getSize(X0,rows_tmp, cols_tmp)
405 24063 : if ( rows_tmp /= lobpcg%spacedim ) then
406 0 : ABI_ERROR("Error X0 spacedim")
407 : endif
408 24063 : if ( cols_tmp /= lobpcg%neigenpairs ) then
409 0 : ABI_ERROR("Error X0 npairs")
410 : endif
411 :
412 24063 : if (isppol==1.and.ikpt==1.and.inonsc==1.and.istep==1) then
413 687 : write(msg,'(a,es16.6)') ' lobpcg%tolerance(tolwfr_diago)=',lobpcg%tolerance
414 687 : call wrtout(std_out,msg,'COLL')
415 : end if
416 :
417 24063 : call xg_init(eigenvalues3N,SPACE_R,blockdim3,1, gpu_option=lobpcg%gpu_option)
418 24063 : call xg_setBlock(eigenvalues3N,eigenvaluesN,blockdim,1)
419 24063 : call xg_setBlock(eigenvalues3N,eigenvalues2N,blockdim2,1)
420 :
421 24063 : call xgBlock_reshape(eigen,blockdim,nblock)
422 24063 : call xgBlock_reshape(residu,blockdim,nblock)
423 24063 : call xgBlock_reshape(occ,blockdim,nblock)
424 :
425 24063 : lobpcg%AllX0 = X0
426 :
427 24063 : call xg_init(residu_eff,SPACE_R,blockdim,1,gpu_option=ABI_GPU_DISABLED)
428 :
429 24063 : if ( lobpcg%paral_kgb == 1 ) then
430 21564 : 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 21564 : gpu_option=lobpcg%gpu_option,gpu_thread_limit=lobpcg%gpu_thread_limit)
434 : call xgTransposer_copyConstructor(lobpcg%xgTransposerAX,lobpcg%xgTransposerX,&
435 21564 : lobpcg%AX,lobpcg%AXColsRows,STATE_LINALG)
436 : call xgTransposer_copyConstructor(lobpcg%xgTransposerBX,lobpcg%xgTransposerX,&
437 21564 : lobpcg%BX,lobpcg%BXColsRows,STATE_LINALG)
438 :
439 : call xgTransposer_copyConstructor(lobpcg%xgTransposerW,lobpcg%xgTransposerX,&
440 21564 : lobpcg%W,lobpcg%WColsRows,STATE_LINALG)
441 : call xgTransposer_copyConstructor(lobpcg%xgTransposerAW,lobpcg%xgTransposerX,&
442 21564 : lobpcg%AW,lobpcg%AWColsRows,STATE_LINALG)
443 : call xgTransposer_copyConstructor(lobpcg%xgTransposerBW,lobpcg%xgTransposerX,&
444 21564 : lobpcg%BW,lobpcg%BWColsRows,STATE_LINALG)
445 21564 : 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 163793 : do iblock = 1, nblock
457 : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_BLOCK)
458 139730 : nrestart = 0
459 :
460 139730 : call lobpcg_getX0(lobpcg,iblock)
461 139730 : call xgBlock_setBlock(residu,residuBlock,blockdim,1,fcol=iblock)
462 139730 : call xgBlock_setBlock(occ, occBlock, blockdim,1,fcol=iblock)
463 :
464 139730 : 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 139730 : if (lobpcg%paral_kgb == 1) then
472 132392 : call timab(tim_transpose,1,tsec)
473 132392 : call xgTransposer_transpose(lobpcg%xgTransposerX,STATE_COLSROWS)
474 132392 : lobpcg%xgTransposerAX%state=STATE_COLSROWS
475 132392 : lobpcg%xgTransposerBX%state=STATE_COLSROWS
476 132392 : call timab(tim_transpose,2,tsec)
477 : end if
478 : ! Initialize some quantitites (AX and BX)
479 139730 : call timab(tim_ax_bx,1,tsec)
480 139730 : call getAX_BX(lobpcg%XColsRows,lobpcg%AXColsRows,lobpcg%BXColsRows)
481 139730 : call xgBlock_zero_im_g0(lobpcg%AXColsRows)
482 139730 : call xgBlock_zero_im_g0(lobpcg%BXColsRows)
483 139730 : call timab(tim_ax_bx,2,tsec)
484 139730 : if (lobpcg%paral_kgb == 1) then
485 132392 : call timab(tim_transpose,1,tsec)
486 132392 : call xgTransposer_transpose(lobpcg%xgTransposerX,STATE_LINALG)
487 132392 : call xgTransposer_transpose(lobpcg%xgTransposerAX,STATE_LINALG)
488 132392 : call xgTransposer_transpose(lobpcg%xgTransposerBX,STATE_LINALG)
489 132392 : call timab(tim_transpose,2,tsec)
490 : end if
491 :
492 : ! B-orthonormalize X, BX and AX
493 139730 : 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 139730 : call xg_RayleighRitz(lobpcg%X,lobpcg%AX,lobpcg%BX,eigenvaluesN,ierr,lobpcg%prtvol,tim_RR_X,lobpcg%gpu_option)
497 :
498 139730 : compute_residu = .true.
499 :
500 696724 : do iline = 1, nline
501 : ABI_NVTX_START_RANGE(NVTX_LOBPCG2_LINE)
502 :
503 567522 : 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 567522 : call lobpcg_getResidu(lobpcg,eigenvaluesN)
512 :
513 : ! Compute residu norm here !
514 567522 : call timab(tim_maxres,1,tsec)
515 567522 : call xgBlock_colwiseNorm2(lobpcg%W,residuBlock)
516 567522 : call timab(tim_maxres,2,tsec)
517 :
518 : ! Apply preconditioner
519 567522 : call timab(tim_pcond,1,tsec)
520 567522 : call xgBlock_apply_diag(lobpcg%W,pcond,nspinor)
521 567522 : call timab(tim_pcond,2,tsec)
522 :
523 567522 : call timab(tim_nbdbuf,1,tsec)
524 567522 : if (nbdbuf>=0) then
525 : ! There is a transfer from GPU to CPU in this copy
526 566154 : call xgBlock_copy(residuBlock,residu_eff%self)
527 566154 : iband_min = 1 + blockdim*(iblock-1)
528 566154 : iband_max = blockdim*iblock
529 566154 : if (iband_max<=nband_eff) then ! all bands of this block are below nband_eff
530 536865 : 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 567522 : call timab(tim_nbdbuf,2,tsec)
546 567522 : 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 556994 : if ( iblock > 1 ) then
554 445151 : call lobpcg_orthoXwrtBlocks(lobpcg,lobpcg%W,iblock)
555 : end if
556 :
557 556994 : if (lobpcg%paral_kgb == 1) then
558 526296 : call timab(tim_transpose,1,tsec)
559 526296 : call xgTransposer_transpose(lobpcg%xgTransposerW,STATE_COLSROWS)
560 526296 : lobpcg%xgTransposerAW%state=STATE_COLSROWS
561 526296 : lobpcg%xgTransposerBW%state=STATE_COLSROWS
562 526296 : call timab(tim_transpose,2,tsec)
563 : end if
564 : ! Apply A and B on W
565 556994 : call timab(tim_ax_bx,1,tsec)
566 556994 : call getAX_BX(lobpcg%WColsRows,lobpcg%AWColsRows,lobpcg%BWColsRows)
567 556994 : call xgBlock_zero_im_g0(lobpcg%AWColsRows)
568 556994 : call xgBlock_zero_im_g0(lobpcg%BWColsRows)
569 556994 : call timab(tim_ax_bx,2,tsec)
570 556994 : if (lobpcg%paral_kgb == 1) then
571 526296 : call timab(tim_transpose,1,tsec)
572 526296 : call xgTransposer_transpose(lobpcg%xgTransposerW,STATE_LINALG)
573 526296 : call xgTransposer_transpose(lobpcg%xgTransposerAW,STATE_LINALG)
574 526296 : call xgTransposer_transpose(lobpcg%xgTransposerBW,STATE_LINALG)
575 526296 : 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 686196 : if ( iline == 1 .or. minResidu < 1e-27) then
582 : ! Do RR on XW to get the eigen vectors
583 133894 : call xg_Borthonormalize(lobpcg%XW,lobpcg%BXW,ierr,tim_Bortho_XW,lobpcg%gpu_option,AX=lobpcg%AXW) ! Do rotate AW
584 133894 : call xgBlock_zero(lobpcg%P)
585 133894 : call xgBlock_zero(lobpcg%AP)
586 133894 : call xgBlock_zero(lobpcg%BP)
587 133894 : 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 133894 : & AWP=lobpcg%AWP,BWP=lobpcg%BWP)
594 133894 : 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 423100 : 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 423100 : 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 423092 : & AWP=lobpcg%AWP,BWP=lobpcg%BWP,XWP=lobpcg%XWP%self)
607 423092 : 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 139730 : if ( compute_residu ) then
638 : ! Recompute AX-Lambda*BX for the last time
639 129202 : call lobpcg_getResidu(lobpcg,eigenvaluesN)
640 : ! Recompute residu norm here !
641 129202 : call timab(tim_maxres,1,tsec)
642 129202 : call xgBlock_colwiseNorm2(lobpcg%W,residuBlock)
643 129202 : call timab(tim_maxres,2,tsec)
644 : ! Apply preconditioner
645 129202 : call timab(tim_pcond,1,tsec)
646 129202 : call xgBlock_apply_diag(lobpcg%W,pcond,nspinor)
647 129202 : call timab(tim_pcond,2,tsec)
648 :
649 129202 : call timab(tim_nbdbuf,1,tsec)
650 129202 : if(lobpcg%gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(residuBlock)
651 129202 : if (nbdbuf>=0) then
652 129047 : call xgBlock_copy(residuBlock,residu_eff%self)
653 129047 : iband_min = 1 + blockdim*(iblock-1)
654 129047 : iband_max = blockdim*iblock
655 129047 : if (iband_max<=nband_eff) then ! all bands of this block are below nband_eff
656 126246 : 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 129202 : call timab(tim_nbdbuf,2,tsec)
672 : end if
673 :
674 139730 : 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 139730 : call timab(tim_copy,1,tsec)
682 139730 : call xgBlock_setBlock(eigen,eigenBlock,blockdim,1,fcol=iblock)
683 139730 : call xgBlock_copy(eigenvaluesN,eigenBlock)
684 139730 : call timab(tim_copy,2,tsec)
685 :
686 : ! Save new X in X0
687 139730 : call lobpcg_setX0(lobpcg,iblock)
688 :
689 : ! Copy previous BX into BX0 for previous block
690 163793 : 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 24063 : call xgBlock_reshape(eigen,blockdim*nblock,1)
698 24063 : call xgBlock_reshape(residu,blockdim*nblock,1)
699 24063 : call xgBlock_reshape(occ,blockdim*nblock,1)
700 :
701 24063 : call xg_free(eigenvalues3N)
702 24063 : call xg_free(residu_eff)
703 :
704 24063 : skip = .false.
705 24063 : 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 24063 : 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 24063 : if ( lobpcg%paral_kgb == 1 ) then
755 21564 : call xgTransposer_free(lobpcg%xgTransposerX)
756 21564 : call xgTransposer_free(lobpcg%xgTransposerAX)
757 21564 : call xgTransposer_free(lobpcg%xgTransposerBX)
758 21564 : call xgTransposer_free(lobpcg%xgTransposerW)
759 21564 : call xgTransposer_free(lobpcg%xgTransposerAW)
760 21564 : call xgTransposer_free(lobpcg%xgTransposerBW)
761 : end if
762 :
763 : ! call timab(tim_run,2,tsec)
764 :
765 24063 : end subroutine lobpcg_run
766 :
767 :
768 139730 : 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 139730 : call timab(tim_copy,1,tsec)
777 :
778 139730 : blockdim = lobpcg%blockdim
779 139730 : spacedim = lobpcg%spacedim
780 :
781 : !lobpcg%XWP(:,X+1:X+blockdim) = lobpcg%X0(:,(iblock-1)*blockdim+1:iblock*blockdim)
782 139730 : call xgBlock_setBlock(lobpcg%AllX0,lobpcg%X0,spacedim,blockdim,fcol=(iblock-1)*blockdim+1)
783 139730 : call xgBlock_copy(lobpcg%X0,lobpcg%X)
784 :
785 139730 : call timab(tim_copy,2,tsec)
786 :
787 139730 : 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 696724 : 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 696724 : 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 696724 : call xgBlock_colwiseCymax(lobpcg%W,eigenvalues,lobpcg%BX,lobpcg%AX)
852 : ABI_NVTX_END_RANGE()
853 696724 : call timab(tim_maxres,2,tsec)
854 696724 : end subroutine lobpcg_getResidu
855 :
856 139730 : 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 139730 : call timab(tim_copy,1,tsec)
866 139730 : blockdim = lobpcg%blockdim
867 139730 : spacedim = lobpcg%spacedim
868 :
869 : !X0(:,(iblock-1)*blockdim+1:iblock*blockdim) = lobpcg%XWP(:,lobpcg%X+1:lobpcg%X+blockdim)
870 139730 : call xgBlock_setBlock(lobpcg%AllX0,Xtmp,spacedim,blockdim,fcol=(iblock-1)*blockdim+1)
871 139730 : call xgBlock_copy(lobpcg%X,Xtmp)
872 139730 : call timab(tim_copy,2,tsec)
873 :
874 139730 : 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 48126 : subroutine lobpcg_free(lobpcg)
907 :
908 : type(lobpcg_t), intent(inout) :: lobpcg
909 :
910 48126 : call xg_free(lobpcg%XWP)
911 48126 : call xg_free(lobpcg%AXWP)
912 48126 : call xg_free(lobpcg%BXWP)
913 48126 : call xg_free(lobpcg%AllAX0)
914 48126 : call xg_free(lobpcg%AllBX0)
915 48126 : end subroutine lobpcg_free
916 :
917 0 : end module m_lobpcg2
918 : !!***
|
