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