Branch data Line data Source code
1 : : !!****m* ABINIT/m_xg_ortho_RR
2 : : !! NAME
3 : : !! m_xg_ortho_RR
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2024-2026 ABINIT group (J. Bieder, L. Baguet)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! NOTES
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : #include "nvtx_macros.h"
24 : :
25 : : module m_xg_ortho_RR
26 : :
27 : : use m_errors
28 : : use m_abicore
29 : : use defs_basis
30 : : use m_time, only : timab,abi_wtime
31 : : use m_xmpi
32 : :
33 : : use m_xg
34 : : use m_xg_nonlop
35 : : use m_xgScalapack
36 : :
37 : : #if defined(HAVE_GPU)
38 : : use m_gpu_toolbox
39 : : #endif
40 : :
41 : : #if defined(HAVE_GPU_MARKERS)
42 : : use m_nvtx_data
43 : : #endif
44 : :
45 : : implicit none
46 : :
47 : : private
48 : :
49 : : integer, parameter :: VAR_X = 1000
50 : : integer, parameter :: VAR_XW = 1010
51 : : integer, parameter :: VAR_XWP = 1100
52 : :
53 : : integer, parameter :: EIGENVX = 1
54 : : integer, parameter :: EIGENVD = 2
55 : : integer, parameter :: EIGENV = 3
56 : : integer, parameter :: EIGENPVX = 4
57 : : integer, parameter :: EIGENPVD = 5
58 : : integer, parameter :: EIGENPV = 6
59 : : integer, parameter :: EIGENEVD = 7
60 : : integer, parameter :: EIGENEV = 8
61 : : integer, parameter :: EIGENPEVD = 9
62 : : integer, parameter :: EIGENPEV = 10
63 : : integer, parameter :: EIGENSLK = 11
64 : : integer, parameter :: EIGENDROP = 12
65 : : logical, parameter :: EIGPACK(12) = &
66 : : (/ .false.,.false.,.false., &
67 : : .true. ,.true. ,.true. ,&
68 : : .false.,.false.,&
69 : : .true. ,.true., .false., .false. /)
70 : :
71 : : integer, parameter :: tim_RR_diago = 1795
72 : : integer, parameter :: tim_RR_gemm_1 = 1796
73 : : integer, parameter :: tim_RR_gemm_2 = 1797
74 : :
75 : : public :: xg_Borthonormalize
76 : : public :: xg_Borthonormalize_cprj
77 : : public :: xg_RayleighRitz
78 : : public :: xg_RayleighRitz_cprj
79 : : ! Number of regularized Gram matrices; reset by the caller (scfcv_core, before each vtorho)
80 : : integer,public,save :: xg_bortho_reg_counter = 0
81 : :
82 : : contains
83 : : !!***
84 : :
85 : : !!****f* m_xg_ortho_RR/xg_Borthonormalize
86 : : !!
87 : : !! NAME
88 : : !! xg_Borthonormalize
89 : 1430310 : subroutine xg_Borthonormalize(X,BX,info,timer,gpu_option,AX,allow_regularize,regularized)
90 : :
91 : : integer , intent(in ) :: timer
92 : : integer , intent(in ) :: gpu_option
93 : : type(xgBlock_t), intent(inout) :: X
94 : : type(xgBlock_t), intent(inout) :: BX
95 : : type(xgBlock_t), intent(inout),optional :: AX
96 : : logical , intent(in ),optional :: allow_regularize
97 : : logical , intent(out ),optional :: regularized
98 : : integer , intent(out ) :: info
99 : : logical :: do_regularize
100 : : integer :: space_buf,k,nn,ierr,nreg
101 : : type(xg_t) :: buffer,evals,evec
102 : : double precision :: tol
103 : : double precision :: tsec(2)
104 : 715155 : double precision,pointer :: ev(:,:), vmat(:,:)
105 : : character(len=500) :: msg
106 : :
107 : : ! Initialize optional arguments
108 : 715155 : do_regularize=.false.
109 : 134038 : if (present(allow_regularize)) do_regularize=allow_regularize
110 [ + + ]: 715155 : if (present(regularized)) regularized=.false.
111 : :
112 : : ABI_NVTX_START_RANGE(NVTX_B_ORTHO)
113 : :
114 : 715155 : call timab(timer,1,tsec)
115 : :
116 : 715155 : call xgBlock_check(X,BX)
117 [ + - ]: 715155 : if (present(AX)) then
118 : 715155 : call xgBlock_check(X,AX)
119 : : end if
120 : :
121 : 715155 : nn = cols(X)
122 [ + + ]: 715155 : if (space(X)/=SPACE_CR) then
123 : 658368 : space_buf = SPACE(X)
124 : : else
125 : 56787 : space_buf = SPACE_R
126 : : end if
127 : 715155 : call xg_init(buffer,space_buf,nn,nn,comm(X),gpu_option=gpu_option)
128 : :
129 : : ! If space(X)==SPACE_CR : set imaginary part of G=0 component to zero to improve numerical stability
130 : 715155 : call xgBlock_zero_im_g0(X)
131 : 715155 : call xgBlock_zero_im_g0(BX)
132 [ + - ]: 715155 : if (present(AX)) then
133 : 715155 : call xgBlock_zero_im_g0(AX)
134 : : end if
135 : :
136 : : ! Compute X^TBX
137 : 715155 : call xgBlock_gemm('t','n',1.d0,X,BX,0.d0,buffer%self,comm=comm(X))
138 : :
139 : : ! Compute Cholesky decomposition (Upper part)
140 : 715155 : call xgBlock_potrf(buffer%self,'u',info)
141 : :
142 [ + + ]: 715155 : if ( info /= 0 ) then
143 [ + + ]: 25 : if (.not.do_regularize) then
144 : 19 : ABI_COMMENT("Cholesky decomposition did not work. Orthonormalization not done")
145 : 19 : call xg_free(buffer)
146 : : ABI_NVTX_END_RANGE()
147 : 19 : return
148 : : end if
149 : :
150 : : ! Cholesky failed. Only regularize when |lam_min| ~ 0 (Gram is singular but
151 : : ! valid). A large negative lam_min means the Gram is not valid, so do nothing.
152 : : ! Rebuild G = X^TBX since potrf overwrote the upper triangle.
153 : 6 : call xgBlock_gemm('t','n',1.d0,X,BX,0.d0,buffer%self,comm=comm(X))
154 : : ! Diagonalize G: evec <- eigenvectors V, evals <- eigenvalues lam (ascending)
155 : 6 : call xg_init(evals,SPACE_R,nn,1,comm(X),gpu_option=gpu_option)
156 : 6 : call xg_init(evec,space_buf,nn,nn,comm(X),gpu_option=gpu_option)
157 : 6 : call xgBlock_copy(buffer%self,evec%self)
158 : 6 : call xgBlock_heevd('v','u',evec%self,evals%self,ierr)
159 [ - + ]: 6 : if (ierr/=0) then
160 : 0 : ABI_COMMENT("Regularization failed. Orthonormalization not done")
161 : 0 : call xg_free(evals)
162 : 0 : call xg_free(evec)
163 : 0 : call xg_free(buffer)
164 : : ABI_NVTX_END_RANGE()
165 : 0 : return
166 : : end if
167 [ - + ]: 6 : if(gpu_option==ABI_GPU_OPENMP) then
168 : 0 : call xgBlock_copy_from_gpu(evals%self)
169 : 0 : call xgBlock_copy_from_gpu(evec%self)
170 : : end if
171 : 6 : call xgBlock_reverseMap(evals%self,ev)
172 : 6 : call xgBlock_reverseMap(evec%self,vmat)
173 : : ! Defining a threshold relatively to the matrix amplitude. The threshold must be
174 : : ! small compared to the highest eigenvalue.
175 : 6 : tol=tol12*ev(nn,1)
176 : : ! Relative threshold: regularize only if |lam_min| <= tol (round-off zero).
177 : : ! A NaN or infinite eigenvalue would pass a plain "abs(lam_min)>tol" test: refuse it explicitly.
178 [ + + + - : 204 : if (any(ev(1:nn,1)/=ev(1:nn,1)) .or. any(abs(ev(1:nn,1))>huge(one)) .or. abs(ev(1,1))>tol) then
+ + + - +
- - + ]
179 : 0 : write(msg,'(a,2(a,es10.3),a)') 'Cholesky failed and the Gram matrix cannot be regularized',&
180 : 0 : & ' (lam_min=',ev(1,1),', lam_max=',ev(nn,1),'). Orthonormalization not done'
181 : 0 : ABI_COMMENT(msg)
182 : 0 : call xg_free(evals)
183 : 0 : call xg_free(evec)
184 : 0 : call xg_free(buffer)
185 : : ABI_NVTX_END_RANGE()
186 : 0 : return
187 : : end if
188 : : ! Flooring the eigenvalues to the threshold and prepare
189 : : ! the vectors to rebuild the Gram matrix.
190 : 6 : nreg=0
191 [ + + ]: 102 : do k=1,nn
192 [ + + ]: 96 : if (ev(k,1)<tol) nreg=nreg+1
193 [ + + ]: 3174 : vmat(:,k)=vmat(:,k)*sqrt(max(ev(k,1),tol))
194 : : end do
195 : 6 : write(msg,'(a,i0,a,i0,3(a,es10.3))') 'Cholesky failed: regularizing the Gram matrix, n=',nn,&
196 : 12 : & ', floored=',nreg,', lam_min=',ev(1,1),', lam_max=',ev(nn,1),', tol=',tol
197 : 6 : ABI_COMMENT(msg)
198 [ - + ]: 6 : if(gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_to_gpu(evec%self)
199 : : ! Rebuild the regularized Gram = V.diag(max(lam,tol)).V^T.
200 : : ! V is the same on all processes: no MPI sum.
201 : 6 : call xgBlock_gemm('n','t',1.d0,evec%self,evec%self,0.d0,buffer%self)
202 : 6 : call xg_free(evec)
203 : 6 : call xg_free(evals)
204 : 6 : call xgBlock_potrf(buffer%self,'u',info)
205 [ - + ]: 6 : if (info/=0) then
206 : 0 : ABI_COMMENT("Cholesky failed on regularized Gram. Orthonormalization not done")
207 : 0 : call xg_free(buffer)
208 : : ABI_NVTX_END_RANGE()
209 : 0 : return
210 : : end if
211 : 6 : xg_bortho_reg_counter=xg_bortho_reg_counter+1
212 [ + - ]: 6 : if (present(regularized)) regularized=.true.
213 : : end if
214 : :
215 : : ! Solve YU=X
216 : 715136 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,X)
217 : :
218 : : ! Solve BYU=BX
219 : 715136 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,BX)
220 : :
221 [ + - ]: 715136 : if (present(AX)) then
222 : : ! Solve AYU=AX
223 : 715136 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,AX)
224 : : end if
225 : :
226 : 715136 : call xg_free(buffer)
227 : :
228 : : ABI_NVTX_END_RANGE()
229 : 715136 : call timab(timer,2,tsec)
230 : :
231 [ + + ]: 715155 : end subroutine xg_Borthonormalize
232 : : !!***
233 : :
234 : : !!****f* m_xg_ortho_RR/xg_Borthonormalize_cprj
235 : : !!
236 : : !! NAME
237 : : !! xg_Borthonormalize_cprj
238 : 158372 : subroutine xg_Borthonormalize_cprj(xg_nonlop,X,cprjX,info,timer,gpu_option,&
239 : : & blockdim_cprj,AX,allow_regularize,regularized)
240 : :
241 : : integer , intent(in ) :: timer,gpu_option
242 : : integer , intent(out ) :: info
243 : : type(xg_nonlop_t), intent(in ) :: xg_nonlop
244 : : type(xgBlock_t) , intent(inout) :: X
245 : : type(xgBlock_t) , intent(inout) :: cprjX
246 : : integer , intent(in ),optional :: blockdim_cprj
247 : : type(xgBlock_t) , intent(inout),optional :: AX
248 : : logical , intent(in ),optional :: allow_regularize
249 : : logical , intent( out),optional :: regularized
250 : : logical :: do_regularize
251 : : type(xg_t) :: buffer,cprj_work,evals,evec
252 : : type(xgBlock_t) :: cprjX_spinor,cprj_work_spinor
253 : : integer :: blockdim_cprj_,space_buf,k,ierr,nreg
254 : : integer :: spacecom,ncols_cprj,nn,nspinor
255 : : double precision :: tol
256 : : double precision :: tsec(2)
257 : 79186 : double precision,pointer :: ev(:,:), vmat(:,:)
258 : : character(len=500) :: msg
259 : :
260 : : ! Initialize optional arguments
261 : 79186 : do_regularize=.false.
262 : 8968 : if (present(allow_regularize)) do_regularize=allow_regularize
263 [ + + ]: 79186 : if (present(regularized)) regularized=.false.
264 : :
265 : 79186 : call timab(timer,1,tsec)
266 : :
267 [ - + ]: 79186 : if (gpu_option/=ABI_GPU_DISABLED) then
268 : 0 : ABI_ERROR('Not implemented for GPU')
269 : : end if
270 : 79186 : call xgBlock_check_gpu_option(X,cprjX)
271 [ + - ]: 79186 : if (present(AX)) then
272 : 79186 : call xgBlock_check_gpu_option(X,AX)
273 : : end if
274 : :
275 : 79186 : nn = cols(X)
276 : 79186 : ncols_cprj = cols(cprjX)
277 : 79186 : blockdim_cprj_ = ncols_cprj
278 [ + + ]: 79186 : if (present(blockdim_cprj)) then
279 : 58573 : blockdim_cprj_ = blockdim_cprj
280 : : end if
281 : 79186 : nspinor = xg_nonlop%nspinor
282 : :
283 : 79186 : spacecom = comm(X)
284 : :
285 [ + + ]: 79186 : if (space(X)/=SPACE_CR) then
286 : 40917 : space_buf = SPACE(X)
287 : : else
288 : 38269 : space_buf = SPACE_R
289 : : end if
290 : 79186 : call xg_init(buffer,space_buf,nn,nn,spacecom)
291 : :
292 : : ! Compute X^TX
293 : 79186 : call xgBlock_gemm('t','n',1.d0,X,X,0.d0,buffer%self,comm=spacecom)
294 : 79186 : call xg_init(cprj_work,space(cprjX),rows(cprjX),ncols_cprj,spacecom)
295 [ + + ]: 79186 : if (xg_nonlop%paw) then
296 : 65499 : call xg_nonlop_getXSY(xg_nonlop,cprjX,cprjX,cprj_work%self,buffer%self,blocksize=blockdim_cprj_)
297 : : end if
298 : :
299 : : ! Compute Cholesky decomposition (Upper part)
300 : 79186 : call xgBlock_potrf(buffer%self,'u',info)
301 : :
302 [ + + ]: 79186 : if ( info /= 0 ) then
303 [ + + ]: 21 : if (.not.do_regularize) then
304 : 13 : ABI_COMMENT("Cholesky decomposition did not work. Orthonormalization not done")
305 : 13 : call xg_free(cprj_work)
306 : 13 : call xg_free(buffer)
307 : 13 : return
308 : : end if
309 : :
310 : : ! Cholesky failed. Only regularize when |lam_min| ~ 0 (Gram is singular but
311 : : ! valid). A large negative lam_min means the Gram is not valid, so do nothing.
312 : : ! Rebuild G = X^T.X + cprj^H.Sij.cprj since potrf overwrote the upper triangle.
313 : 8 : call xgBlock_gemm('t','n',1.d0,X,X,0.d0,buffer%self,comm=spacecom)
314 [ + - ]: 8 : if (xg_nonlop%paw) then
315 : 8 : call xg_nonlop_getXSY(xg_nonlop,cprjX,cprjX,cprj_work%self,buffer%self,blocksize=blockdim_cprj_)
316 : : end if
317 : : ! Diagonalize G: evec <- eigenvectors V, evals <- eigenvalues lam (ascending)
318 : 8 : call xg_init(evals,SPACE_R,nn,1,spacecom)
319 : 8 : call xg_init(evec,space_buf,nn,nn,spacecom)
320 : 8 : call xgBlock_copy(buffer%self,evec%self)
321 : 8 : call xgBlock_heevd('v','u',evec%self,evals%self,ierr)
322 [ - + ]: 8 : if (ierr/=0) then
323 : 0 : ABI_COMMENT("Regularization failed. Orthonormalization not done")
324 : 0 : call xg_free(evals)
325 : 0 : call xg_free(evec)
326 : 0 : call xg_free(cprj_work)
327 : 0 : call xg_free(buffer)
328 : 0 : return
329 : : end if
330 : 8 : call xgBlock_reverseMap(evals%self,ev)
331 : 8 : call xgBlock_reverseMap(evec%self,vmat)
332 : : ! Defining a threshold relatively to the matrix amplitude. The threshold must be
333 : : ! small compared to the highest eigenvalue.
334 : 8 : tol=tol12*ev(nn,1)
335 : : ! Relative threshold: regularize only if |lam_min| <= tol (round-off zero).
336 : : ! A NaN or infinite eigenvalue would pass a plain "abs(lam_min)>tol" test: refuse it explicitly.
337 [ + + + - : 272 : if (any(ev(1:nn,1)/=ev(1:nn,1)) .or. any(abs(ev(1:nn,1))>huge(one)) .or. abs(ev(1,1))>tol) then
+ + + - +
- - + ]
338 : 0 : write(msg,'(a,2(a,es10.3),a)') 'Cholesky failed and the Gram matrix cannot be regularized',&
339 : 0 : & ' (lam_min=',ev(1,1),', lam_max=',ev(nn,1),'). Orthonormalization not done'
340 : 0 : ABI_COMMENT(msg)
341 : 0 : call xg_free(evals)
342 : 0 : call xg_free(evec)
343 : 0 : call xg_free(cprj_work)
344 : 0 : call xg_free(buffer)
345 : 0 : return
346 : : end if
347 : : ! Flooring the eigenvalues to the threshold and prepare
348 : : ! the vectors to rebuild the Gram matrix.
349 : 8 : nreg=0
350 [ + + ]: 136 : do k=1,nn
351 [ + + ]: 128 : if (ev(k,1)<tol) nreg=nreg+1
352 [ + + ]: 4232 : vmat(:,k)=vmat(:,k)*sqrt(max(ev(k,1),tol))
353 : : end do
354 : 8 : write(msg,'(a,i0,a,i0,3(a,es10.3))') 'Cholesky failed: regularizing the Gram matrix, n=',nn,&
355 : 16 : & ', floored=',nreg,', lam_min=',ev(1,1),', lam_max=',ev(nn,1),', tol=',tol
356 : 8 : ABI_COMMENT(msg)
357 : : ! Rebuild the regularized Gram = V.diag(max(lam,tol)).V^T.
358 : : ! V is the same on all processes: no MPI sum.
359 : 8 : call xgBlock_gemm('n','t',1.d0,evec%self,evec%self,0.d0,buffer%self)
360 : 8 : call xg_free(evec)
361 : 8 : call xg_free(evals)
362 : 8 : call xgBlock_potrf(buffer%self,'u',info)
363 [ - + ]: 8 : if ( info /= 0 ) then
364 : 0 : ABI_COMMENT("Cholesky failed on regularized Gram. Orthonormalization not done")
365 : 0 : call xg_free(cprj_work)
366 : 0 : call xg_free(buffer)
367 : 0 : return
368 : : end if
369 : 8 : xg_bortho_reg_counter=xg_bortho_reg_counter+1
370 [ + - ]: 8 : if (present(regularized)) regularized=.true.
371 : : end if
372 : :
373 : : ! Solve YU=X
374 : 79173 : call xgBlock_trsm('r','u','n','n',1.d0,buffer%self,X)
375 : :
376 [ + - ]: 79173 : if (present(AX)) then
377 : : ! Solve AYU=AX
378 : 79173 : call xgBlock_trsm('r','u','n','n',1.d0,buffer%self,AX)
379 : : end if
380 : :
381 : : ! Solve (cprjY)U=(cprjX)
382 : 79173 : call xgBlock_reshape_spinor(cprjX,cprjX_spinor,nspinor,COLS2ROWS)
383 [ + + ]: 79173 : if (ncols_cprj==nspinor*nn) then
384 : 28091 : call xgBlock_trsm('r','u','n','n',1.d0,buffer%self,cprjX_spinor)
385 : : else
386 : 51082 : call xgBlock_invert_tri('u','n',buffer%self)
387 : 51082 : call xgBlock_zerotri(buffer%self,'u')
388 : 51082 : call xgBlock_zero(cprj_work%self)
389 : 51082 : call xgBlock_reshape_spinor(cprj_work%self,cprj_work_spinor,nspinor,COLS2ROWS)
390 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjX_spinor,buffer%self,cprj_work_spinor,&
391 : 51082 : & xg_nonlop%me_band,blocksize=blockdim_cprj_/nspinor,comm=xg_nonlop%comm_band)
392 : 51082 : call xgBlock_copy(cprj_work%self,cprjX)
393 : : end if
394 : :
395 : 79173 : call xg_free(cprj_work)
396 : :
397 : 79173 : call xg_free(buffer)
398 : :
399 : 79173 : call timab(timer,2,tsec)
400 : :
401 [ + + ]: 79186 : end subroutine xg_Borthonormalize_cprj
402 : : !!***
403 : :
404 : : !!****f* m_xg_ortho_RR/xg_RayleighRitz
405 : : !!
406 : : !! NAME
407 : : !! xg_RayleighRitz
408 : 726392 : subroutine xg_RayleighRitz(X,AX,BX,eigenvalues,info,prtvol,timer,gpu_option,&
409 : : & tolerance,XW,AW,BW,P,AP,BP,WP,AWP,BWP,XWP,solve_ax_bx,drop_null)
410 : :
411 : : integer , intent(in ) :: timer
412 : : integer , intent(in ) :: gpu_option
413 : : integer , intent(in ) :: prtvol
414 : : type(xgBlock_t), intent(inout) :: eigenvalues
415 : : type(xgBlock_t), intent(inout) :: X
416 : : type(xgBlock_t), intent(inout) :: AX
417 : : type(xgBlock_t), intent(inout) :: BX
418 : : integer , intent( out) :: info
419 : : double precision, optional, intent(in) :: tolerance
420 : : ! LOBPCG only :
421 : : type(xgBlock_t), intent(inout),optional :: XW
422 : : type(xgBlock_t), intent(inout),optional :: AW
423 : : type(xgBlock_t), intent(inout),optional :: BW
424 : : type(xgBlock_t), intent(inout),optional :: P
425 : : type(xgBlock_t), intent(inout),optional :: AP
426 : : type(xgBlock_t), intent(inout),optional :: BP
427 : : type(xgBlock_t), intent(inout),optional :: WP
428 : : type(xgBlock_t), intent(inout),optional :: AWP
429 : : type(xgBlock_t), intent(inout),optional :: BWP
430 : : type(xgBlock_t), intent(inout),optional :: XWP
431 : : ! End LOBPCG only
432 : : logical, intent(in),optional :: solve_ax_bx
433 : : logical, intent(in),optional :: drop_null
434 : : integer :: var
435 : : integer :: spacedim
436 : : integer :: blockdim
437 : : integer :: space_buf
438 : : integer :: subdim
439 : : integer :: spacecom
440 : : integer :: eigenSolver
441 : : double precision :: abstol
442 : : #ifdef HAVE_LINALG_SCALAPACK
443 : : logical :: use_slk
444 : : #endif
445 : : type(xg_t) :: vec
446 : : type(xg_t) :: subA
447 : : type(xg_t) :: subB
448 : : type(xgBlock_t) :: subsub
449 : : type(xgBlock_t) :: Cwp
450 : : type(xgScalapack_t) :: scalapack
451 : : double precision :: tsec(2)
452 : : logical :: solve_ax_bx_,drop_null_
453 : :
454 : 726392 : call timab(timer , 1, tsec)
455 : :
456 : 726392 : solve_ax_bx_ = .false.
457 [ + + ]: 726392 : if (present(solve_ax_bx)) then
458 : 11256 : solve_ax_bx_ = solve_ax_bx
459 : : end if
460 : 726392 : drop_null_ = .false.
461 [ + + ]: 726392 : if (present(drop_null)) then
462 : 134038 : drop_null_ = drop_null
463 : : end if
464 : :
465 : 726392 : var = VAR_X
466 : 726392 : call xgBlock_check(X,AX)
467 : 726392 : call xgBlock_check(X,BX)
468 [ + + ]: 726392 : if (.not.solve_ax_bx_) then
469 : : eigenSolver = EIGENEVD
470 : : else
471 : 11256 : eigenSolver = EIGENVD
472 : : end if
473 : 726392 : spacedim = rows(X)
474 : 726392 : blockdim = cols(X)
475 : 726392 : space_buf = space(X)
476 [ + + ]: 726392 : if (space(X) == space_CR) then
477 : 60454 : space_buf = space_R
478 : : end if
479 : 726392 : spacecom = comm(X)
480 : 726392 : subdim = blockdim
481 : :
482 [ + + ]: 726392 : if (present(XW)) then
483 [ - + ]: 557470 : if (solve_ax_bx_) then
484 : 0 : ABI_ERROR('solve_ax_bx is not implemented in that case')
485 : : end if
486 : 557470 : var = VAR_XW
487 : 557470 : eigenSolver = EIGENVD
488 : 557470 : call xgBlock_check(X,AW)
489 : 557470 : call xgBlock_check(X,BW)
490 : 557470 : call xgBlock_check(X,P)
491 : 557470 : call xgBlock_check(X,AP)
492 : 557470 : call xgBlock_check(X,BP)
493 : 557470 : call xgBlock_check(X,XW,fact_col=2)
494 : 557470 : call xgBlock_check(X,WP,fact_col=2)
495 : 557470 : call xgBlock_check(X,AWP,fact_col=2)
496 : 557470 : call xgBlock_check(X,BWP,fact_col=2)
497 : 557470 : subdim = 2*blockdim
498 [ + + ]: 557470 : if (present(XWP)) then
499 : 423432 : var = VAR_XWP
500 : 423432 : call xgBlock_check(X,XWP,fact_col=3)
501 : 423432 : subdim = 3*blockdim
502 : : end if
503 : : end if
504 : :
505 : : #ifdef HAVE_LINALG_SCALAPACK
506 : : call xgScalapack_init(scalapack,spacecom,subdim,prtvol-2,(gpu_option/=ABI_GPU_DISABLED),use_slk)
507 : : if ( use_slk) then
508 : : eigenSolver = EIGENSLK
509 : : end if
510 : : #endif
511 : :
512 : : ! drop_null: use xg_hegv_drop (LAPACK on every rank, no ScaLAPACK), which drops the null directions of the Gram
513 : : ! matrix subB. subB only exists for the generalized problem, i.e. not for X alone without solve_ax_bx.
514 [ + + ]: 726392 : if (drop_null_) then
515 [ - + ]: 6 : if (var == VAR_X .and. .not.solve_ax_bx_) then
516 : 0 : ABI_BUG("drop_null needs the generalized eigenproblem")
517 : : end if
518 : : if ( eigenSolver == EIGENSLK ) then
519 : : call xgScalapack_free(scalapack)
520 : : end if
521 : : eigenSolver = EIGENDROP
522 : : end if
523 : :
524 : : ! Select diago algorithm
525 : :
526 [ + + ]: 726392 : abstol = 0d0 ; if ( present(tolerance) ) abstol = tolerance
527 : :
528 : 726392 : call xg_init(subA,space_buf,subdim,subdim,spacecom,gpu_option=gpu_option)
529 [ + + ]: 726392 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
530 : 568726 : call xg_init(subB,space_buf,subdim,subdim,spacecom,gpu_option=gpu_option)
531 : : end if
532 : :
533 : : if ( eigenSolver == EIGENVX .or. eigenSolver == EIGENPVX ) then
534 : : call xg_init(vec,space_buf,subdim,blockdim,gpu_option=gpu_option)
535 [ - + ]: 726392 : else if ( EIGPACK(eigenSolver) ) then
536 : 0 : call xg_init(vec,space_buf,subdim,subdim,gpu_option=gpu_option)
537 : : else
538 : 726392 : call xg_setBlock(subA,vec%self,subdim,blockdim)
539 : : endif
540 : :
541 : : ! Compute subA and subB by part
542 : : !--- begin
543 : : ! | E | XAW | XAP | | I | XBW | XBP |
544 : : ! | * | WAW | WAP | | * | I | WBP |
545 : : ! | * | * | PAP | | * | * | I |
546 : :
547 : 726392 : call timab(tim_RR_gemm_1,1,tsec)
548 : : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_1)
549 : :
550 : : ! If space(X)==SPACE_CR : set imaginary part of G=0 component to zero to improve numerical stability
551 [ + + ]: 726392 : if (var == VAR_X) call xgBlock_zero_im_g0(X)
552 [ + + ]: 726392 : if (var == VAR_XW) call xgBlock_zero_im_g0(XW)
553 [ + + ]: 726392 : if (var == VAR_XWP) call xgBlock_zero_im_g0(XWP)
554 : 726392 : call xgBlock_zero_im_g0(AX)
555 : 726392 : call xgBlock_zero_im_g0(BX)
556 : :
557 : 726392 : call xg_setBlock(subA,subsub,blockdim,blockdim)
558 : 726392 : call xgBlock_gemm('t','n',1.0d0,X,AX,0.d0,subsub,comm=spacecom)
559 : :
560 [ + + ]: 726392 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
561 : 568726 : call xg_setBlock(subB,subsub,blockdim,blockdim)
562 : 568726 : call xgBlock_gemm('t','n',1.0d0,X,BX,0.d0,subsub,comm=spacecom)
563 : : endif
564 : :
565 [ + + ]: 726392 : if ( var == VAR_XW .or. var == VAR_XWP ) then
566 : :
567 : : ! If space(X)==SPACE_CR : set imaginary part of G=0 component to zero to improve numerical stability
568 : 557470 : call xgBlock_zero_im_g0(AW)
569 : 557470 : call xgBlock_zero_im_g0(BW)
570 : :
571 : : ! subA
572 : 557470 : call xg_setBlock(subA,subsub,2*blockdim,blockdim,fcol=blockdim+1)
573 : 557470 : call xgBlock_gemm('t','n',1.0d0,XW,AW,0.d0,subsub,comm=spacecom)
574 : :
575 : : ! subB
576 : 557470 : call xg_setBlock(subB,subsub,2*blockdim,blockdim,fcol=blockdim+1)
577 : 557470 : call xgBlock_gemm('t','n',1.0d0,XW,BW,0.d0,subsub,comm=spacecom)
578 : :
579 : : end if
580 : :
581 [ + + ]: 557470 : if ( var == VAR_XWP ) then
582 : :
583 : : ! If space(X)==SPACE_CR : set imaginary part of G=0 component to zero to improve numerical stability
584 : 423432 : call xgBlock_zero_im_g0(AP)
585 : 423432 : call xgBlock_zero_im_g0(BP)
586 : :
587 : : ! subA
588 : 423432 : call xg_setBlock(subA,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
589 : 423432 : call xgBlock_gemm('t','n',1.0d0,XWP,AP,0.d0,subsub,comm=spacecom)
590 : :
591 : : ! subB
592 : 423432 : call xg_setBlock(subB,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
593 : 423432 : call xgBlock_gemm('t','n',1.0d0,XWP,BP,0.d0,subsub,comm=spacecom)
594 : :
595 : : end if
596 : :
597 : 726392 : call timab(tim_RR_gemm_1,2,tsec)
598 : : ABI_NVTX_END_RANGE()
599 : :
600 [ - + ]: 726392 : if ( EIGPACK(eigenSolver) ) then
601 : 0 : call xgBlock_pack(subA%self,subA%self,'u')
602 [ # # ]: 0 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
603 : 0 : call xgBlock_pack(subB%self,subB%self,'u')
604 : : end if
605 : : end if
606 : :
607 : 726392 : call timab(tim_RR_diago,1,tsec)
608 : 726392 : tsec(2) = abi_wtime()
609 [ + + ]: 726392 : if ( .not.solve_ax_bx_ .and. var == VAR_X ) then
610 : : ABI_NVTX_START_RANGE(NVTX_RR_HEEV)
611 : : ! Solve Hermitian eigen problem
612 : 157666 : select case (eigenSolver)
613 : : case (EIGENEVD)
614 [ + + ]: 157666 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heevd"
615 : 157666 : call xgBlock_heevd('v','u',subA%self,eigenvalues,info) ! work with GPU
616 : : case (EIGENEV)
617 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heev"
618 : 0 : call xgBlock_heev('v','u',subA%self,eigenvalues,info)
619 : : case (EIGENPEVD)
620 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpevd"
621 : 0 : call xgBlock_hpevd('v','u',subA%self,eigenvalues,vec%self,info)
622 : : case (EIGENPEV)
623 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpev"
624 : 0 : call xgBlock_hpev('v','u',subA%self,eigenvalues,vec%self,info)
625 : : case (EIGENSLK)
626 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using pheev"
627 : 0 : call xgScalapack_heev(scalapack,subA%self,eigenvalues,gpu_option=gpu_option) ! work with GPU
628 : 0 : info = 0 ! No error code returned for the moment
629 : : case default
630 [ + - - - : 315332 : ABI_ERROR("Error for Eigen Solver HEEV")
- - ]
631 : : end select
632 : : else
633 : : ABI_NVTX_START_RANGE(NVTX_RR_HEGV)
634 : : ! Solve Hermitian general eigen problem only for first blockdim eigenvalues
635 : : select case (eigenSolver)
636 : : case (EIGENVX)
637 : : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvx"
638 : : call xgBlock_hegvx(1,'v','i','u',subA%self,subB%self,0.d0,0.d0,1,blockdim,abstol,&
639 : 568720 : eigenvalues,vec%self,info)
640 : : case (EIGENVD)
641 [ + + ]: 568720 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvd"
642 : 568720 : call xgBlock_hegvd(1,'v','u',subA%self,subB%self,eigenvalues,info) ! work with GPU
643 : : case (EIGENV)
644 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegv"
645 : 0 : call xgBlock_hegv(1,'v','u',subA%self,subB%self,eigenvalues,info)
646 : : case (EIGENPVX)
647 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgvx"
648 : : call xgBlock_hpgvx(1,'v','i','u',subA%self,subB%self,0.d0,0.d0,1,blockdim,abstol,&
649 : 0 : eigenvalues,vec%self,info)
650 : : case (EIGENPVD)
651 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgvd"
652 : 0 : call xgBlock_hpgvd(1,'v','u',subA%self,subB%self,eigenvalues,vec%self,info)
653 : : case (EIGENPV)
654 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgv"
655 : 0 : call xgBlock_hpgv(1,'v','u',subA%self,subB%self,eigenvalues,vec%self,info)
656 : : case (EIGENSLK)
657 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using phegv"
658 : 0 : call xgScalapack_hegv(scalapack,subA%self,subB%self,eigenvalues,gpu_option=gpu_option) ! work with GPU
659 : 0 : info = 0 ! No error code returned for the moment
660 : : case (EIGENDROP)
661 [ - + ]: 6 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegv with dropping"
662 : 6 : call xg_hegv_drop(subA%self,subB%self,eigenvalues,blockdim,gpu_option,info) ! work with GPU
663 : : case default
664 [ + - - - : 568732 : ABI_ERROR("Error for Eigen Solver HEGV")
- - + - ]
665 : : end select
666 : : end if
667 : : if ( eigenSolver == EIGENSLK ) then
668 : : call xgScalapack_free(scalapack)
669 : : end if
670 : 726392 : tsec(2) = abi_wtime() - tsec(2)
671 [ + + ]: 726392 : if ( prtvol == 4 ) write(std_out,*) tsec(2)
672 : :
673 : 726392 : call timab(tim_RR_diago,2,tsec)
674 : : ABI_NVTX_END_RANGE()
675 : :
676 [ - + ]: 726392 : if ( eigenSolver == EIGENVX .or. EIGPACK(eigenSolver)) then
677 : 0 : call xg_free(subA)
678 : : end if
679 : 726392 : call xg_free(subB)
680 : :
681 : 726392 : call timab(tim_RR_gemm_2,1,tsec)
682 : : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_2)
683 : :
684 : : !FIXME Avoid those transfers
685 [ + - ]: 726392 : if ( info == 0 ) then
686 : 726392 : call xg_init(subB,space(X),spacedim,blockdim,comm=comm(X),me_g0=me_g0(X),gpu_option=gpu_option)
687 : :
688 : : !/* Easy basic solution */
689 : : !/* Compute first part of X here */
690 : : ! Use subB as buffer
691 : : !lobpcg%XWP (:,X+1:X+blockdim) = matmul(lobpcg%XWP (:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
692 : 726392 : call xgBlock_setBlock(vec%self,Cwp,blockdim,blockdim)
693 : 726392 : call xgBlock_gemm('n','n',1.0d0,X,Cwp,0.d0,subB%self)
694 : 726392 : call xgBlock_copy(subB%self,X)
695 : :
696 : : !lobpcg%AXWP(:,X+1:X+blockdim) = matmul(lobpcg%AXWP(:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
697 : 726392 : call xgBlock_gemm('n','n',1.0d0,AX,Cwp,0.d0,subB%self)
698 : 726392 : call xgBlock_copy(subB%self,AX)
699 : :
700 : : !lobpcg%BXWP(:,X+1:X+blockdim) = matmul(lobpcg%BXWP(:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
701 : 726392 : call xgBlock_gemm('n','n',1.0d0,BX,Cwp,0.d0,subB%self)
702 : 726392 : call xgBlock_copy(subB%self,BX)
703 : :
704 [ + + ]: 726392 : if ( var /= VAR_X ) then
705 : : ! Cost to pay to avoid temporary array in xgemm
706 [ - + ]: 557470 : if(gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(vec%self) !FIXME Avoid that transfer
707 : 557470 : call xgBlock_cshift(vec%self,blockdim,1) ! Bottom 2*blockdim lines are now at the top
708 [ - + ]: 557470 : if(gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_to_gpu(vec%self) !FIXME Avoid that transfer
709 : 557470 : call xgBlock_setBlock(vec%self,Cwp,subdim-blockdim,blockdim)
710 : :
711 : : !lobpcg%XWP (:,P+1:P+blockdim) = matmul(lobpcg%XWP (:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
712 : 557470 : call xgBlock_gemm('n','n',1.0d0,WP,Cwp,0.d0,subB%self)
713 : 557470 : call xgBlock_copy(subB%self,P)
714 : :
715 : : !lobpcg%AXWP(:,P+1:P+blockdim) = matmul(lobpcg%AXWP(:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
716 : 557470 : call xgBlock_gemm('n','n',1.0d0,AWP,Cwp,0.d0,subB%self)
717 : 557470 : call xgBlock_copy(subB%self,AP)
718 : :
719 : : !lobpcg%BXWP(:,P+1:P+blockdim) = matmul(lobpcg%BXWP(:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
720 : 557470 : call xgBlock_gemm('n','n',1.0d0,BWP,Cwp,0.d0,subB%self)
721 : 557470 : call xgBlock_copy(subB%self,BP)
722 : :
723 : : !/* Maybe faster solution
724 : : ! * Sum previous contribution plus P direction
725 : : ! */
726 : 557470 : call xgBlock_add(X,P)
727 : 557470 : call xgBlock_add(AX,AP)
728 : 557470 : call xgBlock_add(BX,BP)
729 : : end if
730 : : end if
731 : :
732 : 726392 : call timab(tim_RR_gemm_2,2,tsec)
733 : :
734 : : ! Doing free on an already free object does not do anything
735 : 726392 : call xg_free(vec)
736 : 726392 : call xg_free(subA)
737 : 726392 : call xg_free(subB)
738 : :
739 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
740 : : if (gpu_option==ABI_GPU_KOKKOS) then
741 : : call gpu_device_synchronize()
742 : : end if
743 : : #endif
744 : :
745 : : ABI_NVTX_END_RANGE()
746 : 726392 : call timab(timer , 2, tsec)
747 : :
748 [ + + ]: 2179176 : end subroutine xg_RayleighRitz
749 : : !!***
750 : :
751 : : !!****f* m_xg_ortho_RR/xg_RayleighRitz_cprj
752 : : !!
753 : : !! NAME
754 : : !! xg_RayleighRitz_cprj
755 : 94541 : subroutine xg_RayleighRitz_cprj(xg_nonlop,X,cprjX,AX,eigenvalues,info,prtvol,timer,gpu_option,&
756 : : tolerance,XW,W,cprjXW,cprjW,AW,P,cprjP,AP,WP,cprjWP,AWP,XWP,cprjXWP,blockdim_cprj,solve_ax_bx,add_Anl,drop_null)
757 : :
758 : : integer , intent(in ) :: timer,gpu_option
759 : : integer , intent(in ) :: prtvol
760 : : type(xgBlock_t) , intent(inout) :: eigenvalues
761 : : type(xgBlock_t) , intent(inout) :: X
762 : : type(xgBlock_t) , intent(inout) :: cprjX
763 : : type(xgBlock_t) , intent(inout) :: AX
764 : : type(xg_nonlop_t), intent(in ) :: xg_nonlop
765 : : integer , intent( out) :: info
766 : : double precision, optional, intent(in) :: tolerance
767 : : ! LOBPCG only :
768 : : type(xgBlock_t), intent(inout),optional :: XW
769 : : type(xgBlock_t), intent(inout),optional :: W
770 : : type(xgBlock_t), intent(inout),optional :: cprjXW
771 : : type(xgBlock_t), intent(inout),optional :: cprjW
772 : : type(xgBlock_t), intent(inout),optional :: AW
773 : : type(xgBlock_t), intent(inout),optional :: P
774 : : type(xgBlock_t), intent(inout),optional :: cprjP
775 : : type(xgBlock_t), intent(inout),optional :: AP
776 : : type(xgBlock_t), intent(inout),optional :: WP
777 : : type(xgBlock_t), intent(inout),optional :: cprjWP
778 : : type(xgBlock_t), intent(inout),optional :: AWP
779 : : type(xgBlock_t), intent(inout),optional :: XWP
780 : : type(xgBlock_t), intent(inout),optional :: cprjXWP
781 : : integer,intent(in),optional :: blockdim_cprj
782 : : ! End LOBPCG only
783 : : logical, intent(in),optional :: solve_ax_bx
784 : : logical, intent(in),optional :: add_Anl
785 : : logical, intent(in),optional :: drop_null
786 : : integer :: var
787 : : integer :: spacedim
788 : : integer :: blockdim,blockdim_cprj_
789 : : integer :: space_buf
790 : : integer :: subdim
791 : : integer :: spacecom
792 : : integer :: eigenSolver
793 : : integer :: nrows_B, ncols_B
794 : : integer :: nspinor,cprjdim,ncols_cprj
795 : : !integer :: neigen
796 : : double precision :: abstol
797 : : double precision :: cond
798 : : #ifdef HAVE_LINALG_SCALAPACK
799 : : logical :: use_slk
800 : : #endif
801 : : type(xg_t) :: vec
802 : : type(xg_t) :: subA
803 : : type(xg_t) :: subB
804 : : type(xg_t) :: Xwork
805 : : type(xg_t) :: cprjXwork
806 : : type(xg_t) :: cprj_work
807 : : type(xgBlock_t) :: cprj_workX,cprj_workX_spinor,cprjX_spinor,cprjW_spinor,cprjWP_spinor
808 : : type(xgBlock_t) :: subsub
809 : : type(xgBlock_t) :: Cwp
810 : : type(xgScalapack_t) :: scalapack
811 : : double precision :: tsec(2)
812 : : logical :: solve_ax_bx_,add_Anl_,drop_null_
813 : :
814 : 94541 : call timab(timer , 1, tsec)
815 : :
816 [ - + ]: 94541 : if (gpu_option/=ABI_GPU_DISABLED) then
817 : 0 : ABI_ERROR('Not implemented for GPU')
818 : : end if
819 : 94541 : call xgBlock_check_gpu_option(X,cprjX)
820 : 94541 : call xgBlock_check_gpu_option(X,AX)
821 : :
822 : 94541 : solve_ax_bx_ = .false.
823 [ + + ]: 94541 : if (present(solve_ax_bx)) then
824 : 15368 : solve_ax_bx_ = solve_ax_bx
825 : : end if
826 : 94541 : add_Anl_ = .false.
827 [ + + ]: 94541 : if (present(add_Anl)) then
828 : 79173 : add_Anl_ = add_Anl
829 : : end if
830 : 94541 : drop_null_ = .false.
831 [ + + ]: 94541 : if (present(drop_null)) then
832 : 8968 : drop_null_ = drop_null
833 : : end if
834 : :
835 : 94541 : call timab(timer , 1, tsec)
836 : :
837 : 94541 : var = VAR_X
838 : 94541 : call xgBlock_check(X,AX)
839 [ + + ]: 94541 : if (.not.solve_ax_bx_) then
840 : 79173 : eigenSolver = EIGENEVD
841 : : else
842 : 15368 : eigenSolver = EIGENVD
843 : : end if
844 : 94541 : spacedim = rows(X)
845 : 94541 : blockdim = cols(X)
846 : 94541 : ncols_cprj = cols(cprjX)
847 : 94541 : blockdim_cprj_ = ncols_cprj
848 [ + + ]: 94541 : if (present(blockdim_cprj)) then
849 : 58560 : blockdim_cprj_ = blockdim_cprj
850 : : end if
851 : 94541 : space_buf = space(X)
852 [ + + ]: 94541 : if (space(X)==SPACE_CR) then
853 : 44525 : space_buf = SPACE_R
854 : : end if
855 : 94541 : spacecom = comm(X)
856 : 94541 : subdim = blockdim
857 : 94541 : nspinor = xg_nonlop%nspinor
858 : 94541 : cprjdim = xg_nonlop%cprjdim
859 : :
860 [ + + ]: 94541 : if (present(XW)) then
861 [ - + ]: 58560 : if (solve_ax_bx_) then
862 : 0 : ABI_ERROR('solve_ax_bx is not implemented in that case')
863 : : end if
864 : 58560 : var = VAR_XW
865 : 58560 : eigenSolver = EIGENVD
866 : : !TODO Do checks on other optional arguments
867 : 58560 : call xgBlock_check(X,AW)
868 : 58560 : call xgBlock_check_gpu_option(X,AW)
869 : 58560 : call xgBlock_check(X,P)
870 : 58560 : call xgBlock_check_gpu_option(X,P)
871 : 58560 : call xgBlock_check(X,AP)
872 : 58560 : call xgBlock_check_gpu_option(X,AP)
873 : 58560 : call xgBlock_check(X,XW,fact_col=2)
874 : 58560 : call xgBlock_check_gpu_option(X,XW)
875 : 58560 : call xgBlock_check(X,WP,fact_col=2)
876 : 58560 : call xgBlock_check_gpu_option(X,WP)
877 : 58560 : call xgBlock_check(X,AWP,fact_col=2)
878 : 58560 : call xgBlock_check_gpu_option(X,AWP)
879 : 58560 : subdim = 2*blockdim
880 [ + + ]: 58560 : if (present(XWP)) then
881 : 49592 : var = VAR_XWP
882 : 49592 : call xgBlock_check(X,XWP,fact_col=3)
883 : 49592 : subdim = 3*blockdim
884 : : end if
885 : : end if
886 : :
887 : : #ifdef HAVE_LINALG_SCALAPACK
888 : : call xgScalapack_init(scalapack,spacecom,subdim,prtvol-2,(gpu_option/=ABI_GPU_DISABLED),use_slk)
889 : : if (use_slk) then
890 : : eigenSolver = EIGENSLK
891 : : end if
892 : : #endif
893 : :
894 : : ! drop_null: use xg_hegv_drop (LAPACK on every rank, no ScaLAPACK), which drops the null directions of the Gram
895 : : ! matrix subB. subB only exists for the generalized problem, i.e. not for X alone without solve_ax_bx.
896 [ + + ]: 94541 : if (drop_null_) then
897 [ - + ]: 8 : if (var == VAR_X .and. .not.solve_ax_bx_) then
898 : 0 : ABI_BUG("drop_null needs the generalized eigenproblem")
899 : : end if
900 [ - + ]: 8 : if ( eigenSolver == EIGENSLK ) then
901 : 0 : call xgScalapack_free(scalapack)
902 : : end if
903 : 8 : eigenSolver = EIGENDROP
904 : : end if
905 : :
906 [ + + ]: 94541 : abstol = 0d0 ; if ( present(tolerance) ) abstol = tolerance
907 : :
908 : 94541 : call xg_init(subA,space_buf,subdim,subdim,spacecom)
909 [ + + ]: 94541 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
910 : 73928 : call xg_init(subB,space_buf,subdim,subdim,spacecom)
911 : : end if
912 : :
913 [ - + ]: 94541 : if ( eigenSolver == EIGENVX .or. eigenSolver == EIGENPVX ) then
914 : 0 : call xg_init(vec,space_buf,subdim,blockdim)
915 [ - + ]: 94541 : else if ( EIGPACK(eigenSolver) ) then
916 : 0 : call xg_init(vec,space_buf,subdim,subdim)
917 : : else
918 : 94541 : call xg_setBlock(subA,vec%self,subdim,blockdim)
919 : : endif
920 : :
921 : : ! Compute subA and subB by part
922 : : !--- begin
923 : : ! | E | XAW | XAP | | I | XBW | XBP |
924 : : ! | * | WAW | WAP | | * | I | WBP |
925 : : ! | * | * | PAP | | * | * | I |
926 : :
927 : 94541 : call timab(tim_RR_gemm_1,1,tsec)
928 : : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_1)
929 : :
930 : : ! Compute XAX
931 : 94541 : call xg_setBlock(subA,subsub,blockdim,blockdim)
932 : :
933 : 94541 : call xgBlock_gemm('t','n',1.0d0,X,AX,0.d0,subsub,comm=spacecom)
934 : :
935 : : ! Add the nonlocal part (H)
936 : 94541 : call xg_init(cprj_work,space(cprjX),rows(cprjX),cols(cprjX),spacecom)
937 [ + + ]: 94541 : if (add_Anl_) then
938 : 79173 : call xg_nonlop_getXHY(xg_nonlop,cprjX,cprjX,cprj_work%self,subsub,blocksize=blockdim_cprj_)
939 : : end if
940 : :
941 [ + + ]: 94541 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
942 : : ! Compute XBX
943 : 73928 : call xg_setBlock(subB,subsub,blockdim,blockdim)
944 : 73928 : call xgBlock_gemm('t','n',1.0d0,X,X,0.d0,subsub,comm=spacecom)
945 : : ! Add the nonlocal part (S)
946 [ + + ]: 73928 : if (xg_nonlop%paw) then
947 : 61873 : call xg_nonlop_getXSY(xg_nonlop,cprjX,cprjX,cprj_work%self,subsub,blocksize=blockdim_cprj_)
948 : : end if
949 : : end if
950 : :
951 [ + + ]: 94541 : if ( var == VAR_XW .or. var == VAR_XWP ) then
952 : :
953 : : ! Compute XAW and WAW
954 : 58560 : call xg_setBlock(subA,subsub,2*blockdim,blockdim,fcol=blockdim+1)
955 : 58560 : call xgBlock_gemm('t','n',1.0d0,XW,AW,0.d0,subsub,comm=spacecom)
956 : :
957 : : ! Add the nonlocal part (H)
958 [ + - ]: 58560 : if (add_Anl_) then
959 : 58560 : call xg_nonlop_getXHY(xg_nonlop,cprjXW,cprjW,cprj_work%self,subsub,blocksize=blockdim_cprj_)
960 : : end if
961 : :
962 : : ! Compute XBW and WBW
963 : 58560 : call xg_setBlock(subB,subsub,2*blockdim,blockdim,fcol=blockdim+1)
964 : 58560 : call xgBlock_gemm('t','n',1.0d0,XW,W,0.d0,subsub,comm=spacecom)
965 : : ! Add the nonlocal part (S)
966 [ + + ]: 58560 : if (xg_nonlop%paw) then
967 : 48953 : call xg_nonlop_getXSY(xg_nonlop,cprjXW,cprjW,cprj_work%self,subsub,blocksize=blockdim_cprj_)
968 : : end if
969 : :
970 : : end if
971 : :
972 [ + + ]: 58560 : if ( var == VAR_XWP ) then
973 : : ! Compute XAP, WAP and PAP
974 : 49592 : call xg_setBlock(subA,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
975 : 49592 : call xgBlock_gemm('t','n',1.0d0,XWP,AP,0.d0,subsub,comm=spacecom)
976 : :
977 : : ! Add the nonlocal part (H)
978 [ + - ]: 49592 : if (add_Anl_) then
979 : 49592 : call xg_nonlop_getXHY(xg_nonlop,cprjXWP,cprjP,cprj_work%self,subsub,blocksize=blockdim_cprj_)
980 : : end if
981 : :
982 : 49592 : call xg_setBlock(subB,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
983 : 49592 : call xgBlock_gemm('t','n',1.0d0,XWP,P,0.d0,subsub,comm=spacecom)
984 : : ! Add the nonlocal part (S)
985 [ + + ]: 49592 : if (xg_nonlop%paw) then
986 : 41333 : call xg_nonlop_getXSY(xg_nonlop,cprjXWP,cprjP,cprj_work%self,subsub,blocksize=blockdim_cprj_)
987 : : end if
988 : : end if
989 : :
990 : 94541 : call xg_free(cprj_work)
991 : :
992 : 94541 : call timab(tim_RR_gemm_1,2,tsec)
993 : : ABI_NVTX_END_RANGE()
994 : :
995 [ - + ]: 94541 : if ( EIGPACK(eigenSolver) ) then
996 : 0 : call xgBlock_pack(subA%self,subA%self,'u')
997 : : end if
998 : :
999 : 94541 : call timab(tim_RR_diago,1,tsec)
1000 : 94541 : tsec(2) = abi_wtime()
1001 [ + + ]: 94541 : if ( .not.solve_ax_bx_ .and. var == VAR_X ) then
1002 : : ! Solve Hermitian eigen problem
1003 : 20613 : select case (eigenSolver)
1004 : : case (EIGENEVD)
1005 [ - + ]: 20613 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heevd"
1006 : 20613 : call xgBlock_heevd('v','u',subA%self,eigenvalues,info)
1007 : : case (EIGENEV)
1008 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heev"
1009 : 0 : call xgBlock_heev('v','u',subA%self,eigenvalues,info)
1010 : : case (EIGENPEVD)
1011 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpevd"
1012 : 0 : call xgBlock_hpevd('v','u',subA%self,eigenvalues,vec%self,info)
1013 : : case (EIGENPEV)
1014 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpev"
1015 : 0 : call xgBlock_hpev('v','u',subA%self,eigenvalues,vec%self,info)
1016 : : case (EIGENSLK)
1017 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using pheev"
1018 : 0 : call xgScalapack_heev(scalapack,subA%self,eigenvalues)
1019 : 0 : info = 0 ! No error code returned for the moment
1020 : : case default
1021 [ + - - - : 41226 : ABI_ERROR("Error for Eigen Solver HEEV")
- - ]
1022 : : end select
1023 : : else
1024 : : ! Solve Hermitian general eigen problem only for first blockdim eigenvalues
1025 : 0 : select case (eigenSolver)
1026 : : case (EIGENVX)
1027 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvx"
1028 : : call xgBlock_hegvx(1,'v','i','u',subA%self,subB%self,0.d0,0.d0,1,blockdim,abstol,&
1029 : 0 : eigenvalues,vec%self,info)
1030 : : case (EIGENVD)
1031 : : !ITEST
1032 [ - + ]: 73920 : if (prtvol == 15015015) then
1033 : 0 : write(902,*) 'eigenSolver', eigenSolver
1034 : 0 : write(902,*) 'Using hegvd'
1035 : 0 : nrows_B = rows(subB%self)
1036 : 0 : ncols_B = cols(subB%self)
1037 : 0 : write(902,*) 'computing cond for Hermitian pd matrix', nrows_B, ncols_B
1038 : : !# Validation
1039 : : !nrows_B = 2
1040 : : !#
1041 : 0 : call xgBlock_hermitian_pd_cond(subB%self, nrows_B, cond)
1042 : 0 : write(902,*) 'cond(B)=', cond
1043 : : !# Validation: should be 2.076578056
1044 : 0 : flush(902)
1045 : : end if
1046 : : !ITEST
1047 [ - + ]: 73920 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvd"
1048 : 73920 : call xgBlock_hegvd(1,'v','u',subA%self,subB%self,eigenvalues,info)
1049 : : case (EIGENV)
1050 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegv"
1051 : 0 : call xgBlock_hegv(1,'v','u',subA%self,subB%self,eigenvalues,info)
1052 : : case (EIGENPVX)
1053 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgvx"
1054 : : call xgBlock_hpgvx(1,'v','i','u',subA%self,subB%self,0.d0,0.d0,1,blockdim,abstol,&
1055 : 0 : eigenvalues,vec%self,info)
1056 : : case (EIGENPVD)
1057 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgvd"
1058 : 0 : call xgBlock_hpgvd(1,'v','u',subA%self,subB%self,eigenvalues,vec%self,info)
1059 : : case (EIGENPV)
1060 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgv"
1061 : 0 : call xgBlock_hpgv(1,'v','u',subA%self,subB%self,eigenvalues,vec%self,info)
1062 : : case (EIGENSLK)
1063 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using phegv"
1064 : 0 : call xgScalapack_hegv(scalapack,subA%self,subB%self,eigenvalues)
1065 : 0 : info = 0 ! No error code returned for the moment
1066 : : case (EIGENDROP)
1067 [ - + ]: 8 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegv with dropping"
1068 : 8 : call xg_hegv_drop(subA%self,subB%self,eigenvalues,blockdim,gpu_option,info)
1069 : : case default
1070 [ - + - - : 73936 : ABI_ERROR("Error for Eigen Solver HEGV")
- - - +
- ]
1071 : : end select
1072 : : end if
1073 [ - + ]: 94541 : if ( eigenSolver == EIGENSLK ) then
1074 : 0 : call xgScalapack_free(scalapack)
1075 : : end if
1076 : 94541 : tsec(2) = abi_wtime() - tsec(2)
1077 [ - + ]: 94541 : if ( prtvol == 4 ) write(std_out,*) tsec(2)
1078 : 94541 : call timab(tim_RR_diago,2,tsec)
1079 : :
1080 [ + - - + ]: 94541 : if ( eigenSolver == EIGENVX .or. EIGPACK(eigenSolver)) then
1081 : 0 : call xg_free(subA)
1082 : : end if
1083 : 94541 : call xg_free(subB)
1084 : :
1085 : 94541 : call timab(tim_RR_gemm_2,1,tsec)
1086 : : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_2)
1087 : :
1088 [ + - ]: 94541 : if ( info == 0 ) then
1089 : 94541 : call xg_init(Xwork,space(X),spacedim,blockdim,me_g0=me_g0(X))
1090 : 94541 : call xg_init(cprjXwork,space(cprjX),cprjdim,ncols_cprj,comm(cprjX))
1091 : 94541 : call xg_setBlock(cprjXwork,cprj_workX,cprjdim,ncols_cprj)
1092 : 94541 : call xgBlock_reshape_spinor(cprj_workX,cprj_workX_spinor,nspinor,COLS2ROWS)
1093 : :
1094 : 94541 : call xgBlock_setBlock(vec%self,Cwp,blockdim,blockdim)
1095 : :
1096 : : !/* Easy basic solution */
1097 : : !/* Compute first part of X here */
1098 : : !XWP (:,X+1:X+blockdim) = matmul(XWP (:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
1099 : 94541 : call xgBlock_gemm('n','n',1.0d0,X,Cwp,0.d0,Xwork%self)
1100 : 94541 : call xgBlock_copy(Xwork%self,X)
1101 : :
1102 : : !AXWP(:,X+1:X+blockdim) = matmul(AXWP(:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
1103 : 94541 : call xgBlock_gemm('n','n',1.0d0,AX,Cwp,0.d0,Xwork%self)
1104 : 94541 : call xgBlock_copy(Xwork%self,AX)
1105 : :
1106 : : !cprjXWP(:,X+1:X+blockdim) = matmul(cprjXWP(:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
1107 : 94541 : call xgBlock_reshape_spinor(cprjX,cprjX_spinor,nspinor,COLS2ROWS)
1108 [ + + ]: 94541 : if (ncols_cprj==nspinor*cols(X)) then
1109 : 32307 : call xgBlock_gemm('n','n',1.0d0,cprjX_spinor,Cwp,0.d0,cprj_workX_spinor)
1110 : : else
1111 : 62234 : call xgBlock_zero(cprj_workX)
1112 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjX_spinor,Cwp,cprj_workX_spinor,&
1113 : 62234 : & xg_nonlop%me_band,blocksize=blockdim_cprj_/nspinor,comm=xg_nonlop%comm_band)
1114 : : end if
1115 : 94541 : call xgBlock_copy(cprj_workX,cprjX)
1116 : :
1117 [ + + ]: 94541 : if ( var /= VAR_X ) then
1118 : : ! Cost to pay to avoid temporary array in xgemm
1119 : 58560 : call xgBlock_cshift(vec%self,blockdim,1) ! Bottom 2*blockdim lines are now at the top
1120 : 58560 : call xgBlock_setBlock(vec%self,Cwp,subdim-blockdim,blockdim)
1121 : :
1122 : : !XWP (:,P+1:P+blockdim) = matmul(XWP (:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
1123 : 58560 : call xgBlock_gemm('n','n',1.0d0,WP,Cwp,0.d0,Xwork%self)
1124 : 58560 : call xgBlock_copy(Xwork%self,P)
1125 : :
1126 : : !AXWP(:,P+1:P+blockdim) = matmul(AXWP(:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
1127 : 58560 : call xgBlock_gemm('n','n',1.0d0,AWP,Cwp,0.d0,Xwork%self)
1128 : 58560 : call xgBlock_copy(Xwork%self,AP)
1129 : :
1130 : : !cprjXWP(:,P+1:P+blockdim) = matmul(cprjXWP(:,W+1:X+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
1131 : 58560 : call xgBlock_reshape_spinor(cprjWP,cprjWP_spinor,nspinor,COLS2ROWS)
1132 [ + + ]: 58560 : if (ncols_cprj==cols(WP)) then
1133 : 1863 : call xgBlock_gemm('n','n',1.0d0,cprjWP_spinor,Cwp,0.d0,cprj_workX_spinor)
1134 : : else
1135 : 56697 : call xgBlock_zero(cprj_workX)
1136 [ + + ]: 56697 : if ( var==VAR_XW ) then
1137 : 8685 : call xgBlock_reshape_spinor(cprjW,cprjW_spinor,nspinor,COLS2ROWS)
1138 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjW_spinor,Cwp,cprj_workX_spinor,&
1139 : 8685 : & xg_nonlop%me_band,blocksize=blockdim_cprj_/nspinor,comm=xg_nonlop%comm_band)
1140 : : else if ( var==VAR_XWP ) then
1141 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjWP_spinor,Cwp,cprj_workX_spinor,&
1142 : 48012 : & xg_nonlop%me_band,blocksize=blockdim_cprj_/nspinor,comm=xg_nonlop%comm_band)
1143 : : else
1144 : : ABI_ERROR('not implemented')
1145 : : end if
1146 : : end if
1147 : 58560 : call xgBlock_copy(cprj_workX,cprjP)
1148 : :
1149 : : !/* Maybe faster solution
1150 : : ! * Sum previous contribution plus P direction
1151 : : ! */
1152 : 58560 : call xgBlock_add(X,P)
1153 : 58560 : call xgBlock_add(AX,AP)
1154 : 58560 : call xgBlock_add(cprjX,cprjP)
1155 : : end if
1156 : :
1157 : : end if
1158 : :
1159 : 94541 : call timab(tim_RR_gemm_2,2,tsec)
1160 : : ABI_NVTX_END_RANGE()
1161 : :
1162 : : ! Doing free on an already free object does not do anything
1163 : 94541 : call xg_free(vec)
1164 : 94541 : call xg_free(subA)
1165 : 94541 : call xg_free(Xwork)
1166 : 94541 : call xg_free(cprjXwork)
1167 : :
1168 : 94541 : call timab(timer , 2, tsec)
1169 : :
1170 [ + + ]: 283623 : end subroutine xg_RayleighRitz_cprj
1171 : : !!***
1172 : :
1173 : : !!****f* m_xg_ortho_RR/xg_hegv_drop
1174 : : !!
1175 : : !! NAME
1176 : : !! xg_hegv_drop
1177 : : !!
1178 : : !! FUNCTION
1179 : : !! Solve the generalized eigenproblem A.c = lambda.B.c of a Rayleigh-Ritz step whose basis was
1180 : : !! B-orthonormalized with a regularized (floored) Gram matrix, by dropping the null directions of B.
1181 : : !! After the regularization, B is the identity except in the floored directions, where it is the
1182 : : !! ratio of the true Gram eigenvalue to the floor: a small but accurate value for a real direction,
1183 : : !! and rounding noise for a direction that was null. The noise directions have no meaningful A, and
1184 : : !! would give arbitrary Ritz values: they are the directions with an eigenvalue of B below
1185 : : !! tol_drop*mu_max, and are dropped. The standard eigenproblem is solved in the other directions.
1186 : : !! Only the upper triangles of A and B are read, as in hegvd. The nsel lowest eigenpairs are
1187 : : !! returned as hegvd does: eigenvalues in the first nsel rows of eigenvalues, B-orthonormal
1188 : : !! eigenvectors in the first nsel columns of A. B is destroyed.
1189 : : !! info/=0 if B is not positive semi-definite, if fewer than nsel directions are kept,
1190 : : !! or if the reduced eigenproblem fails.
1191 : : !!
1192 : : !! SOURCE
1193 : 28 : subroutine xg_hegv_drop(subA,subB,eigenvalues,nsel,gpu_option,info)
1194 : :
1195 : : type(xgBlock_t), intent(inout) :: subA
1196 : : type(xgBlock_t), intent(inout) :: subB
1197 : : type(xgBlock_t), intent(inout) :: eigenvalues
1198 : : integer , intent(in ) :: nsel
1199 : : integer , intent(in ) :: gpu_option
1200 : : integer , intent(out ) :: info
1201 : : ! The rounding noise left in a null direction is ~1e-16 relative to mu_max, a real floored
1202 : : ! direction is above ~1e-4 (Gram eigenvalue above its own rounding level, floor at tol12)
1203 : : double precision, parameter :: tol_drop = tol12
1204 : : integer :: nn,nkeep,space_buf,ii,jj,kk
1205 : : logical :: valid
1206 : : double precision :: mu_max,mu_drop
1207 : 14 : double precision, pointer :: mu(:,:),lam(:,:),ev(:,:),aR(:,:),vR(:,:)
1208 : 14 : complex(dp), pointer :: aC(:,:)
1209 : : type(xg_t) :: evals_B,evals_A,work,subAr
1210 : : type(xgBlock_t) :: tmat,ymat,cmat,cview
1211 : : character(len=500) :: msg
1212 : :
1213 : 14 : nn = cols(subA)
1214 : 14 : nkeep = 0
1215 : 14 : space_buf = space(subA)
1216 [ + - + - : 14 : if (rows(subA)/=nn .or. rows(subB)/=nn .or. cols(subB)/=nn .or. rows(eigenvalues)<nsel) then
- + ]
1217 : 0 : ABI_BUG("xg_hegv_drop: inconsistent sizes")
1218 : : end if
1219 : :
1220 : : ! Fill the lower triangle of A, which is used in full in the products below
1221 [ - + ]: 14 : if (gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(subA)
1222 : 0 : select case (space_buf)
1223 : : case (SPACE_R)
1224 : 0 : call xgBlock_reverseMap(subA,aR)
1225 [ # # ]: 0 : do jj=1,nn
1226 [ # # ]: 0 : do ii=1,jj-1
1227 : 0 : aR(jj,ii) = aR(ii,jj)
1228 : : end do
1229 : : end do
1230 : : case (SPACE_C)
1231 : 14 : call xgBlock_reverseMap(subA,aC)
1232 [ + + ]: 238 : do jj=1,nn
1233 [ + + ]: 1904 : do ii=1,jj-1
1234 : 1904 : aC(jj,ii) = conjg(aC(ii,jj))
1235 : : end do
1236 : 238 : aC(jj,jj) = cmplx(dble(aC(jj,jj)),zero,kind=dp)
1237 : : end do
1238 : : case default
1239 [ - + - ]: 14 : ABI_BUG("xg_hegv_drop: space should be SPACE_R or SPACE_C")
1240 : : end select
1241 [ - + ]: 14 : if (gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_to_gpu(subA)
1242 : :
1243 : : ! Diagonalize B: subB <- eigenvectors V, evals_B <- eigenvalues mu (ascending)
1244 : 14 : call xg_init(evals_B,SPACE_R,nn,1,gpu_option=gpu_option)
1245 : 14 : call xgBlock_heevd('v','u',subB,evals_B%self,info)
1246 : :
1247 [ + - ]: 14 : if (info==0) then
1248 [ - + ]: 14 : if (gpu_option==ABI_GPU_OPENMP) then
1249 : 0 : call xgBlock_copy_from_gpu(evals_B%self)
1250 : 0 : call xgBlock_copy_from_gpu(subB)
1251 : : end if
1252 : 14 : call xgBlock_reverseMap(evals_B%self,mu)
1253 : 14 : mu_max = mu(nn,1)
1254 : : ! A valid B is positive semi-definite: a NaN, or an eigenvalue clearly below zero, means that B
1255 : : ! is not an overlap matrix. Then nothing is dropped and info/=0, as hegvd would fail.
1256 [ + + + - : 476 : valid = all(mu(1:nn,1)==mu(1:nn,1)) .and. all(abs(mu(1:nn,1))<=huge(one))
+ + + - ]
1257 [ + - + - : 14 : if (valid) valid = (mu_max>zero .and. mu(1,1)>=-tol_drop*mu_max)
+ - ]
1258 : : if (valid) then
1259 : : ! The eigenvalues are ascending: the kept directions are the last nkeep ones.
1260 : : nkeep = 0
1261 [ + + ]: 238 : do kk=1,nn
1262 [ + + ]: 238 : if (mu(kk,1)>tol_drop*mu_max) nkeep = nkeep+1
1263 : : end do
1264 : 14 : mu_drop = zero
1265 [ + - ]: 14 : if (nkeep<nn) mu_drop = mu(nn-nkeep,1)
1266 : 14 : write(msg,'(a,i0,a,i0,4(a,es10.3))') 'Rayleigh-Ritz after regularization: n=',nn,&
1267 : 14 : & ', dropped=',nn-nkeep,', mu_min=',mu(1,1),', mu_max=',mu_max,&
1268 : 28 : & ', largest dropped=',mu_drop,', smallest kept=',mu(nn-nkeep+1,1)
1269 : 14 : ABI_COMMENT(msg)
1270 [ - + ]: 14 : if (nkeep<nsel) then
1271 : 0 : write(msg,'(2(a,i0))') 'Too few directions left after dropping: kept=',nkeep,', needed=',nsel
1272 : 0 : ABI_COMMENT(msg)
1273 : 0 : info = 1
1274 : : end if
1275 : : else
1276 : 0 : write(msg,'(2(a,es10.3),a)') 'Rayleigh-Ritz: B is not positive semi-definite (mu_min=',mu(1,1),&
1277 : 0 : & ', mu_max=',mu_max,'). Eigenproblem not solved'
1278 : 0 : ABI_COMMENT(msg)
1279 : 0 : info = 1
1280 : : end if
1281 : : end if
1282 : :
1283 [ + - ]: 14 : if (info==0) then
1284 : : ! T = V(:,kept).diag(mu(kept))^(-1/2), in place in subB
1285 : 14 : call xgBlock_reverseMap(subB,vR)
1286 [ + + ]: 176 : do kk=nn-nkeep+1,nn
1287 [ + + ]: 5360 : vR(:,kk) = vR(:,kk)/sqrt(mu(kk,1))
1288 : : end do
1289 [ - + ]: 14 : if (gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_to_gpu(subB)
1290 : 14 : call xgBlock_setBlock(subB,tmat,nn,nkeep,fcol=nn-nkeep+1)
1291 : : ! Reduced standard eigenproblem: T^H.A.T
1292 : 14 : call xg_init(work,space_buf,nn,nkeep,gpu_option=gpu_option)
1293 : 14 : call xgBlock_gemm('n','n',1.0d0,subA,tmat,0.d0,work%self)
1294 : 14 : call xg_init(subAr,space_buf,nkeep,nkeep,gpu_option=gpu_option)
1295 : 14 : call xgBlock_gemm('t','n',1.0d0,tmat,work%self,0.d0,subAr%self)
1296 : 14 : call xg_init(evals_A,SPACE_R,nkeep,1,gpu_option=gpu_option)
1297 : 14 : call xgBlock_heevd('v','u',subAr%self,evals_A%self,info)
1298 : : end if
1299 : :
1300 [ + - ]: 14 : if (info==0) then
1301 [ - + ]: 14 : if (gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(evals_A%self)
1302 : 14 : call xgBlock_reverseMap(evals_A%self,lam)
1303 [ + + + - : 252 : if (any(lam(1:nsel,1)/=lam(1:nsel,1)) .or. any(abs(lam(1:nsel,1))>huge(one))) then
+ + + - -
+ ]
1304 : 0 : ABI_COMMENT("Non-finite eigenvalues after dropping the null directions")
1305 : 0 : info = 1
1306 : : end if
1307 : : end if
1308 : :
1309 [ + - ]: 14 : if (info==0) then
1310 : : ! Eigenvectors in the original basis, C = T.Y(:,1:nsel), in the first nsel columns of A
1311 : 14 : call xgBlock_setBlock(subAr%self,ymat,nkeep,nsel)
1312 : 14 : call xgBlock_setBlock(work%self,cmat,nn,nsel)
1313 : 14 : call xgBlock_gemm('n','n',1.0d0,tmat,ymat,0.d0,cmat)
1314 : 14 : call xgBlock_setBlock(subA,cview,nn,nsel)
1315 : 14 : call xgBlock_copy(cmat,cview)
1316 [ - + ]: 14 : if (gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(eigenvalues)
1317 : 14 : call xgBlock_reverseMap(eigenvalues,ev)
1318 [ + - + + : 238 : ev(1:nsel,1) = lam(1:nsel,1)
+ + ]
1319 [ - + ]: 14 : if (gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_to_gpu(eigenvalues)
1320 : : end if
1321 : :
1322 : 14 : call xg_free(evals_B)
1323 : 14 : call xg_free(evals_A)
1324 : 14 : call xg_free(work)
1325 : 14 : call xg_free(subAr)
1326 : :
1327 [ + - ]: 14 : end subroutine xg_hegv_drop
1328 : : !!***
1329 : :
1330 : : end module m_xg_ortho_RR
1331 : : !!***
|
