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 1429086 : 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 714543 : call timab(timer,1,tsec)
101 :
102 714543 : call xgBlock_check(X,BX)
103 714543 : if (present(AX)) then
104 714543 : call xgBlock_check(X,AX)
105 : end if
106 :
107 714543 : if (space(X)/=SPACE_CR) then
108 657756 : space_buf = SPACE(X)
109 : else
110 56787 : space_buf = SPACE_R
111 : end if
112 714543 : 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 714543 : call xgBlock_zero_im_g0(X)
116 714543 : call xgBlock_zero_im_g0(BX)
117 714543 : if (present(AX)) then
118 714543 : call xgBlock_zero_im_g0(AX)
119 : end if
120 :
121 : ! Compute X^TBX
122 714543 : call xgBlock_gemm('t','n',1.d0,X,BX,0.d0,buffer%self,comm=comm(X))
123 :
124 : ! Compute Cholesky decomposition (Upper part)
125 714543 : call xgBlock_potrf(buffer%self,'u',info)
126 :
127 714543 : 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 714535 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,X)
136 :
137 : ! Solve BYU=BX
138 714535 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,BX)
139 :
140 714535 : if (present(AX)) then
141 : ! Solve AYU=AX
142 714535 : call xgBlock_trsm('r','u',buffer%normal,'n',1.d0,buffer%self,AX)
143 : end if
144 :
145 714535 : call xg_free(buffer)
146 :
147 : ABI_NVTX_END_RANGE()
148 714535 : call timab(timer,2,tsec)
149 :
150 714543 : 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 723559 : 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 723559 : call timab(timer , 1, tsec)
297 :
298 723559 : solve_ax_bx_ = .false.
299 723559 : if (present(solve_ax_bx)) then
300 9024 : solve_ax_bx_ = solve_ax_bx
301 : end if
302 :
303 723559 : var = VAR_X
304 723559 : call xgBlock_check(X,AX)
305 723559 : call xgBlock_check(X,BX)
306 723559 : if (.not.solve_ax_bx_) then
307 : eigenSolver = EIGENEVD
308 : else
309 9024 : eigenSolver = EIGENVD
310 : end if
311 723559 : spacedim = rows(X)
312 723559 : blockdim = cols(X)
313 723559 : space_buf = space(X)
314 723559 : if (space(X) == space_CR) then
315 60323 : space_buf = space_R
316 : end if
317 723559 : spacecom = comm(X)
318 723559 : subdim = blockdim
319 :
320 723559 : if (present(XW)) then
321 556994 : if (solve_ax_bx_) then
322 0 : ABI_ERROR('solve_ax_bx is not implemented in that case')
323 : end if
324 556994 : var = VAR_XW
325 556994 : eigenSolver = EIGENVD
326 556994 : call xgBlock_check(X,AW)
327 556994 : call xgBlock_check(X,BW)
328 556994 : call xgBlock_check(X,P)
329 556994 : call xgBlock_check(X,AP)
330 556994 : call xgBlock_check(X,BP)
331 556994 : call xgBlock_check(X,XW,fact_col=2)
332 556994 : call xgBlock_check(X,WP,fact_col=2)
333 556994 : call xgBlock_check(X,AWP,fact_col=2)
334 556994 : call xgBlock_check(X,BWP,fact_col=2)
335 556994 : subdim = 2*blockdim
336 556994 : if (present(XWP)) then
337 423092 : var = VAR_XWP
338 423092 : call xgBlock_check(X,XWP,fact_col=3)
339 423092 : 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 723559 : abstol = 0d0 ; if ( present(tolerance) ) abstol = tolerance
353 :
354 723559 : call xg_init(subA,space_buf,subdim,subdim,spacecom,gpu_option=gpu_option)
355 723559 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
356 566018 : 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 723559 : else if ( EIGPACK(eigenSolver) ) then
362 0 : call xg_init(vec,space_buf,subdim,subdim,gpu_option=gpu_option)
363 : else
364 723559 : 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 723559 : 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 723559 : if (var == VAR_X) call xgBlock_zero_im_g0(X)
378 723559 : if (var == VAR_XW) call xgBlock_zero_im_g0(XW)
379 723559 : if (var == VAR_XWP) call xgBlock_zero_im_g0(XWP)
380 723559 : call xgBlock_zero_im_g0(AX)
381 723559 : call xgBlock_zero_im_g0(BX)
382 :
383 723559 : call xg_setBlock(subA,subsub,blockdim,blockdim)
384 723559 : call xgBlock_gemm('t','n',1.0d0,X,AX,0.d0,subsub,comm=spacecom)
385 :
386 723559 : if ( solve_ax_bx_ .or. var /= VAR_X ) then
387 566018 : call xg_setBlock(subB,subsub,blockdim,blockdim)
388 566018 : call xgBlock_gemm('t','n',1.0d0,X,BX,0.d0,subsub,comm=spacecom)
389 : endif
390 :
391 723559 : 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 556994 : call xgBlock_zero_im_g0(AW)
395 556994 : call xgBlock_zero_im_g0(BW)
396 :
397 : ! subA
398 556994 : call xg_setBlock(subA,subsub,2*blockdim,blockdim,fcol=blockdim+1)
399 556994 : call xgBlock_gemm('t','n',1.0d0,XW,AW,0.d0,subsub,comm=spacecom)
400 :
401 : ! subB
402 556994 : call xg_setBlock(subB,subsub,2*blockdim,blockdim,fcol=blockdim+1)
403 556994 : call xgBlock_gemm('t','n',1.0d0,XW,BW,0.d0,subsub,comm=spacecom)
404 :
405 : end if
406 :
407 556994 : 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 423092 : call xgBlock_zero_im_g0(AP)
411 423092 : call xgBlock_zero_im_g0(BP)
412 :
413 : ! subA
414 423092 : call xg_setBlock(subA,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
415 423092 : call xgBlock_gemm('t','n',1.0d0,XWP,AP,0.d0,subsub,comm=spacecom)
416 :
417 : ! subB
418 423092 : call xg_setBlock(subB,subsub,3*blockdim,blockdim,fcol=2*blockdim+1)
419 423092 : call xgBlock_gemm('t','n',1.0d0,XWP,BP,0.d0,subsub,comm=spacecom)
420 :
421 : end if
422 :
423 723559 : call timab(tim_RR_gemm_1,2,tsec)
424 : ABI_NVTX_END_RANGE()
425 :
426 723559 : 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 723559 : call timab(tim_RR_diago,1,tsec)
434 723559 : tsec(2) = abi_wtime()
435 723559 : if ( .not.solve_ax_bx_ .and. var == VAR_X ) then
436 : ABI_NVTX_START_RANGE(NVTX_RR_HEEV)
437 : ! Solve Hermitian eigen problem
438 157541 : select case (eigenSolver)
439 : case (EIGENEVD)
440 157541 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using heevd"
441 157541 : 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 315082 : 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 566018 : eigenvalues,vec%self,info)
466 : case (EIGENVD)
467 566018 : if ( prtvol == 4 ) write(std_out,'(A,1x)',advance="no") "Using hegvd"
468 566018 : 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 566018 : 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 723559 : tsec(2) = abi_wtime() - tsec(2)
494 723559 : if ( prtvol == 4 ) write(std_out,*) tsec(2)
495 :
496 723559 : call timab(tim_RR_diago,2,tsec)
497 : ABI_NVTX_END_RANGE()
498 :
499 723559 : if ( eigenSolver == EIGENVX .or. EIGPACK(eigenSolver)) then
500 0 : call xg_free(subA)
501 : end if
502 723559 : call xg_free(subB)
503 :
504 723559 : call timab(tim_RR_gemm_2,1,tsec)
505 : ABI_NVTX_START_RANGE(NVTX_RR_GEMM_2)
506 :
507 : !FIXME Avoid those transfers
508 723559 : if ( info == 0 ) then
509 723559 : 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 723559 : call xgBlock_setBlock(vec%self,Cwp,blockdim,blockdim)
516 723559 : call xgBlock_gemm('n','n',1.0d0,X,Cwp,0.d0,subB%self)
517 723559 : 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 723559 : call xgBlock_gemm('n','n',1.0d0,AX,Cwp,0.d0,subB%self)
521 723559 : 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 723559 : call xgBlock_gemm('n','n',1.0d0,BX,Cwp,0.d0,subB%self)
525 723559 : call xgBlock_copy(subB%self,BX)
526 :
527 723559 : if ( var /= VAR_X ) then
528 : ! Cost to pay to avoid temporary array in xgemm
529 556994 : if(gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_from_gpu(vec%self) !FIXME Avoid that transfer
530 556994 : call xgBlock_cshift(vec%self,blockdim,1) ! Bottom 2*blockdim lines are now at the top
531 556994 : if(gpu_option==ABI_GPU_OPENMP) call xgBlock_copy_to_gpu(vec%self) !FIXME Avoid that transfer
532 556994 : 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 556994 : call xgBlock_gemm('n','n',1.0d0,WP,Cwp,0.d0,subB%self)
536 556994 : 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 556994 : call xgBlock_gemm('n','n',1.0d0,AWP,Cwp,0.d0,subB%self)
540 556994 : 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 556994 : call xgBlock_gemm('n','n',1.0d0,BWP,Cwp,0.d0,subB%self)
544 556994 : call xgBlock_copy(subB%self,BP)
545 :
546 : !/* Maybe faster solution
547 : ! * Sum previous contribution plus P direction
548 : ! */
549 556994 : call xgBlock_add(X,P)
550 556994 : call xgBlock_add(AX,AP)
551 556994 : call xgBlock_add(BX,BP)
552 : end if
553 : end if
554 :
555 723559 : call timab(tim_RR_gemm_2,2,tsec)
556 :
557 : ! Doing free on an already free object does not do anything
558 723559 : call xg_free(vec)
559 723559 : call xg_free(subA)
560 723559 : 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 723559 : call timab(timer , 2, tsec)
570 :
571 2170677 : 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 : !!***
|
