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 : : logical, parameter :: EIGPACK(11) = &
65 : : (/ .false.,.false.,.false., &
66 : : .true. ,.true. ,.true. ,&
67 : : .false.,.false.,&
68 : : .true. ,.true., .false. /)
69 : :
70 : : integer, parameter :: tim_RR_diago = 1795
71 : : integer, parameter :: tim_RR_gemm_1 = 1796
72 : : integer, parameter :: tim_RR_gemm_2 = 1797
73 : :
74 : : public :: xg_Borthonormalize
75 : : public :: xg_Borthonormalize_cprj
76 : : public :: xg_RayleighRitz
77 : : public :: xg_RayleighRitz_cprj
78 : :
79 : : contains
80 : : !!***
81 : :
82 : : !!****f* m_xg_ortho_RR/xg_Borthonormalize
83 : : !!
84 : : !! NAME
85 : : !! xg_Borthonormalize
86 : 1453838 : subroutine xg_Borthonormalize(X,BX,info,timer,gpu_option,AX)
87 : :
88 : : integer , intent(in ) :: timer
89 : : integer , intent(in ) :: gpu_option
90 : : type(xgBlock_t), intent(inout) :: X
91 : : type(xgBlock_t), intent(inout) :: BX
92 : : type(xgBlock_t), intent(inout),optional :: AX
93 : : integer , intent( out) :: info
94 : : integer :: space_buf
95 : : type(xg_t) :: buffer
96 : : double precision :: tsec(2)
97 : :
98 : : ABI_NVTX_START_RANGE(NVTX_B_ORTHO)
99 : :
100 : 726919 : call timab(timer,1,tsec)
101 : :
102 : 726919 : call xgBlock_check(X,BX)
103 [ + - ]: 726919 : if (present(AX)) then
104 : 726919 : call xgBlock_check(X,AX)
105 : : end if
106 : :
107 [ + + ]: 726919 : if (space(X)/=SPACE_CR) then
108 : 670132 : space_buf = SPACE(X)
109 : : else
110 : 56787 : space_buf = SPACE_R
111 : : end if
112 : 726919 : call xg_init(buffer,space_buf,cols(X),cols(X),comm(X),gpu_option=gpu_option)
113 : :
114 : : ! If space(X)==SPACE_CR : set imaginary part of G=0 component to zero to improve numerical stability
115 : 726919 : call xgBlock_zero_im_g0(X)
116 : 726919 : call xgBlock_zero_im_g0(BX)
117 [ + - ]: 726919 : if (present(AX)) then
118 : 726919 : call xgBlock_zero_im_g0(AX)
119 : : end if
120 : :
121 : : ! Compute X^TBX
122 : 726919 : call xgBlock_gemm('t','n',1.d0,X,BX,0.d0,buffer%self,comm=comm(X))
123 : :
124 : : ! Compute Cholesky decomposition (Upper part)
125 : 726919 : call xgBlock_potrf(buffer%self,'u',info)
126 : :
127 [ + + ]: 726919 : if ( info /= 0 ) then
128 : 8 : ABI_COMMENT("Cholesky decomposition did not work. Orthonormalization not done")
129 : 8 : call xg_free(buffer)
130 : : ABI_NVTX_END_RANGE()
131 : 8 : return
132 : : end if
133 : :
134 : : ! Solve YU=X
135 : 726911 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,X)
136 : :
137 : : ! Solve BYU=BX
138 : 726911 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,BX)
139 : :
140 [ + - ]: 726911 : if (present(AX)) then
141 : : ! Solve AYU=AX
142 : 726911 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,AX)
143 : : end if
144 : :
145 : 726911 : call xg_free(buffer)
146 : :
147 : : ABI_NVTX_END_RANGE()
148 : 726911 : call timab(timer,2,tsec)
149 : :
150 : 726919 : end subroutine xg_Borthonormalize
151 : : !!***
152 : :
153 : : !!****f* m_xg_ortho_RR/xg_Borthonormalize_cprj
154 : : !!
155 : : !! NAME
156 : : !! xg_Borthonormalize_cprj
157 : 158200 : subroutine xg_Borthonormalize_cprj(xg_nonlop,X,cprjX,info,timer,gpu_option,blockdim_cprj,AX)
158 : :
159 : : integer , intent(in ) :: timer,gpu_option
160 : : integer , intent( out) :: info
161 : : type(xg_nonlop_t), intent(in ) :: xg_nonlop
162 : : type(xgBlock_t) , intent(inout) :: X
163 : : type(xgBlock_t) , intent(inout) :: cprjX
164 : : integer , intent(in ),optional :: blockdim_cprj
165 : : type(xgBlock_t) , intent(inout),optional :: AX
166 : :
167 : : type(xg_t) :: buffer,cprj_work
168 : : type(xgBlock_t) :: cprjX_spinor,cprj_work_spinor
169 : : integer :: blockdim_cprj_,space_buf
170 : : integer :: spacecom,ncols_cprj,nn,nspinor
171 : : double precision :: tsec(2)
172 : :
173 : 79100 : call timab(timer,1,tsec)
174 : :
175 [ - + ]: 79100 : if (gpu_option/=ABI_GPU_DISABLED) then
176 : 0 : ABI_ERROR('Not implemented for GPU')
177 : : end if
178 : 79100 : call xgBlock_check_gpu_option(X,cprjX)
179 [ + - ]: 79100 : if (present(AX)) then
180 : 79100 : call xgBlock_check_gpu_option(X,AX)
181 : : end if
182 : :
183 : 79100 : nn = cols(X)
184 : 79100 : ncols_cprj = cols(cprjX)
185 : 79100 : blockdim_cprj_ = ncols_cprj
186 [ + + ]: 79100 : if (present(blockdim_cprj)) then
187 : 58504 : blockdim_cprj_ = blockdim_cprj
188 : : end if
189 : 79100 : nspinor = xg_nonlop%nspinor
190 : :
191 : 79100 : spacecom = comm(X)
192 : :
193 [ + + ]: 79100 : if (space(X)/=SPACE_CR) then
194 : 40831 : space_buf = SPACE(X)
195 : : else
196 : 38269 : space_buf = SPACE_R
197 : : end if
198 : 79100 : call xg_init(buffer,space_buf,nn,nn,spacecom)
199 : :
200 : : ! Compute X^TX
201 : 79100 : call xgBlock_gemm('t','n',1.d0,X,X,0.d0,buffer%self,comm=spacecom)
202 : 79100 : call xg_init(cprj_work,space(cprjX),rows(cprjX),ncols_cprj,spacecom)
203 [ + + ]: 79100 : if (xg_nonlop%paw) then
204 : 65413 : call xg_nonlop_getXSY(xg_nonlop,cprjX,cprjX,cprj_work%self,buffer%self,blocksize=blockdim_cprj_)
205 : : end if
206 : :
207 : : ! Compute Cholesky decomposition (Upper part)
208 : 79100 : call xgBlock_potrf(buffer%self,'u',info)
209 : :
210 [ - + ]: 79100 : if ( info /= 0 ) then
211 : 0 : ABI_COMMENT("Cholesky decomposition did not work. Orthonormalization not done")
212 : 0 : call xg_free(buffer)
213 : 0 : return
214 : : end if
215 : :
216 : : ! Solve YU=X
217 : 79100 : call xgBlock_trsm('r','u','n','n',1.d0,buffer%self,X)
218 : :
219 [ + - ]: 79100 : if (present(AX)) then
220 : : ! Solve AYU=AX
221 : 79100 : call xgBlock_trsm('r','u','n','n',1.d0,buffer%self,AX)
222 : : end if
223 : :
224 : : ! Solve (cprjY)U=(cprjX)
225 : 79100 : call xgBlock_reshape_spinor(cprjX,cprjX_spinor,nspinor,COLS2ROWS)
226 [ + + ]: 79100 : if (ncols_cprj==nspinor*nn) then
227 : 28018 : call xgBlock_trsm('r','u','n','n',1.d0,buffer%self,cprjX_spinor)
228 : : else
229 : 51082 : call xgBlock_invert_tri('u','n',buffer%self)
230 : 51082 : call xgBlock_zerotri(buffer%self,'u')
231 : 51082 : call xgBlock_zero(cprj_work%self)
232 : 51082 : call xgBlock_reshape_spinor(cprj_work%self,cprj_work_spinor,nspinor,COLS2ROWS)
233 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjX_spinor,buffer%self,cprj_work_spinor,&
234 : 51082 : & xg_nonlop%me_band,blocksize=blockdim_cprj_/nspinor,comm=xg_nonlop%comm_band)
235 : 51082 : call xgBlock_copy(cprj_work%self,cprjX)
236 : : end if
237 : :
238 : 79100 : call xg_free(cprj_work)
239 : :
240 : 79100 : call xg_free(buffer)
241 : :
242 : 79100 : call timab(timer,2,tsec)
243 : :
244 : 79100 : end subroutine xg_Borthonormalize_cprj
245 : : !!***
246 : :
247 : : !!****f* m_xg_ortho_RR/xg_RayleighRitz
248 : : !!
249 : : !! NAME
250 : : !! xg_RayleighRitz
251 : 738167 : subroutine xg_RayleighRitz(X,AX,BX,eigenvalues,info,prtvol,timer,gpu_option,&
252 : : & tolerance,XW,AW,BW,P,AP,BP,WP,AWP,BWP,XWP,solve_ax_bx)
253 : :
254 : : integer , intent(in ) :: timer
255 : : integer , intent(in ) :: gpu_option
256 : : integer , intent(in ) :: prtvol
257 : : type(xgBlock_t), intent(inout) :: eigenvalues
258 : : type(xgBlock_t), intent(inout) :: X
259 : : type(xgBlock_t), intent(inout) :: AX
260 : : type(xgBlock_t), intent(inout) :: BX
261 : : integer , intent( out) :: info
262 : : double precision, optional, intent(in) :: tolerance
263 : : ! LOBPCG only :
264 : : type(xgBlock_t), intent(inout),optional :: XW
265 : : type(xgBlock_t), intent(inout),optional :: AW
266 : : type(xgBlock_t), intent(inout),optional :: BW
267 : : type(xgBlock_t), intent(inout),optional :: P
268 : : type(xgBlock_t), intent(inout),optional :: AP
269 : : type(xgBlock_t), intent(inout),optional :: BP
270 : : type(xgBlock_t), intent(inout),optional :: WP
271 : : type(xgBlock_t), intent(inout),optional :: AWP
272 : : type(xgBlock_t), intent(inout),optional :: BWP
273 : : type(xgBlock_t), intent(inout),optional :: XWP
274 : : ! End LOBPCG only
275 : : logical, intent(in),optional :: solve_ax_bx
276 : : integer :: var
277 : : integer :: spacedim
278 : : integer :: blockdim
279 : : integer :: space_buf
280 : : integer :: subdim
281 : : integer :: spacecom
282 : : integer :: eigenSolver
283 : : double precision :: abstol
284 : : #ifdef HAVE_LINALG_SCALAPACK
285 : : logical :: use_slk
286 : : #endif
287 : : type(xg_t) :: vec
288 : : type(xg_t) :: subA
289 : : type(xg_t) :: subB
290 : : type(xgBlock_t) :: subsub
291 : : type(xgBlock_t) :: Cwp
292 : : type(xgScalapack_t) :: scalapack
293 : : double precision :: tsec(2)
294 : : logical :: solve_ax_bx_
295 : :
296 : 738167 : call timab(timer , 1, tsec)
297 : :
298 : 738167 : solve_ax_bx_ = .false.
299 [ + + ]: 738167 : if (present(solve_ax_bx)) then
300 : 11256 : solve_ax_bx_ = solve_ax_bx
301 : : end if
302 : :
303 : 738167 : var = VAR_X
304 : 738167 : call xgBlock_check(X,AX)
305 : 738167 : call xgBlock_check(X,BX)
306 [ + + ]: 738167 : if (.not.solve_ax_bx_) then
307 : : eigenSolver = EIGENEVD
308 : : else
309 : 11256 : eigenSolver = EIGENVD
310 : : end if
311 : 738167 : spacedim = rows(X)
312 : 738167 : blockdim = cols(X)
313 : 738167 : space_buf = space(X)
314 [ + + ]: 738167 : if (space(X) == space_CR) then
315 : 60454 : space_buf = space_R
316 : : end if
317 : 738167 : spacecom = comm(X)
318 : 738167 : subdim = blockdim
319 : :
320 [ + + ]: 738167 : if (present(XW)) then
321 [ - + ]: 566854 : if (solve_ax_bx_) then
322 : 0 : ABI_ERROR('solve_ax_bx is not implemented in that case')
323 : : end if
324 : 566854 : var = VAR_XW
325 : 566854 : eigenSolver = EIGENVD
326 : 566854 : call xgBlock_check(X,AW)
327 : 566854 : call xgBlock_check(X,BW)
328 : 566854 : call xgBlock_check(X,P)
329 : 566854 : call xgBlock_check(X,AP)
330 : 566854 : call xgBlock_check(X,BP)
331 : 566854 : call xgBlock_check(X,XW,fact_col=2)
332 : 566854 : call xgBlock_check(X,WP,fact_col=2)
333 : 566854 : call xgBlock_check(X,AWP,fact_col=2)
334 : 566854 : call xgBlock_check(X,BWP,fact_col=2)
335 : 566854 : subdim = 2*blockdim
336 [ + + ]: 566854 : if (present(XWP)) then
337 : 430484 : var = VAR_XWP
338 : 430484 : call xgBlock_check(X,XWP,fact_col=3)
339 : 430484 : subdim = 3*blockdim
340 : : end if
341 : : end if
342 : :
343 : : #ifdef HAVE_LINALG_SCALAPACK
344 : : call xgScalapack_init(scalapack,spacecom,subdim,prtvol-2,(gpu_option/=ABI_GPU_DISABLED),use_slk)
345 : : if ( use_slk) then
346 : : eigenSolver = EIGENSLK
347 : : end if
348 : : #endif
349 : :
350 : : ! Select diago algorithm
351 : :
352 [ + + ]: 738167 : abstol = 0d0 ; if ( present(tolerance) ) abstol = tolerance
353 : :
354 : 738167 : call xg_init(subA,space_buf,subdim,subdim,spacecom,gpu_option=gpu_option)
355 [ + + ]: 738167 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
356 : 578110 : call xg_init(subB,space_buf,subdim,subdim,spacecom,gpu_option=gpu_option)
357 : : end if
358 : :
359 : : if ( eigenSolver == EIGENVX .or. eigenSolver == EIGENPVX ) then
360 : : call xg_init(vec,space_buf,subdim,blockdim,gpu_option=gpu_option)
361 [ - + ]: 738167 : else if ( EIGPACK(eigenSolver) ) then
362 : 0 : call xg_init(vec,space_buf,subdim,subdim,gpu_option=gpu_option)
363 : : else
364 : 738167 : call xg_setBlock(subA,vec%self,subdim,blockdim)
365 : : endif
366 : :
367 : : ! Compute subA and subB by part
368 : : !--- begin
369 : : ! | E | XAW | XAP | | I | XBW | XBP |
370 : : ! | * | WAW | WAP | | * | I | WBP |
371 : : ! | * | * | PAP | | * | * | I |
372 : :
373 : 738167 : call timab(tim_RR_gemm_1,1,tsec)
374 : : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_1)
375 : :
376 : : ! If space(X)==SPACE_CR : set imaginary part of G=0 component to zero to improve numerical stability
377 [ + + ]: 738167 : if (var == VAR_X) call xgBlock_zero_im_g0(X)
378 [ + + ]: 738167 : if (var == VAR_XW) call xgBlock_zero_im_g0(XW)
379 [ + + ]: 738167 : if (var == VAR_XWP) call xgBlock_zero_im_g0(XWP)
380 : 738167 : call xgBlock_zero_im_g0(AX)
381 : 738167 : call xgBlock_zero_im_g0(BX)
382 : :
383 : 738167 : call xg_setBlock(subA,subsub,blockdim,blockdim)
384 : 738167 : call xgBlock_gemm('t','n',1.0d0,X,AX,0.d0,subsub,comm=spacecom)
385 : :
386 [ + + ]: 738167 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
387 : 578110 : call xg_setBlock(subB,subsub,blockdim,blockdim)
388 : 578110 : call xgBlock_gemm('t','n',1.0d0,X,BX,0.d0,subsub,comm=spacecom)
389 : : endif
390 : :
391 [ + + ]: 738167 : if ( var == VAR_XW .or. var == VAR_XWP ) then
392 : :
393 : : ! If space(X)==SPACE_CR : set imaginary part of G=0 component to zero to improve numerical stability
394 : 566854 : call xgBlock_zero_im_g0(AW)
395 : 566854 : call xgBlock_zero_im_g0(BW)
396 : :
397 : : ! subA
398 : 566854 : call xg_setBlock(subA,subsub,2*blockdim,blockdim,fcol=blockdim+1)
399 : 566854 : call xgBlock_gemm('t','n',1.0d0,XW,AW,0.d0,subsub,comm=spacecom)
400 : :
401 : : ! subB
402 : 566854 : call xg_setBlock(subB,subsub,2*blockdim,blockdim,fcol=blockdim+1)
403 : 566854 : call xgBlock_gemm('t','n',1.0d0,XW,BW,0.d0,subsub,comm=spacecom)
404 : :
405 : : end if
406 : :
407 [ + + ]: 566854 : if ( var == VAR_XWP ) then
408 : :
409 : : ! If space(X)==SPACE_CR : set imaginary part of G=0 component to zero to improve numerical stability
410 : 430484 : call xgBlock_zero_im_g0(AP)
411 : 430484 : call xgBlock_zero_im_g0(BP)
412 : :
413 : : ! subA
414 : 430484 : call xg_setBlock(subA,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
415 : 430484 : call xgBlock_gemm('t','n',1.0d0,XWP,AP,0.d0,subsub,comm=spacecom)
416 : :
417 : : ! subB
418 : 430484 : call xg_setBlock(subB,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
419 : 430484 : call xgBlock_gemm('t','n',1.0d0,XWP,BP,0.d0,subsub,comm=spacecom)
420 : :
421 : : end if
422 : :
423 : 738167 : call timab(tim_RR_gemm_1,2,tsec)
424 : : ABI_NVTX_END_RANGE()
425 : :
426 [ - + ]: 738167 : if ( EIGPACK(eigenSolver) ) then
427 : 0 : call xgBlock_pack(subA%self,subA%self,'u')
428 [ # # ]: 0 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
429 : 0 : call xgBlock_pack(subB%self,subB%self,'u')
430 : : end if
431 : : end if
432 : :
433 : 738167 : call timab(tim_RR_diago,1,tsec)
434 : 738167 : tsec(2) = abi_wtime()
435 [ + + ]: 738167 : if ( .not.solve_ax_bx_ .and. var == VAR_X ) then
436 : : ABI_NVTX_START_RANGE(NVTX_RR_HEEV)
437 : : ! Solve Hermitian eigen problem
438 : 160057 : select case (eigenSolver)
439 : : case (EIGENEVD)
440 [ + + ]: 160057 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heevd"
441 : 160057 : call xgBlock_heevd('v','u',subA%self,eigenvalues,info) ! work with GPU
442 : : case (EIGENEV)
443 : : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heev"
444 : : call xgBlock_heev('v','u',subA%self,eigenvalues,info)
445 : : case (EIGENPEVD)
446 : : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpevd"
447 : : call xgBlock_hpevd('v','u',subA%self,eigenvalues,vec%self,info)
448 : : case (EIGENPEV)
449 : : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpev"
450 : : call xgBlock_hpev('v','u',subA%self,eigenvalues,vec%self,info)
451 : : case (EIGENSLK)
452 : : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using pheev"
453 : : call xgScalapack_heev(scalapack,subA%self,eigenvalues,gpu_option=gpu_option) ! work with GPU
454 : 0 : info = 0 ! No error code returned for the moment
455 : : case default
456 [ + - ]: 320114 : ABI_ERROR("Error for Eigen Solver HEEV")
457 : : end select
458 : : else
459 : : ABI_NVTX_START_RANGE(NVTX_RR_HEGV)
460 : : ! Solve Hermitian general eigen problem only for first blockdim eigenvalues
461 : : select case (eigenSolver)
462 : : case (EIGENVX)
463 : : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvx"
464 : : call xgBlock_hegvx(1,'v','i','u',subA%self,subB%self,0.d0,0.d0,1,blockdim,abstol,&
465 : 578110 : eigenvalues,vec%self,info)
466 : : case (EIGENVD)
467 [ + + ]: 578110 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvd"
468 : 578110 : call xgBlock_hegvd(1,'v','u',subA%self,subB%self,eigenvalues,info) ! work with GPU
469 : : case (EIGENV)
470 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegv"
471 : 0 : call xgBlock_hegv(1,'v','u',subA%self,subB%self,eigenvalues,info)
472 : : case (EIGENPVX)
473 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgvx"
474 : : call xgBlock_hpgvx(1,'v','i','u',subA%self,subB%self,0.d0,0.d0,1,blockdim,abstol,&
475 : 0 : eigenvalues,vec%self,info)
476 : : case (EIGENPVD)
477 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgvd"
478 : 0 : call xgBlock_hpgvd(1,'v','u',subA%self,subB%self,eigenvalues,vec%self,info)
479 : : case (EIGENPV)
480 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgv"
481 : 0 : call xgBlock_hpgv(1,'v','u',subA%self,subB%self,eigenvalues,vec%self,info)
482 : : case (EIGENSLK)
483 : : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using phegv"
484 : : call xgScalapack_hegv(scalapack,subA%self,subB%self,eigenvalues,gpu_option=gpu_option) ! work with GPU
485 : 0 : info = 0 ! No error code returned for the moment
486 : : case default
487 [ + - - - : 578110 : ABI_ERROR("Error for Eigen Solver HEGV")
- - ]
488 : : end select
489 : : end if
490 : : if ( eigenSolver == EIGENSLK ) then
491 : : call xgScalapack_free(scalapack)
492 : : end if
493 : 738167 : tsec(2) = abi_wtime() - tsec(2)
494 [ + + ]: 738167 : if ( prtvol == 4 ) write(std_out,*) tsec(2)
495 : :
496 : 738167 : call timab(tim_RR_diago,2,tsec)
497 : : ABI_NVTX_END_RANGE()
498 : :
499 [ - + ]: 738167 : if ( eigenSolver == EIGENVX .or. EIGPACK(eigenSolver)) then
500 : 0 : call xg_free(subA)
501 : : end if
502 : 738167 : call xg_free(subB)
503 : :
504 : 738167 : call timab(tim_RR_gemm_2,1,tsec)
505 : : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_2)
506 : :
507 : : !FIXME Avoid those transfers
508 [ + - ]: 738167 : if ( info == 0 ) then
509 : 738167 : call xg_init(subB,space(X),spacedim,blockdim,comm=comm(X),me_g0=me_g0(X),gpu_option=gpu_option)
510 : :
511 : : !/* Easy basic solution */
512 : : !/* Compute first part of X here */
513 : : ! Use subB as buffer
514 : : !lobpcg%XWP (:,X+1:X+blockdim) = matmul(lobpcg%XWP (:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
515 : 738167 : call xgBlock_setBlock(vec%self,Cwp,blockdim,blockdim)
516 : 738167 : call xgBlock_gemm('n','n',1.0d0,X,Cwp,0.d0,subB%self)
517 : 738167 : call xgBlock_copy(subB%self,X)
518 : :
519 : : !lobpcg%AXWP(:,X+1:X+blockdim) = matmul(lobpcg%AXWP(:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
520 : 738167 : call xgBlock_gemm('n','n',1.0d0,AX,Cwp,0.d0,subB%self)
521 : 738167 : call xgBlock_copy(subB%self,AX)
522 : :
523 : : !lobpcg%BXWP(:,X+1:X+blockdim) = matmul(lobpcg%BXWP(:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
524 : 738167 : call xgBlock_gemm('n','n',1.0d0,BX,Cwp,0.d0,subB%self)
525 : 738167 : call xgBlock_copy(subB%self,BX)
526 : :
527 [ + + ]: 738167 : if ( var /= VAR_X ) then
528 : : ! Cost to pay to avoid temporary array in xgemm
529 [ - + ]: 566854 : if(gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(vec%self) !FIXME Avoid that transfer
530 : 566854 : call xgBlock_cshift(vec%self,blockdim,1) ! Bottom 2*blockdim lines are now at the top
531 [ - + ]: 566854 : if(gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_to_gpu(vec%self) !FIXME Avoid that transfer
532 : 566854 : call xgBlock_setBlock(vec%self,Cwp,subdim-blockdim,blockdim)
533 : :
534 : : !lobpcg%XWP (:,P+1:P+blockdim) = matmul(lobpcg%XWP (:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
535 : 566854 : call xgBlock_gemm('n','n',1.0d0,WP,Cwp,0.d0,subB%self)
536 : 566854 : call xgBlock_copy(subB%self,P)
537 : :
538 : : !lobpcg%AXWP(:,P+1:P+blockdim) = matmul(lobpcg%AXWP(:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
539 : 566854 : call xgBlock_gemm('n','n',1.0d0,AWP,Cwp,0.d0,subB%self)
540 : 566854 : call xgBlock_copy(subB%self,AP)
541 : :
542 : : !lobpcg%BXWP(:,P+1:P+blockdim) = matmul(lobpcg%BXWP(:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
543 : 566854 : call xgBlock_gemm('n','n',1.0d0,BWP,Cwp,0.d0,subB%self)
544 : 566854 : call xgBlock_copy(subB%self,BP)
545 : :
546 : : !/* Maybe faster solution
547 : : ! * Sum previous contribution plus P direction
548 : : ! */
549 : 566854 : call xgBlock_add(X,P)
550 : 566854 : call xgBlock_add(AX,AP)
551 : 566854 : call xgBlock_add(BX,BP)
552 : : end if
553 : : end if
554 : :
555 : 738167 : call timab(tim_RR_gemm_2,2,tsec)
556 : :
557 : : ! Doing free on an already free object does not do anything
558 : 738167 : call xg_free(vec)
559 : 738167 : call xg_free(subA)
560 : 738167 : call xg_free(subB)
561 : :
562 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
563 : : if (gpu_option==ABI_GPU_KOKKOS) then
564 : : call gpu_device_synchronize()
565 : : end if
566 : : #endif
567 : :
568 : : ABI_NVTX_END_RANGE()
569 : 738167 : call timab(timer , 2, tsec)
570 : :
571 [ + + ]: 2214501 : end subroutine xg_RayleighRitz
572 : : !!***
573 : :
574 : : !!****f* m_xg_ortho_RR/xg_RayleighRitz_cprj
575 : : !!
576 : : !! NAME
577 : : !! xg_RayleighRitz_cprj
578 : 94468 : subroutine xg_RayleighRitz_cprj(xg_nonlop,X,cprjX,AX,eigenvalues,info,prtvol,timer,gpu_option,&
579 : : tolerance,XW,W,cprjXW,cprjW,AW,P,cprjP,AP,WP,cprjWP,AWP,XWP,cprjXWP,blockdim_cprj,solve_ax_bx,add_Anl)
580 : :
581 : : integer , intent(in ) :: timer,gpu_option
582 : : integer , intent(in ) :: prtvol
583 : : type(xgBlock_t) , intent(inout) :: eigenvalues
584 : : type(xgBlock_t) , intent(inout) :: X
585 : : type(xgBlock_t) , intent(inout) :: cprjX
586 : : type(xgBlock_t) , intent(inout) :: AX
587 : : type(xg_nonlop_t), intent(in ) :: xg_nonlop
588 : : integer , intent( out) :: info
589 : : double precision, optional, intent(in) :: tolerance
590 : : ! LOBPCG only :
591 : : type(xgBlock_t), intent(inout),optional :: XW
592 : : type(xgBlock_t), intent(inout),optional :: W
593 : : type(xgBlock_t), intent(inout),optional :: cprjXW
594 : : type(xgBlock_t), intent(inout),optional :: cprjW
595 : : type(xgBlock_t), intent(inout),optional :: AW
596 : : type(xgBlock_t), intent(inout),optional :: P
597 : : type(xgBlock_t), intent(inout),optional :: cprjP
598 : : type(xgBlock_t), intent(inout),optional :: AP
599 : : type(xgBlock_t), intent(inout),optional :: WP
600 : : type(xgBlock_t), intent(inout),optional :: cprjWP
601 : : type(xgBlock_t), intent(inout),optional :: AWP
602 : : type(xgBlock_t), intent(inout),optional :: XWP
603 : : type(xgBlock_t), intent(inout),optional :: cprjXWP
604 : : integer,intent(in),optional :: blockdim_cprj
605 : : ! End LOBPCG only
606 : : logical, intent(in),optional :: solve_ax_bx
607 : : logical, intent(in),optional :: add_Anl
608 : : integer :: var
609 : : integer :: spacedim
610 : : integer :: blockdim,blockdim_cprj_
611 : : integer :: space_buf
612 : : integer :: subdim
613 : : integer :: spacecom
614 : : integer :: eigenSolver
615 : : integer :: nrows_B, ncols_B
616 : : integer :: nspinor,cprjdim,ncols_cprj
617 : : !integer :: neigen
618 : : double precision :: abstol
619 : : double precision :: cond
620 : : #ifdef HAVE_LINALG_SCALAPACK
621 : : logical :: use_slk
622 : : #endif
623 : : type(xg_t) :: vec
624 : : type(xg_t) :: subA
625 : : type(xg_t) :: subB
626 : : type(xg_t) :: Xwork
627 : : type(xg_t) :: cprjXwork
628 : : type(xg_t) :: cprj_work
629 : : type(xgBlock_t) :: cprj_workX,cprj_workX_spinor,cprjX_spinor,cprjW_spinor,cprjWP_spinor
630 : : type(xgBlock_t) :: subsub
631 : : type(xgBlock_t) :: Cwp
632 : : type(xgScalapack_t) :: scalapack
633 : : double precision :: tsec(2)
634 : : logical :: solve_ax_bx_,add_Anl_
635 : :
636 : 94468 : call timab(timer , 1, tsec)
637 : :
638 [ - + ]: 94468 : if (gpu_option/=ABI_GPU_DISABLED) then
639 : 0 : ABI_ERROR('Not implemented for GPU')
640 : : end if
641 : 94468 : call xgBlock_check_gpu_option(X,cprjX)
642 : 94468 : call xgBlock_check_gpu_option(X,AX)
643 : :
644 : 94468 : solve_ax_bx_ = .false.
645 [ + + ]: 94468 : if (present(solve_ax_bx)) then
646 : 15368 : solve_ax_bx_ = solve_ax_bx
647 : : end if
648 : 94468 : add_Anl_ = .false.
649 [ + + ]: 94468 : if (present(add_Anl)) then
650 : 79100 : add_Anl_ = add_Anl
651 : : end if
652 : :
653 : 94468 : call timab(timer , 1, tsec)
654 : :
655 : 94468 : var = VAR_X
656 : 94468 : call xgBlock_check(X,AX)
657 [ + + ]: 94468 : if (.not.solve_ax_bx_) then
658 : 79100 : eigenSolver = EIGENEVD
659 : : else
660 : 15368 : eigenSolver = EIGENVD
661 : : end if
662 : 94468 : spacedim = rows(X)
663 : 94468 : blockdim = cols(X)
664 : 94468 : ncols_cprj = cols(cprjX)
665 : 94468 : blockdim_cprj_ = ncols_cprj
666 [ + + ]: 94468 : if (present(blockdim_cprj)) then
667 : 58504 : blockdim_cprj_ = blockdim_cprj
668 : : end if
669 : 94468 : space_buf = space(X)
670 [ + + ]: 94468 : if (space(X)==SPACE_CR) then
671 : 44525 : space_buf = SPACE_R
672 : : end if
673 : 94468 : spacecom = comm(X)
674 : 94468 : subdim = blockdim
675 : 94468 : nspinor = xg_nonlop%nspinor
676 : 94468 : cprjdim = xg_nonlop%cprjdim
677 : :
678 [ + + ]: 94468 : if (present(XW)) then
679 [ - + ]: 58504 : if (solve_ax_bx_) then
680 : 0 : ABI_ERROR('solve_ax_bx is not implemented in that case')
681 : : end if
682 : 58504 : var = VAR_XW
683 : 58504 : eigenSolver = EIGENVD
684 : : !TODO Do checks on other optional arguments
685 : 58504 : call xgBlock_check(X,AW)
686 : 58504 : call xgBlock_check_gpu_option(X,AW)
687 : 58504 : call xgBlock_check(X,P)
688 : 58504 : call xgBlock_check_gpu_option(X,P)
689 : 58504 : call xgBlock_check(X,AP)
690 : 58504 : call xgBlock_check_gpu_option(X,AP)
691 : 58504 : call xgBlock_check(X,XW,fact_col=2)
692 : 58504 : call xgBlock_check_gpu_option(X,XW)
693 : 58504 : call xgBlock_check(X,WP,fact_col=2)
694 : 58504 : call xgBlock_check_gpu_option(X,WP)
695 : 58504 : call xgBlock_check(X,AWP,fact_col=2)
696 : 58504 : call xgBlock_check_gpu_option(X,AWP)
697 : 58504 : subdim = 2*blockdim
698 [ + + ]: 58504 : if (present(XWP)) then
699 : 49566 : var = VAR_XWP
700 : 49566 : call xgBlock_check(X,XWP,fact_col=3)
701 : 49566 : subdim = 3*blockdim
702 : : end if
703 : : end if
704 : :
705 : : #ifdef HAVE_LINALG_SCALAPACK
706 : : call xgScalapack_init(scalapack,spacecom,subdim,prtvol-2,(gpu_option/=ABI_GPU_DISABLED),use_slk)
707 : : if (use_slk) then
708 : : eigenSolver = EIGENSLK
709 : : end if
710 : : #endif
711 : :
712 [ + + ]: 94468 : abstol = 0d0 ; if ( present(tolerance) ) abstol = tolerance
713 : :
714 : 94468 : call xg_init(subA,space_buf,subdim,subdim,spacecom)
715 [ + + ]: 94468 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
716 : 73872 : call xg_init(subB,space_buf,subdim,subdim,spacecom)
717 : : end if
718 : :
719 [ - + ]: 94468 : if ( eigenSolver == EIGENVX .or. eigenSolver == EIGENPVX ) then
720 : 0 : call xg_init(vec,space_buf,subdim,blockdim)
721 [ - + ]: 94468 : else if ( EIGPACK(eigenSolver) ) then
722 : 0 : call xg_init(vec,space_buf,subdim,subdim)
723 : : else
724 : 94468 : call xg_setBlock(subA,vec%self,subdim,blockdim)
725 : : endif
726 : :
727 : : ! Compute subA and subB by part
728 : : !--- begin
729 : : ! | E | XAW | XAP | | I | XBW | XBP |
730 : : ! | * | WAW | WAP | | * | I | WBP |
731 : : ! | * | * | PAP | | * | * | I |
732 : :
733 : 94468 : call timab(tim_RR_gemm_1,1,tsec)
734 : : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_1)
735 : :
736 : : ! Compute XAX
737 : 94468 : call xg_setBlock(subA,subsub,blockdim,blockdim)
738 : :
739 : 94468 : call xgBlock_gemm('t','n',1.0d0,X,AX,0.d0,subsub,comm=spacecom)
740 : :
741 : : ! Add the nonlocal part (H)
742 : 94468 : call xg_init(cprj_work,space(cprjX),rows(cprjX),cols(cprjX),spacecom)
743 [ + + ]: 94468 : if (add_Anl_) then
744 : 79100 : call xg_nonlop_getXHY(xg_nonlop,cprjX,cprjX,cprj_work%self,subsub,blocksize=blockdim_cprj_)
745 : : end if
746 : :
747 [ + + ]: 94468 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
748 : : ! Compute XBX
749 : 73872 : call xg_setBlock(subB,subsub,blockdim,blockdim)
750 : 73872 : call xgBlock_gemm('t','n',1.0d0,X,X,0.d0,subsub,comm=spacecom)
751 : : ! Add the nonlocal part (S)
752 [ + + ]: 73872 : if (xg_nonlop%paw) then
753 : 61817 : call xg_nonlop_getXSY(xg_nonlop,cprjX,cprjX,cprj_work%self,subsub,blocksize=blockdim_cprj_)
754 : : end if
755 : : end if
756 : :
757 [ + + ]: 94468 : if ( var == VAR_XW .or. var == VAR_XWP ) then
758 : :
759 : : ! Compute XAW and WAW
760 : 58504 : call xg_setBlock(subA,subsub,2*blockdim,blockdim,fcol=blockdim+1)
761 : 58504 : call xgBlock_gemm('t','n',1.0d0,XW,AW,0.d0,subsub,comm=spacecom)
762 : :
763 : : ! Add the nonlocal part (H)
764 [ + - ]: 58504 : if (add_Anl_) then
765 : 58504 : call xg_nonlop_getXHY(xg_nonlop,cprjXW,cprjW,cprj_work%self,subsub,blocksize=blockdim_cprj_)
766 : : end if
767 : :
768 : : ! Compute XBW and WBW
769 : 58504 : call xg_setBlock(subB,subsub,2*blockdim,blockdim,fcol=blockdim+1)
770 : 58504 : call xgBlock_gemm('t','n',1.0d0,XW,W,0.d0,subsub,comm=spacecom)
771 : : ! Add the nonlocal part (S)
772 [ + + ]: 58504 : if (xg_nonlop%paw) then
773 : 48897 : call xg_nonlop_getXSY(xg_nonlop,cprjXW,cprjW,cprj_work%self,subsub,blocksize=blockdim_cprj_)
774 : : end if
775 : :
776 : : end if
777 : :
778 [ + + ]: 58504 : if ( var == VAR_XWP ) then
779 : : ! Compute XAP, WAP and PAP
780 : 49566 : call xg_setBlock(subA,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
781 : 49566 : call xgBlock_gemm('t','n',1.0d0,XWP,AP,0.d0,subsub,comm=spacecom)
782 : :
783 : : ! Add the nonlocal part (H)
784 [ + - ]: 49566 : if (add_Anl_) then
785 : 49566 : call xg_nonlop_getXHY(xg_nonlop,cprjXWP,cprjP,cprj_work%self,subsub,blocksize=blockdim_cprj_)
786 : : end if
787 : :
788 : 49566 : call xg_setBlock(subB,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
789 : 49566 : call xgBlock_gemm('t','n',1.0d0,XWP,P,0.d0,subsub,comm=spacecom)
790 : : ! Add the nonlocal part (S)
791 [ + + ]: 49566 : if (xg_nonlop%paw) then
792 : 41307 : call xg_nonlop_getXSY(xg_nonlop,cprjXWP,cprjP,cprj_work%self,subsub,blocksize=blockdim_cprj_)
793 : : end if
794 : : end if
795 : :
796 : 94468 : call xg_free(cprj_work)
797 : :
798 : 94468 : call timab(tim_RR_gemm_1,2,tsec)
799 : : ABI_NVTX_END_RANGE()
800 : :
801 [ - + ]: 94468 : if ( EIGPACK(eigenSolver) ) then
802 : 0 : call xgBlock_pack(subA%self,subA%self,'u')
803 : : end if
804 : :
805 : 94468 : call timab(tim_RR_diago,1,tsec)
806 : 94468 : tsec(2) = abi_wtime()
807 [ + + ]: 94468 : if ( .not.solve_ax_bx_ .and. var == VAR_X ) then
808 : : ! Solve Hermitian eigen problem
809 : 20596 : select case (eigenSolver)
810 : : case (EIGENEVD)
811 [ - + ]: 20596 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heevd"
812 : 20596 : call xgBlock_heevd('v','u',subA%self,eigenvalues,info)
813 : : case (EIGENEV)
814 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heev"
815 : 0 : call xgBlock_heev('v','u',subA%self,eigenvalues,info)
816 : : case (EIGENPEVD)
817 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpevd"
818 : 0 : call xgBlock_hpevd('v','u',subA%self,eigenvalues,vec%self,info)
819 : : case (EIGENPEV)
820 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpev"
821 : 0 : call xgBlock_hpev('v','u',subA%self,eigenvalues,vec%self,info)
822 : : case (EIGENSLK)
823 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using pheev"
824 : 0 : call xgScalapack_heev(scalapack,subA%self,eigenvalues)
825 : 0 : info = 0 ! No error code returned for the moment
826 : : case default
827 [ + - - - : 41192 : ABI_ERROR("Error for Eigen Solver HEEV")
- - ]
828 : : end select
829 : : else
830 : : ! Solve Hermitian general eigen problem only for first blockdim eigenvalues
831 : 0 : select case (eigenSolver)
832 : : case (EIGENVX)
833 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvx"
834 : : call xgBlock_hegvx(1,'v','i','u',subA%self,subB%self,0.d0,0.d0,1,blockdim,abstol,&
835 : 0 : eigenvalues,vec%self,info)
836 : : case (EIGENVD)
837 : : !ITEST
838 [ - + ]: 73872 : if (prtvol == 15015015) then
839 : 0 : write(902,*) 'eigenSolver', eigenSolver
840 : 0 : write(902,*) 'Using hegvd'
841 : 0 : nrows_B = rows(subB%self)
842 : 0 : ncols_B = cols(subB%self)
843 : 0 : write(902,*) 'computing cond for Hermitian pd matrix', nrows_B, ncols_B
844 : : !# Validation
845 : : !nrows_B = 2
846 : : !#
847 : 0 : call xgBlock_hermitian_pd_cond(subB%self, nrows_B, cond)
848 : 0 : write(902,*) 'cond(B)=', cond
849 : : !# Validation: should be 2.076578056
850 : 0 : flush(902)
851 : : end if
852 : : !ITEST
853 [ - + ]: 73872 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvd"
854 : 73872 : call xgBlock_hegvd(1,'v','u',subA%self,subB%self,eigenvalues,info)
855 : : case (EIGENV)
856 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegv"
857 : 0 : call xgBlock_hegv(1,'v','u',subA%self,subB%self,eigenvalues,info)
858 : : case (EIGENPVX)
859 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgvx"
860 : : call xgBlock_hpgvx(1,'v','i','u',subA%self,subB%self,0.d0,0.d0,1,blockdim,abstol,&
861 : 0 : eigenvalues,vec%self,info)
862 : : case (EIGENPVD)
863 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgvd"
864 : 0 : call xgBlock_hpgvd(1,'v','u',subA%self,subB%self,eigenvalues,vec%self,info)
865 : : case (EIGENPV)
866 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hpgv"
867 : 0 : call xgBlock_hpgv(1,'v','u',subA%self,subB%self,eigenvalues,vec%self,info)
868 : : case (EIGENSLK)
869 [ # # ]: 0 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using phegv"
870 : 0 : call xgScalapack_hegv(scalapack,subA%self,subB%self,eigenvalues)
871 : 0 : info = 0 ! No error code returned for the moment
872 : : case default
873 [ - + - - : 73872 : ABI_ERROR("Error for Eigen Solver HEGV")
- - - - ]
874 : : end select
875 : : end if
876 [ - + ]: 94468 : if ( eigenSolver == EIGENSLK ) then
877 : 0 : call xgScalapack_free(scalapack)
878 : : end if
879 : 94468 : tsec(2) = abi_wtime() - tsec(2)
880 [ - + ]: 94468 : if ( prtvol == 4 ) write(std_out,*) tsec(2)
881 : 94468 : call timab(tim_RR_diago,2,tsec)
882 : :
883 [ + - - + ]: 94468 : if ( eigenSolver == EIGENVX .or. EIGPACK(eigenSolver)) then
884 : 0 : call xg_free(subA)
885 : : end if
886 : 94468 : call xg_free(subB)
887 : :
888 : 94468 : call timab(tim_RR_gemm_2,1,tsec)
889 : : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_2)
890 : :
891 [ + - ]: 94468 : if ( info == 0 ) then
892 : 94468 : call xg_init(Xwork,space(X),spacedim,blockdim,me_g0=me_g0(X))
893 : 94468 : call xg_init(cprjXwork,space(cprjX),cprjdim,ncols_cprj,comm(cprjX))
894 : 94468 : call xg_setBlock(cprjXwork,cprj_workX,cprjdim,ncols_cprj)
895 : 94468 : call xgBlock_reshape_spinor(cprj_workX,cprj_workX_spinor,nspinor,COLS2ROWS)
896 : :
897 : 94468 : call xgBlock_setBlock(vec%self,Cwp,blockdim,blockdim)
898 : :
899 : : !/* Easy basic solution */
900 : : !/* Compute first part of X here */
901 : : !XWP (:,X+1:X+blockdim) = matmul(XWP (:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
902 : 94468 : call xgBlock_gemm('n','n',1.0d0,X,Cwp,0.d0,Xwork%self)
903 : 94468 : call xgBlock_copy(Xwork%self,X)
904 : :
905 : : !AXWP(:,X+1:X+blockdim) = matmul(AXWP(:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
906 : 94468 : call xgBlock_gemm('n','n',1.0d0,AX,Cwp,0.d0,Xwork%self)
907 : 94468 : call xgBlock_copy(Xwork%self,AX)
908 : :
909 : : !cprjXWP(:,X+1:X+blockdim) = matmul(cprjXWP(:,X+1:X+blockdim),vec(1:blockdim,1:blockdim))
910 : 94468 : call xgBlock_reshape_spinor(cprjX,cprjX_spinor,nspinor,COLS2ROWS)
911 [ + + ]: 94468 : if (ncols_cprj==nspinor*cols(X)) then
912 : 32234 : call xgBlock_gemm('n','n',1.0d0,cprjX_spinor,Cwp,0.d0,cprj_workX_spinor)
913 : : else
914 : 62234 : call xgBlock_zero(cprj_workX)
915 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjX_spinor,Cwp,cprj_workX_spinor,&
916 : 62234 : & xg_nonlop%me_band,blocksize=blockdim_cprj_/nspinor,comm=xg_nonlop%comm_band)
917 : : end if
918 : 94468 : call xgBlock_copy(cprj_workX,cprjX)
919 : :
920 [ + + ]: 94468 : if ( var /= VAR_X ) then
921 : : ! Cost to pay to avoid temporary array in xgemm
922 : 58504 : call xgBlock_cshift(vec%self,blockdim,1) ! Bottom 2*blockdim lines are now at the top
923 : 58504 : call xgBlock_setBlock(vec%self,Cwp,subdim-blockdim,blockdim)
924 : :
925 : : !XWP (:,P+1:P+blockdim) = matmul(XWP (:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
926 : 58504 : call xgBlock_gemm('n','n',1.0d0,WP,Cwp,0.d0,Xwork%self)
927 : 58504 : call xgBlock_copy(Xwork%self,P)
928 : :
929 : : !AXWP(:,P+1:P+blockdim) = matmul(AXWP(:,W+1:W+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
930 : 58504 : call xgBlock_gemm('n','n',1.0d0,AWP,Cwp,0.d0,Xwork%self)
931 : 58504 : call xgBlock_copy(Xwork%self,AP)
932 : :
933 : : !cprjXWP(:,P+1:P+blockdim) = matmul(cprjXWP(:,W+1:X+subdim-blockdim),vec(1:subdim-blockdim,1:blockdim))
934 : 58504 : call xgBlock_reshape_spinor(cprjWP,cprjWP_spinor,nspinor,COLS2ROWS)
935 [ + + ]: 58504 : if (ncols_cprj==cols(WP)) then
936 : 1863 : call xgBlock_gemm('n','n',1.0d0,cprjWP_spinor,Cwp,0.d0,cprj_workX_spinor)
937 : : else
938 : 56641 : call xgBlock_zero(cprj_workX)
939 [ + + ]: 56641 : if ( var==VAR_XW ) then
940 : 8655 : call xgBlock_reshape_spinor(cprjW,cprjW_spinor,nspinor,COLS2ROWS)
941 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjW_spinor,Cwp,cprj_workX_spinor,&
942 : 8655 : & xg_nonlop%me_band,blocksize=blockdim_cprj_/nspinor,comm=xg_nonlop%comm_band)
943 : : else if ( var==VAR_XWP ) then
944 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjWP_spinor,Cwp,cprj_workX_spinor,&
945 : 47986 : & xg_nonlop%me_band,blocksize=blockdim_cprj_/nspinor,comm=xg_nonlop%comm_band)
946 : : else
947 : : ABI_ERROR('not implemented')
948 : : end if
949 : : end if
950 : 58504 : call xgBlock_copy(cprj_workX,cprjP)
951 : :
952 : : !/* Maybe faster solution
953 : : ! * Sum previous contribution plus P direction
954 : : ! */
955 : 58504 : call xgBlock_add(X,P)
956 : 58504 : call xgBlock_add(AX,AP)
957 : 58504 : call xgBlock_add(cprjX,cprjP)
958 : : end if
959 : :
960 : : end if
961 : :
962 : 94468 : call timab(tim_RR_gemm_2,2,tsec)
963 : : ABI_NVTX_END_RANGE()
964 : :
965 : : ! Doing free on an already free object does not do anything
966 : 94468 : call xg_free(vec)
967 : 94468 : call xg_free(subA)
968 : 94468 : call xg_free(Xwork)
969 : 94468 : call xg_free(cprjXwork)
970 : :
971 : 94468 : call timab(timer , 2, tsec)
972 : :
973 [ + + ]: 283404 : end subroutine xg_RayleighRitz_cprj
974 : : !!***
975 : :
976 : : end module m_xg_ortho_RR
977 : : !!***
|
