Line data Source code
1 : !!****f* ABINIT/m_chebfi2
2 : !! NAME
3 : !! m_chebfi2
4 : !!
5 : !! FUNCTION
6 : !! This module contains the types and routines used to apply the
7 : !! Chebyshev filtering method (2021 implementation using xG abstraction layer)
8 : !! It mainly defines a 'chebfi' datatypes and associated methods.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2018-2026 ABINIT group (BS, L. Baguet, IML)
12 : !! This file is distributed under the terms of the
13 : !! gnu general public license, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !! for the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : ! nvtx related macro definition
26 : #include "nvtx_macros.h"
27 :
28 : module m_chebfi2
29 :
30 : use defs_basis
31 : use defs_abitypes
32 : use m_abicore
33 : use m_errors
34 : use m_time, only : timab
35 :
36 : use m_cgtools
37 : use m_xg
38 : use m_xgTransposer
39 : use m_xg_ortho_RR
40 : use m_polynomial_filter
41 :
42 : use m_xmpi
43 : use m_xomp
44 : #ifdef HAVE_OPENMP
45 : use omp_lib
46 : #endif
47 : use, intrinsic :: iso_c_binding, only: c_size_t
48 :
49 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
50 : use m_gpu_toolbox, only : CPU_DEVICE_ID, gpu_device_synchronize
51 : #endif
52 :
53 : #if defined(HAVE_GPU_MARKERS)
54 : use m_nvtx_data
55 : #endif
56 :
57 : implicit none
58 :
59 : private
60 :
61 : !Several (private) parameters
62 : !-------------------------------------------------
63 :
64 : integer, parameter :: tim_init = 1751
65 : integer, parameter :: tim_free = 1752
66 : ! 1753 is used by chebfi2_nonlop
67 : integer, parameter :: tim_getAX_BX = 1754
68 : integer, parameter :: tim_invovl = 1755
69 : integer, parameter :: tim_residu = 1756
70 : integer, parameter :: tim_RR = 1757
71 : integer, parameter :: tim_transpose = 1758
72 : integer, parameter :: tim_RR_q = 1759
73 : integer, parameter :: tim_postinvovl = 1760
74 : integer, parameter :: tim_swap = 1761
75 : integer, parameter :: tim_amp_f = 1762
76 : integer, parameter :: tim_oracle = 1763
77 : integer, parameter :: tim_barrier = 1764
78 : integer, parameter :: tim_copy = 1765
79 :
80 : !Public 'chebfi' datatype
81 : !-------------------------------------------------
82 :
83 : type, public :: chebfi_t
84 : integer :: space
85 : integer :: spacedim ! Space dimension for one vector
86 : integer :: total_spacedim ! Maybe not needed
87 : integer :: neigenpairs ! Number of eigen values/vectors we want
88 : integer :: ndeg_filter ! Degree of the polynomial filter
89 : integer :: nbdbuf ! Number of bands in the buffer
90 : integer :: spacecom ! Communicator for MPI
91 : integer :: oracle ! Option to compute ndeg_filter from residuals
92 : real(dp) :: tolerance ! Tolerance on the residu to stop the minimization
93 : real(dp) :: ecut ! Ecut for Chebfi oracle
94 : real(dp) :: oracle_factor ! factor used to decrease residuals
95 : real(dp) :: oracle_min_occ ! threshold on occupancies used for nbdbuf=-101
96 :
97 : integer :: paral_kgb ! MPI parallelization variables
98 : integer :: bandpp
99 : integer :: comm_cols
100 : integer :: comm_rows
101 : integer :: me_g0
102 : integer :: me_g0_fft
103 :
104 : logical :: from_linalg ! Transposer allocates representation ColsRows if true, Linalg if false
105 : logical :: paw
106 : integer :: eigenProblem !1 (A*x = (lambda)*B*x), 2 (A*B*x = (lambda)*x), 3 (B*A*x = (lambda)*x)
107 :
108 : ! when GPU with Kokkos is enabled, currently OpenMP is not fully supported, abinit is launched
109 : ! with OMP_NUM_THREADS=1, but we may locally increase the number of OpenMP threads
110 : ! wherever it is safe to do; in that case we use gpu_kokkos_nthrd to specify
111 : ! the number of OpenMP threads. This value is controlled by dtset variable
112 : ! dtset%gpu_kokkos_nthrd
113 : integer :: gpu_option
114 : integer :: gpu_kokkos_nthrd = 1 ! only used if GPU Kokkos is enabled, number of OpenMP threads used
115 : integer :: gpu_thread_limit = 1 ! only used if GPU is enabled, max number of OpenMP threads used in sensitive areas
116 :
117 : !ARRAYS
118 : type(xgBlock_t) :: X
119 :
120 : type(xg_t) :: X_NP
121 : type(xgBlock_t) :: X_next
122 : type(xgBlock_t) :: X_prev
123 :
124 : type(xg_t) :: AX
125 : type(xg_t) :: BX
126 : type(xg_t) :: xAXColsRows_W ! only used if from_linalg is false
127 : type(xg_t) :: xBXColsRows_W ! only used if from linalg is false
128 :
129 : type(xgBlock_t) :: xXColsRows
130 : type(xgBlock_t) :: xAXColsRows
131 : type(xgBlock_t) :: xBXColsRows
132 :
133 : type(xgTransposer_t) :: xgTransposerX
134 : type(xgTransposer_t) :: xgTransposerAX
135 : type(xgTransposer_t) :: xgTransposerBX
136 :
137 : type(xgBlock_t) :: eigenvalues
138 :
139 : !SWAP POINTERS
140 : type(xgBlock_t) :: X_swap
141 : type(xgBlock_t) :: AX_swap
142 : type(xgBlock_t) :: BX_swap
143 :
144 : end type chebfi_t
145 :
146 : ! Partition Info for column vectors X=(X_lock X_active)
147 : ! Function: handles MPI distribution variables and communicators
148 : !-------------------------------------------------
149 : type, private :: bandPartitionInfo_t
150 :
151 : integer :: n_locked ! MPI size in linalg for locked
152 : integer :: n_active ! MPI size in linalg for active
153 : integer :: comm_active ! communicator for processes treating active bands
154 :
155 : end type bandPartitionInfo_t
156 :
157 : ! Partition Data for column vectors X=(X_lock X_active)
158 : ! Function: handles data (e.g. X, AX, BX) in MPI distributions
159 : !-------------------------------------------------
160 : type, private :: bandPartitionData_t
161 :
162 : type(xgBlock_t) :: linalg_active ! active vectors in linalg MPI distribution
163 : type(xgBlock_t) :: linalg_locked ! locked vectors in linalg MPI distribution
164 : type(xgBlock_t) :: colsrows_active ! active vectors in colsrows MPI distribution
165 : type(xgTransposer_t) :: transposer_active ! transposer to switch between MPI distrs for active only
166 :
167 : end type bandPartitionData_t
168 :
169 : !Public methods associated to 'chebfi' datatype
170 : !-------------------------------------------------
171 : public :: chebfi_init
172 : public :: chebfi_free
173 : public :: chebfi_memInfo
174 : public :: chebfi_run
175 : public :: chebfi_runSlice
176 : public :: chebfi_runSubspaceIteration
177 : public :: chebfi_runSubspaceIterationDummy
178 : public :: chebfi_rayleighRitzQuotients
179 : public :: chebfi_computeNextOrderChebfiPolynom
180 : public :: chebfi_swapInnerBuffers
181 :
182 : CONTAINS !========================================================================================
183 : !!***
184 :
185 : !!****f* m_chebfi2/chebfi_init
186 : !! NAME
187 : !! chebfi_init
188 : !!
189 : !! FUNCTION
190 : !! Initialize a 'chebfi' datastructure.
191 : !!
192 : !! INPUTS
193 : !! bandpp= number of 'bands' handled by a processor
194 : !! eigenProblem= type of eigenpb: 1 (A*x = (lambda)*B*x), 2 (A*B*x = (lambda)*x), 3 (B*A*x = (lambda)*x)
195 : !! me_g0= 1 if this processors treats G=0, 0 otherwise
196 : !! me_g0_fft= 1 if this processors treats G=0 in FFT, 0 otherwise
197 : !! neigenpairs= number of requested eigenvectors/eigenvalues
198 : !! ndeg_filter= polynomial degree of the Chebyshev filter (.i.e. number of H applications)
199 : !! comm_rows= "rows" communicator
200 : !! comm_cols= "cols" communicator
201 : !! paral_kgb= flag controlling (k,g,bands) parallelization
202 : !! space= defines in which space we are (columns, rows, etc.)
203 : !! spacecom= MPI communicator
204 : !! spacedim= space dimension for one vector
205 : !! paw= flag. TRUE if current calculation ses the PAW approach
206 : !! ecut= plane-wave cut-off energy
207 : !! tolerance= tolerance criterion on the residu to stop the minimization
208 : !! nbdbuf= number of bands in the buffer
209 : !! oracle= option compute ndeg_filter from residuals
210 : !! oracle_factor= factor used to decrease residuals
211 : !! oracle_min_occ= threshold on occupancies used for nbdbuf=-101
212 : !! gpu_option= flag. Enable GPU if true
213 : !! gpu_kokkos_nthrd= number of OpenMP offloaded threads used
214 : !! gpu_thread_limit= maximum number of OpenMP offloaded threads
215 : !! from_linalg= flag. Transposer allocates representation ColsRows if true, Linalg if false
216 : !!
217 : !! OUTPUT
218 : !!
219 : !! SIDE EFFECTS
220 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
221 : !!
222 : !! SOURCE
223 :
224 9042 : subroutine chebfi_init(chebfi,neigenpairs,spacedim,tolerance,ecut,paral_kgb,bandpp, &
225 : ndeg_filter,nbdbuf,space,eigenProblem,spacecom,me_g0,me_g0_fft,paw,comm_rows,comm_cols, &
226 : oracle,oracle_factor,oracle_min_occ,gpu_option,gpu_kokkos_nthrd,gpu_thread_limit,from_linalg)
227 :
228 : ! Arguments ------------------------------------
229 : integer , intent(in ) :: bandpp
230 : integer , intent(in ) :: eigenProblem
231 : integer , intent(in ) :: me_g0
232 : integer , intent(in ) :: me_g0_fft
233 : integer , intent(in ) :: neigenpairs
234 : integer , intent(in ) :: ndeg_filter
235 : integer , intent(in ) :: nbdbuf
236 : integer , intent(in ) :: comm_cols
237 : integer , intent(in ) :: comm_rows
238 : integer , intent(in ) :: paral_kgb
239 : integer , intent(in ) :: space
240 : integer , intent(in ) :: spacecom
241 : integer , intent(in ) :: spacedim
242 : integer , intent(in ) :: gpu_option
243 : integer , intent(in ) :: oracle
244 : logical , intent(in ) :: paw
245 : real(dp) , intent(in ) :: ecut
246 : real(dp) , intent(in ) :: tolerance
247 : real(dp) , intent(in ) :: oracle_factor
248 : real(dp) , intent(in ) :: oracle_min_occ
249 : type(chebfi_t), intent(inout) :: chebfi
250 : integer , intent(in ), optional :: gpu_kokkos_nthrd
251 : integer , intent(in ), optional :: gpu_thread_limit
252 : logical , intent(in ), optional :: from_linalg
253 :
254 : ! Local variables-------------------------------
255 : real(dp) :: tsec(2)
256 :
257 : ! *********************************************************************
258 :
259 9042 : call timab(tim_init,1,tsec)
260 :
261 9042 : chebfi%space = space
262 9042 : chebfi%neigenpairs = neigenpairs
263 9042 : chebfi%spacedim = spacedim
264 9042 : if (tolerance > 0.0) then
265 3010 : chebfi%tolerance = tolerance
266 : else
267 6032 : chebfi%tolerance = 1.0e-20
268 : end if
269 9042 : chebfi%ecut = ecut
270 9042 : chebfi%paral_kgb = paral_kgb
271 9042 : chebfi%comm_cols = comm_cols
272 9042 : chebfi%bandpp = bandpp
273 9042 : chebfi%comm_rows = comm_rows
274 9042 : chebfi%ndeg_filter = ndeg_filter
275 9042 : chebfi%nbdbuf = nbdbuf
276 9042 : chebfi%spacecom = spacecom
277 9042 : chebfi%eigenProblem = eigenProblem
278 9042 : chebfi%me_g0 = me_g0
279 9042 : chebfi%me_g0_fft = me_g0_fft
280 9042 : chebfi%paw = paw
281 9042 : chebfi%gpu_option = gpu_option
282 9042 : chebfi%oracle = oracle
283 9042 : chebfi%oracle_factor = oracle_factor
284 9042 : chebfi%oracle_min_occ = oracle_min_occ
285 :
286 9042 : chebfi%gpu_kokkos_nthrd = 1
287 9042 : if (present(gpu_kokkos_nthrd)) chebfi%gpu_kokkos_nthrd = gpu_kokkos_nthrd
288 9042 : chebfi%gpu_thread_limit = 0
289 9042 : if (present(gpu_thread_limit)) chebfi%gpu_thread_limit = gpu_thread_limit
290 9042 : chebfi%from_linalg = .true.
291 9042 : if (present(from_linalg)) chebfi%from_linalg = from_linalg
292 :
293 9042 : call chebfi_allocateAll(chebfi)
294 :
295 9042 : call timab(tim_init,2,tsec)
296 :
297 9042 : end subroutine chebfi_init
298 : !!***
299 :
300 : !----------------------------------------------------------------------
301 :
302 : !!****f* m_chebfi2/chebfi_allocateAll
303 : !! NAME
304 : !! chebfi_allocateAll
305 : !!
306 : !! FUNCTION
307 : !! Allocate all memory spaces in a 'chebfi' datastructure.
308 : !!
309 : !! INPUTS
310 : !!
311 : !! OUTPUT
312 : !!
313 : !! SIDE EFFECTS
314 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
315 : !!
316 : !! SOURCE
317 :
318 9042 : subroutine chebfi_allocateAll(chebfi)
319 :
320 : ! Arguments ------------------------------------
321 : type(chebfi_t) , intent(inout) :: chebfi
322 :
323 : ! Local variables-------------------------------
324 : ! scalars
325 : integer :: neigenpairs
326 : integer :: space
327 : integer :: spacedim
328 : integer :: total_spacedim, ierr
329 :
330 : ! *********************************************************************
331 :
332 9042 : space = chebfi%space
333 9042 : spacedim = chebfi%spacedim
334 9042 : neigenpairs = chebfi%neigenpairs
335 :
336 9042 : call chebfi_free(chebfi)
337 :
338 9042 : if (chebfi%paral_kgb == 0) then
339 2722 : chebfi%total_spacedim = spacedim
340 2722 : call xg_init(chebfi%X_NP,space,spacedim,2*neigenpairs,chebfi%spacecom,me_g0=chebfi%me_g0,gpu_option=chebfi%gpu_option) !regular arrays
341 2722 : call xg_setBlock(chebfi%X_NP, chebfi%X_next,spacedim, neigenpairs)
342 2722 : call xg_setBlock(chebfi%X_NP, chebfi%X_prev,spacedim, neigenpairs, fcol=neigenpairs+1)
343 : else
344 6320 : if (chebfi%from_linalg) then
345 6304 : total_spacedim = spacedim
346 6304 : call xmpi_sum(total_spacedim,chebfi%spacecom,ierr)
347 6304 : chebfi%total_spacedim = total_spacedim
348 : else
349 16 : chebfi%total_spacedim = chebfi%spacedim
350 : end if
351 : call xg_init(chebfi%X_NP,space,chebfi%total_spacedim,2*chebfi%bandpp,chebfi%spacecom,me_g0=chebfi%me_g0_fft,&
352 6320 : & gpu_option=chebfi%gpu_option) !transposed arrays
353 6320 : call xg_setBlock(chebfi%X_NP, chebfi%X_next, chebfi%total_spacedim, chebfi%bandpp)
354 6320 : call xg_setBlock(chebfi%X_NP, chebfi%X_prev, chebfi%total_spacedim, chebfi%bandpp, fcol=chebfi%bandpp+1)
355 : end if
356 :
357 : !transposer will handle these arrays automatically
358 9042 : if (chebfi%from_linalg) then
359 9024 : call xg_init(chebfi%AX,space,spacedim,neigenpairs,chebfi%spacecom,me_g0=chebfi%me_g0,gpu_option=chebfi%gpu_option)
360 9024 : call xg_init(chebfi%BX,space,spacedim,neigenpairs,chebfi%spacecom,me_g0=chebfi%me_g0,gpu_option=chebfi%gpu_option)
361 : else
362 : call xg_init(chebfi%xAXColsRows_W,space,chebfi%total_spacedim,chebfi%bandpp,chebfi%comm_rows,&
363 18 : me_g0=chebfi%me_g0_fft,gpu_option=chebfi%gpu_option)
364 : call xg_init(chebfi%xBXColsRows_W,space,chebfi%total_spacedim,chebfi%bandpp,chebfi%comm_rows,&
365 18 : me_g0=chebfi%me_g0_fft,gpu_option=chebfi%gpu_option)
366 18 : chebfi%xAXColsRows = chebfi%xAXColsRows_W%self
367 18 : chebfi%xBXColsRows = chebfi%xBXColsRows_W%self
368 : end if
369 :
370 9042 : end subroutine chebfi_allocateAll
371 : !!***
372 :
373 : !----------------------------------------------------------------------
374 :
375 : !!****f* m_chebfi2/chebfi_free
376 : !! NAME
377 : !! chebfi_free
378 : !!
379 : !! FUNCTION
380 : !! Destroy a 'chebfi' datastructure.
381 : !!
382 : !! INPUTS
383 : !!
384 : !! OUTPUT
385 : !! arraymem(2)= memory information
386 : !!
387 : !! SIDE EFFECTS
388 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
389 : !!
390 : !! SOURCE
391 :
392 18084 : subroutine chebfi_free(chebfi)
393 :
394 : !Arguments ------------------------------------
395 : type(chebfi_t) , intent(inout) :: chebfi
396 :
397 : ! *********************************************************************
398 :
399 18084 : call xg_free(chebfi%X_NP)
400 :
401 18084 : if (chebfi%from_linalg) then
402 18048 : call xg_free(chebfi%AX)
403 18048 : call xg_free(chebfi%BX)
404 : else
405 36 : call xg_free(chebfi%xAXColsRows_W)
406 36 : call xg_free(chebfi%xBXColsRows_W)
407 : end if
408 :
409 18084 : end subroutine chebfi_free
410 : !!***
411 :
412 : !----------------------------------------------------------------------
413 :
414 : !!****f* m_chebfi2/chebfi_memInfo
415 : !! NAME
416 : !! chebfi_memInfo
417 : !!
418 : !! FUNCTION
419 : !! Provides memory information about a 'chebfi' datastructure.
420 : !!
421 : !! INPUTS
422 : !! bandpp= number of 'bands' handled by a processor
423 : !! neigenpairs= number of requested eigenvectors/eigenvalues
424 : !! paral_kgb= flag controlling (k,g,bands) parallelization
425 : !! space= defines in which space we are (columns, rows, etc.)
426 : !! spacedim= dimension of MPI communicator
427 : !! total_spacedim= size of global KGB communicator (typically 'banspinorfft' comm.)
428 : !!
429 : !! OUTPUT
430 : !! arraymem(2)= memory information
431 : !!
432 : !! SIDE EFFECTS
433 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
434 : !!
435 : !! SOURCE
436 :
437 0 : function chebfi_memInfo(neigenpairs,spacedim,space,paral_kgb,total_spacedim,bandpp) result(arraymem)
438 :
439 : !Arguments ------------------------------------
440 : integer, intent(in ) :: bandpp
441 : integer, intent(in ) :: neigenpairs
442 : integer, intent(in ) :: paral_kgb
443 : integer, intent(in ) :: space
444 : integer, intent(in ) :: spacedim
445 : integer, intent(in ) :: total_spacedim
446 :
447 : !Local variables-------------------------------
448 : !scalars
449 : integer(kind=c_size_t) :: memX
450 : integer(kind=c_size_t) :: memX_next
451 : integer(kind=c_size_t) :: memX_prev
452 : integer(kind=c_size_t) :: memAX
453 : integer(kind=c_size_t) :: memBX
454 : !Transposer variables
455 : integer(kind=c_size_t) :: memX_CR
456 : integer(kind=c_size_t) :: memAX_CR
457 : integer(kind=c_size_t) :: memBX_CR
458 : integer(kind=c_size_t) :: mem_sendrecv_CR
459 : !chebfi_rayleighRitz function variables
460 : integer(kind=c_size_t) :: memA_und_X
461 : integer(kind=c_size_t) :: memB_und_X
462 : integer(kind=c_size_t) :: memEigenvalues
463 : integer(kind=c_size_t) :: cplx
464 : !arrays
465 : integer(kind=c_size_t) :: arraymem(2)
466 :
467 : ! *********************************************************************
468 0 : cplx = 1
469 0 : if ( space == SPACE_C ) cplx = 2 !for now only complex
470 :
471 : !Permanent in chebfi
472 0 : memX = int(cplx,c_size_t) * kind(1.d0) * spacedim * neigenpairs
473 :
474 0 : if (paral_kgb == 0) then
475 : memX_next = int(cplx,c_size_t) * kind(1.d0) * spacedim * neigenpairs
476 : memX_prev = int(cplx,c_size_t) * kind(1.d0) * spacedim * neigenpairs
477 : else
478 0 : memX_next = int(cplx,c_size_t) * kind(1.d0) * total_spacedim * bandpp
479 0 : memX_prev = int(cplx,c_size_t) * kind(1.d0) * total_spacedim * bandpp
480 : end if
481 :
482 0 : memAX = int(cplx,c_size_t) * kind(1.d0) * spacedim * neigenpairs
483 0 : memBX = int(cplx,c_size_t) * kind(1.d0) * spacedim * neigenpairs
484 :
485 : !Transposer colrow array
486 0 : if (paral_kgb == 1) then
487 0 : memX_CR = int(cplx,c_size_t) * kind(1.d0) * total_spacedim * bandpp
488 0 : memAX_CR = int(cplx,c_size_t) * kind(1.d0) * total_spacedim * bandpp
489 0 : memBX_CR = int(cplx,c_size_t) * kind(1.d0) * total_spacedim * bandpp
490 0 : mem_sendrecv_CR = int(cplx,c_size_t) * kind(1.d0) * total_spacedim * bandpp
491 : else
492 : memX_CR = 0
493 : memAX_CR = 0
494 : memBX_CR = 0
495 : mem_sendrecv_CR = 0
496 : end if
497 :
498 : !chebfi_rayleighRitz function variables
499 0 : memA_und_X = int(cplx,c_size_t) * kind(1.d0) * neigenpairs * neigenpairs
500 0 : memB_und_X = int(cplx,c_size_t) * kind(1.d0) * neigenpairs * neigenpairs
501 0 : memEigenvalues = int(kind(1.d0),c_size_t) * neigenpairs
502 :
503 : arraymem(1) = memX + memX_next + memX_prev + &
504 0 : memAX + memBX + memX_CR + memAX_CR + memBX_CR + mem_sendrecv_CR
505 0 : arraymem(2) = memA_und_X + memB_und_X + memEigenvalues
506 :
507 0 : end function chebfi_memInfo
508 : !!***
509 :
510 : !----------------------------------------------------------------------
511 :
512 : !!****f* m_chebfi2/chebfi_run
513 : !! NAME
514 : !! chebfi_run
515 : !!
516 : !! FUNCTION
517 : !! Apply the Chebyshev Filtering algorithm on a set of vectors.
518 : !!
519 : !! INPUTS
520 : !! getAX_BX= pointer to the function giving A|X> and B|X>
521 : !! A is typically the Hamiltonian H, and B the overlap operator S
522 : !! getBm1X= pointer to the function giving B^-1|X>
523 : !! B is typically the overlap operator S
524 : !! nspinor= number of spinorial components of the wavefunctions
525 : !!
526 : !! OUTPUT
527 : !!
528 : !! SIDE EFFECTS
529 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
530 : !! eigen= Full eigenvalues (initial values on entry)
531 : !! residu= residuals, i.e. norm of (A-lambdaB)|X>
532 : !! X0= Full set of vectors (initial values on entry) distributed in Linalg representation
533 : !!
534 : !! SOURCE
535 :
536 9024 : subroutine chebfi_run(chebfi,X0,getAX_BX,getBm1X,eigen,occ,residu,nspinor)
537 :
538 : !Arguments ------------------------------------
539 : type(chebfi_t) , intent(inout) :: chebfi
540 : integer, intent(in) :: nspinor
541 : type(xgBlock_t), intent(inout) :: X0
542 : type(xgBlock_t), intent(inout) :: eigen
543 : type(xgBlock_t), intent(in) :: occ
544 : type(xgBlock_t), intent(inout) :: residu
545 : interface
546 : subroutine getAX_BX(X,AX,BX)
547 : use m_xg, only : xgBlock_t
548 : type(xgBlock_t), intent(inout) :: X
549 : type(xgBlock_t), intent(inout) :: AX
550 : type(xgBlock_t), intent(inout) :: BX
551 : end subroutine getAX_BX
552 : end interface
553 : interface
554 : subroutine getBm1X(X,Bm1X)
555 : use m_xg, only : xgBlock_t
556 : type(xgBlock_t), intent(inout) :: X
557 : type(xgBlock_t), intent(inout) :: Bm1X
558 : end subroutine getBm1X
559 : end interface
560 :
561 : !Local variables-------------------------------
562 : !scalars
563 : integer :: spacedim
564 : integer :: space_res
565 : integer :: neigenpairs,bandpp
566 : integer :: ndeg_filter,ndeg_filter_max
567 : integer :: ideg, ierr
568 : real(dp) :: tolerance
569 : real(dp) :: maxeig, maxeig_global
570 : real(dp) :: mineig, mineig_global
571 : real(dp) :: lambda_minus
572 : real(dp) :: lambda_plus
573 : real(dp) :: one_over_r
574 : real(dp) :: two_over_r
575 : real(dp) :: center
576 : real(dp) :: radius
577 : type(xg_t) :: DivResults
578 : !arrays
579 : real(dp) :: tsec(2)
580 : !Pointers similar to old Chebfi
581 9024 : integer,allocatable :: ndeg_filter_bands(:) !Oracle variable
582 :
583 : ! *********************************************************************
584 :
585 : ! call timab(tim_run,1,tsec)
586 :
587 9024 : spacedim = chebfi%spacedim
588 9024 : neigenpairs = chebfi%neigenpairs
589 9024 : bandpp = chebfi%bandpp
590 9024 : ndeg_filter = chebfi%ndeg_filter
591 9024 : chebfi%eigenvalues = eigen
592 :
593 9024 : if (chebfi%space==SPACE_C) then
594 5488 : space_res = SPACE_C
595 3536 : else if (chebfi%space==SPACE_CR) then
596 3536 : space_res = SPACE_R
597 : else
598 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
599 : end if
600 :
601 9024 : if (chebfi%paral_kgb == 0) then
602 8160 : ABI_MALLOC(ndeg_filter_bands,(neigenpairs))
603 2720 : call xg_init(DivResults, space_res, neigenpairs, 1, gpu_option=chebfi%gpu_option)
604 : else
605 18912 : ABI_MALLOC(ndeg_filter_bands,(bandpp))
606 6304 : call xg_init(DivResults, space_res, bandpp, 1, gpu_option=chebfi%gpu_option)
607 : end if
608 :
609 9024 : tolerance = chebfi%tolerance
610 9024 : lambda_plus = chebfi%ecut
611 9024 : chebfi%X = X0
612 :
613 : ! Transpose
614 9024 : if (chebfi%paral_kgb == 1) then
615 :
616 6304 : call timab(tim_transpose,1,tsec)
617 : call xgTransposer_constructor(chebfi%xgTransposerX,chebfi%X,chebfi%xXColsRows,nspinor,&
618 : STATE_LINALG,TRANS_ALL2ALL,chebfi%comm_rows,chebfi%comm_cols,0,0,chebfi%me_g0_fft,&
619 6304 : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit)
620 :
621 6304 : call xgTransposer_copyConstructor(chebfi%xgTransposerAX,chebfi%xgTransposerX,chebfi%AX%self,chebfi%xAXColsRows,STATE_LINALG)
622 6304 : call xgTransposer_copyConstructor(chebfi%xgTransposerBX,chebfi%xgTransposerX,chebfi%BX%self,chebfi%xBXColsRows,STATE_LINALG)
623 :
624 6304 : chebfi%xgTransposerX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
625 6304 : chebfi%xgTransposerAX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
626 6304 : chebfi%xgTransposerBX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
627 :
628 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
629 6304 : call xgTransposer_transpose(chebfi%xgTransposerX,STATE_COLSROWS)
630 6304 : chebfi%xgTransposerAX%state = STATE_COLSROWS
631 6304 : chebfi%xgTransposerBX%state = STATE_COLSROWS
632 : ABI_NVTX_END_RANGE()
633 6304 : call timab(tim_transpose,2,tsec)
634 : else
635 2720 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs) !use xXColsRows instead of X notion
636 2720 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs) !use xAXColsRows instead of AX notion
637 2720 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
638 : end if
639 :
640 9024 : call timab(tim_getAX_BX,1,tsec)
641 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
642 9024 : call getAX_BX(chebfi%xXColsRows,chebfi%xAXColsRows,chebfi%xBXColsRows)
643 9024 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
644 9024 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
645 : ABI_NVTX_END_RANGE()
646 9024 : call timab(tim_getAX_BX,2,tsec)
647 :
648 9024 : if (chebfi%paral_kgb == 1) then
649 6304 : call timab(tim_barrier,1,tsec)
650 6304 : call xmpi_barrier(chebfi%spacecom)
651 6304 : call timab(tim_barrier,2,tsec)
652 : end if
653 :
654 : !********************* Compute Rayleigh quotients for every band, and set lambda equal to the largest one *****
655 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RRQ)
656 :
657 : ! NOTICE : the following lines are kept for reference
658 : ! they are no longer necessary as chebfi_rayleighRitzQuotients can fully run on GPU
659 : ! it is no longer necessary to issue a data prefetch, data are already present on device
660 :
661 : ! if (chebfi%gpu_option == ABI_GPU_KOKKOS) then
662 : ! #if defined(HAVE_GPU_CUDA)
663 : ! call xgBlock_prefetch_async(chebfi%xXColsRows, CPU_DEVICE_ID)
664 : ! call xgBlock_prefetch_async(chebfi%xAXColsRows, CPU_DEVICE_ID)
665 : ! call xgBlock_prefetch_async(chebfi%xBXColsRows, CPU_DEVICE_ID)
666 : ! #endif
667 : ! end if
668 :
669 9024 : call timab(tim_RR_q, 1, tsec)
670 9024 : call chebfi_rayleighRitzQuotients(chebfi, maxeig, mineig, DivResults%self)
671 :
672 9024 : if (chebfi%paral_kgb == 1) then
673 6304 : call xmpi_max(maxeig,maxeig_global,chebfi%spacecom,ierr)
674 6304 : call xmpi_min(mineig,mineig_global,chebfi%spacecom,ierr)
675 : else
676 2720 : maxeig_global = maxeig
677 2720 : mineig_global = mineig
678 : end if
679 9024 : call timab(tim_RR_q, 2, tsec)
680 : ABI_NVTX_END_RANGE()
681 :
682 9024 : lambda_minus = maxeig_global
683 :
684 9024 : call timab(tim_oracle,1,tsec)
685 :
686 : ! ndeg_filter_max limits the reduction of the residual of the smallest eigenvalue (i.e. the most amplified one by the filter) by a factor 1e8.
687 : ! Also, the maximal value of ndeg_filter_max is 40.
688 9024 : ndeg_filter_max = cheb_oracle1(mineig_global, lambda_minus, lambda_plus, 1D-16, 40)
689 9024 : ndeg_filter = MIN(ndeg_filter_max,chebfi%ndeg_filter)
690 9024 : if (chebfi%oracle>0) then
691 2312 : call chebfi_set_ndeg_from_residu(chebfi,lambda_minus,lambda_plus,occ,DivResults%self,ndeg_filter_max,ndeg_filter)
692 : end if
693 53520 : ndeg_filter_bands(:) = ndeg_filter
694 :
695 9024 : call timab(tim_oracle,2,tsec)
696 :
697 9024 : center = (lambda_plus + lambda_minus)*0.5
698 9024 : radius = (lambda_plus - lambda_minus)*0.5
699 :
700 9024 : one_over_r = 1/radius
701 9024 : two_over_r = 2/radius
702 :
703 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
704 72586 : do ideg = 0, ndeg_filter - 1
705 :
706 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
707 63562 : call chebfi_computeNextOrderChebfiPolynom(chebfi, ideg, center, one_over_r, two_over_r, getBm1X)
708 : ABI_NVTX_END_RANGE()
709 :
710 63562 : call timab(tim_swap,1,tsec)
711 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
712 63562 : if (chebfi%paral_kgb == 0) then
713 19050 : call chebfi_swapInnerBuffers(chebfi, spacedim, neigenpairs)
714 : else
715 44512 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, bandpp)
716 : end if
717 : ABI_NVTX_END_RANGE()
718 63562 : call timab(tim_swap,2,tsec)
719 :
720 : !A * Psi
721 63562 : call timab(tim_getAX_BX,1,tsec)
722 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
723 63562 : call getAX_BX(chebfi%xXColsRows,chebfi%xAXColsRows,chebfi%xBXColsRows)
724 63562 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
725 63562 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
726 : ABI_NVTX_END_RANGE()
727 72586 : call timab(tim_getAX_BX,2,tsec)
728 :
729 : end do ! ideg
730 : ABI_NVTX_END_RANGE()
731 :
732 9024 : if (chebfi%paral_kgb == 1) then
733 6304 : call timab(tim_barrier,1,tsec)
734 6304 : call xmpi_barrier(chebfi%spacecom)
735 6304 : call timab(tim_barrier,2,tsec)
736 : end if
737 :
738 9024 : call timab(tim_amp_f,1,tsec)
739 9024 : call chebfi_ampfactor(chebfi, DivResults%self, lambda_minus, lambda_plus, ndeg_filter_bands)
740 9024 : call timab(tim_amp_f,2,tsec)
741 :
742 9024 : call xg_free(DivResults)
743 9024 : ABI_SFREE(ndeg_filter_bands)
744 :
745 9024 : call timab(tim_transpose,1,tsec)
746 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
747 9024 : if (chebfi%paral_kgb == 1) then
748 6304 : call xmpi_barrier(chebfi%spacecom)
749 :
750 6304 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_LINALG)
751 6304 : call xgTransposer_transpose(chebfi%xgTransposerAX,STATE_LINALG)
752 6304 : call xgTransposer_transpose(chebfi%xgTransposerBX,STATE_LINALG)
753 :
754 : !only one MPI proc reset buffers to right addresses (because of X-Xcolwise swaps)
755 6304 : if (xmpi_comm_size(chebfi%spacecom) == 1) then
756 0 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
757 0 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
758 0 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
759 : end if
760 : else
761 2720 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
762 2720 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
763 2720 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
764 : end if
765 : ABI_NVTX_END_RANGE()
766 9024 : call timab(tim_transpose,2,tsec)
767 :
768 : ! Apply Rayleigh-Ritz
769 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RR)
770 : call xg_RayleighRitz(chebfi%X,chebfi%AX%self,chebfi%BX%self,chebfi%eigenvalues,ierr,0,tim_RR,&
771 9024 : & chebfi%gpu_option,solve_ax_bx=.true.)
772 : ABI_NVTX_END_RANGE()
773 :
774 : ! Compute residual
775 9024 : call timab(tim_residu, 1, tsec)
776 9024 : if (chebfi%paw) then
777 7528 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%BX%self,chebfi%AX%self)
778 : else
779 1496 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%X,chebfi%AX%self)
780 : end if
781 :
782 9024 : call xgBlock_colwiseNorm2(chebfi%AX%self, residu)
783 9024 : call timab(tim_residu, 2, tsec)
784 :
785 9024 : call timab(tim_copy, 1, tsec)
786 9024 : call xgBlock_copy(chebfi%X,X0)
787 9024 : call timab(tim_copy, 2, tsec)
788 :
789 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
790 : if (chebfi%gpu_option==ABI_GPU_KOKKOS) then
791 : call gpu_device_synchronize()
792 : end if
793 : #endif
794 :
795 9024 : if (chebfi%paral_kgb == 1) then
796 6304 : call xgTransposer_free(chebfi%xgTransposerX)
797 6304 : call xgTransposer_free(chebfi%xgTransposerAX)
798 6304 : call xgTransposer_free(chebfi%xgTransposerBX)
799 : end if
800 :
801 : ! call timab(tim_run,2,tsec)
802 :
803 18048 : end subroutine chebfi_run
804 : !!***
805 :
806 : !----------------------------------------------------------------------
807 :
808 : !!****f* m_chebfi2/chebfi_rayleighRitzQuotients
809 : !! NAME
810 : !! chebfi_rayleighRitzQuotients
811 : !!
812 : !! FUNCTION
813 : !! Compute the Rayleigh-Ritz quotients.
814 : !!
815 : !! INPUTS
816 : !!
817 : !! OUTPUT
818 : !!
819 : !! SIDE EFFECTS
820 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
821 : !! maxeig= highest eigenvalue
822 : !! mineig= lowest eigenvalue
823 : !! DivResults= Rayleigh-Ritz quotients
824 : !!
825 : !! SOURCE
826 :
827 9024 : subroutine chebfi_rayleighRitzQuotients(chebfi,maxeig,mineig,DivResults)
828 :
829 : !Arguments ------------------------------------
830 : real(dp), intent(inout) :: maxeig
831 : real(dp), intent(inout) :: mineig
832 : type(chebfi_t), intent(inout) :: chebfi
833 : type(xgBlock_t), intent(inout) :: DivResults
834 :
835 : !Local variables-------------------------------
836 : !scalars
837 : type(xg_t)::Results1
838 : type(xg_t)::Results2
839 : !arrays
840 : integer :: maxeig_pos(2)
841 : integer :: mineig_pos(2)
842 : integer :: space_res
843 :
844 : ! *********************************************************************
845 :
846 9024 : if (chebfi%space==SPACE_C) then
847 5488 : space_res = SPACE_C
848 3536 : else if (chebfi%space==SPACE_CR) then
849 3536 : space_res = SPACE_R
850 : else
851 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
852 : end if
853 :
854 : !Doesnt work with npfft (ncols=1 in the formula below) ???
855 9024 : if (chebfi%paral_kgb == 0) then
856 2720 : call xg_init(Results1, space_res, chebfi%neigenpairs, 1, gpu_option=chebfi%gpu_option)
857 2720 : call xg_init(Results2, space_res, chebfi%neigenpairs, 1, gpu_option=chebfi%gpu_option)
858 : else
859 6304 : call xg_init(Results1, space_res, chebfi%bandpp, 1, gpu_option=chebfi%gpu_option)
860 6304 : call xg_init(Results2, space_res, chebfi%bandpp, 1, gpu_option=chebfi%gpu_option)
861 : end if
862 :
863 : ! <Psi|H|Psi>
864 9024 : call xgBlock_colwiseDotProduct(chebfi%xXColsRows, chebfi%xAXColsRows, Results1%self, comm_loc=xmpi_comm_null)
865 :
866 : ! <Psi|S|Psi>
867 9024 : call xgBlock_colwiseDotProduct(chebfi%xXColsRows, chebfi%xBXColsRows, Results2%self, comm_loc=xmpi_comm_null)
868 :
869 : ! eig = <Psi|H|Psi> / <Psi|S|Psi>
870 : call xgBlock_colwiseDivision(Results1%self, Results2%self, DivResults, &
871 9024 : & maxeig, maxeig_pos, mineig, mineig_pos)
872 :
873 9024 : call xg_free(Results1)
874 9024 : call xg_free(Results2)
875 :
876 9024 : end subroutine chebfi_rayleighRitzQuotients
877 : !!***
878 :
879 : !----------------------------------------------------------------------
880 :
881 : !!****f* m_chebfi2/chebfi_computeNextOrderChebfiPolynom
882 : !! NAME
883 : !! chebfi_computeNextOrderChebfiPolynom
884 : !!
885 : !! FUNCTION
886 : !! From P_n(B-^1.A)|X> (where P_n is the Chebyshev polynom of order n),
887 : !! computes P_n+1(B-^1.A)|X>
888 : !!
889 : !! INPUTS
890 : !! ideg=current degree of polynom
891 : !! center=filter center
892 : !! one_over_r,two_over_r=1/R, 2/R, R being the radius of the filter
893 : !! getBm1X= pointer to the function giving B^-1|X>
894 : !! B is typically the overlap operator S
895 : !!
896 : !! OUTPUT
897 : !!
898 : !! SIDE EFFECTS
899 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
900 : !!
901 : !! SOURCE
902 :
903 63742 : subroutine chebfi_computeNextOrderChebfiPolynom(chebfi,ideg,center,one_over_r,two_over_r,getBm1X)
904 :
905 : !Arguments ------------------------------------
906 : real(dp) , intent(in) :: center
907 : integer , intent(in) :: ideg
908 : real(dp) , intent(in) :: one_over_r
909 : real(dp) , intent(in) :: two_over_r
910 : type(chebfi_t) , intent(inout) :: chebfi
911 : interface
912 : subroutine getBm1X(X,Bm1X)
913 : use m_xg, only : xgBlock_t
914 : type(xgBlock_t), intent(inout) :: X
915 : type(xgBlock_t), intent(inout) :: Bm1X
916 : end subroutine getBm1X
917 : end interface
918 :
919 : !Local variables-------------------------------
920 : real(dp) :: tsec(2)
921 :
922 : ! *********************************************************************
923 :
924 63742 : if (chebfi%paw) then
925 51898 : call timab(tim_invovl, 1, tsec)
926 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_BM1X)
927 51898 : call getBm1X(chebfi%xAXColsRows, chebfi%X_next)
928 : ABI_NVTX_END_RANGE()
929 51898 : call timab(tim_invovl, 2, tsec)
930 : else
931 11844 : call timab(tim_copy, 1, tsec)
932 11844 : call xgBlock_copy(chebfi%xAXColsRows,chebfi%X_next)
933 11844 : call timab(tim_copy, 2, tsec)
934 : end if
935 :
936 63742 : call timab(tim_postinvovl, 1, tsec)
937 : ABI_NVTX_START_RANGE(NVTX_INVOVL_POST3)
938 63742 : call xgBlock_scale(chebfi%xXColsRows, center, 1) !scale by center
939 :
940 : !(B-1 * A * Psi^i-1 - c * Psi^i-1)
941 63742 : call xgBlock_saxpy(chebfi%X_next, dble(-1.0), chebfi%xXColsRows)
942 :
943 : !Psi^i-1 = 1/c * Psi^i-1
944 63742 : call xgBlock_scale(chebfi%xXColsRows, dble(1.0)/center, 1) !counter scale by 1/center
945 :
946 63742 : if (ideg == 0) then
947 8658 : call xgBlock_scale(chebfi%X_next, one_over_r, 1)
948 : else
949 55084 : call xgBlock_scale(chebfi%X_next, two_over_r, 1)
950 :
951 55084 : call xgBlock_saxpy(chebfi%X_next, dble(-1.0), chebfi%X_prev)
952 : end if
953 :
954 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
955 : if (chebfi%gpu_option==ABI_GPU_KOKKOS) then
956 : call gpu_device_synchronize()
957 : end if
958 : #endif
959 : ABI_NVTX_END_RANGE()
960 63742 : call timab(tim_postinvovl, 2, tsec)
961 :
962 63742 : end subroutine chebfi_computeNextOrderChebfiPolynom
963 : !!***
964 :
965 : !----------------------------------------------------------------------
966 :
967 : !!****f* m_chebfi2/chebfi_swapInnerBuffers
968 : !! NAME
969 : !! chebfi_swapInnerBuffers
970 : !!
971 : !! FUNCTION
972 : !! Swap buffers inside a 'chebfi' datastructure.
973 : !!
974 : !! INPUTS
975 : !! neigenpairs= number of requested eigenvectors/eigenvalues
976 : !! spacedim= space dimension for one vector
977 : !!
978 : !! OUTPUT
979 : !!
980 : !! SIDE EFFECTS
981 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
982 : !!
983 : !! SOURCE
984 :
985 63742 : subroutine chebfi_swapInnerBuffers(chebfi,spacedim,neigenpairs)
986 :
987 : ! Arguments ------------------------------------
988 : integer , intent(in ) :: spacedim
989 : integer , intent(in ) :: neigenpairs
990 : type(chebfi_t) , intent(inout) :: chebfi
991 :
992 : ! *********************************************************************
993 :
994 63742 : call xgBlock_setBlock(chebfi%X_prev, chebfi%X_swap, spacedim, neigenpairs) !X_swap = X_prev
995 63742 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X_prev, spacedim, neigenpairs) !X_prev = xXColsRows
996 63742 : call xgBlock_setBlock(chebfi%X_next, chebfi%xXColsRows, spacedim, neigenpairs) !xXColsRows = X_next
997 63742 : call xgBlock_setBlock(chebfi%X_swap, chebfi%X_next, spacedim, neigenpairs) !X_next = X_swap
998 :
999 63742 : end subroutine chebfi_swapInnerBuffers
1000 : !!***
1001 :
1002 : !----------------------------------------------------------------------
1003 :
1004 : !!****f* m_chebfi2/chebfi_prepareAmpfactor
1005 : !! NAME
1006 : !! chebfi_prepAmpfactor
1007 : !!
1008 : !! FUNCTION
1009 : !! Utility function to distribute correctly eigenvalues to MPI procs.
1010 : !! Performs MPI communication
1011 : !!
1012 : !! INPUTS
1013 : !! chebfi=
1014 : !! eigen= eigenvalues not distributed
1015 : !!
1016 : !! OUTPUT
1017 : !! DivResults= prepared eigenvalues same array as in chebfi_run
1018 : !!
1019 : !! SOURCE
1020 :
1021 0 : subroutine chebfi_prepAmpfactor(chebfi, eigen, DivResults)
1022 :
1023 : implicit none
1024 :
1025 : ! Arguments ------------------------------------
1026 : type(xg_t), intent(inout) :: DivResults
1027 : type(xgBlock_t), intent(inout) :: eigen
1028 : type(chebfi_t), intent(inout) :: chebfi
1029 :
1030 : ! Local variables-------------------------------
1031 : ! scalars
1032 : integer :: space_res
1033 : integer :: my_rank, num_proc, shift, ierr
1034 : type(xgBlock_t) :: eigen_block
1035 : ! Arrays
1036 0 : integer, allocatable, target :: allbandpp(:)
1037 0 : real(dp), allocatable, target :: theta_reshaped(:,:)
1038 : integer, pointer :: allbandpp_ptr(:) => null()
1039 : real(dp), pointer :: theta_reshaped_ptr(:,:) => null()
1040 : real(dp), pointer :: theta(:,:) => null()
1041 :
1042 : ! *********************************************************************
1043 :
1044 0 : if (chebfi%space==SPACE_C) then
1045 0 : space_res = SPACE_C
1046 0 : else if (chebfi%space==SPACE_CR) then
1047 0 : space_res = SPACE_R
1048 : else
1049 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
1050 : end if
1051 :
1052 0 : if (chebfi%paral_kgb == 0) then
1053 0 : call xg_init(DivResults, space_res, rows=chebfi%neigenpairs, cols=1, gpu_option=chebfi%gpu_option)
1054 : ! Fill DivResults with full eigenvalues
1055 : ! TODO fix for workaround copy between space_res and SPACE_R
1056 0 : if (space_res==SPACE_R) then
1057 0 : call xgBlock_copy(eigen, DivResults%self)
1058 : else
1059 : ! workaround to copy from SPACE_R to SPACE_C
1060 0 : ABI_MALLOC_IFNOT(theta_reshaped,(2,chebfi%neigenpairs))
1061 0 : theta_reshaped_ptr => theta_reshaped
1062 0 : call xgBlock_reverseMap(eigen, theta, rows=1, cols=chebfi%neigenpairs)
1063 0 : theta_reshaped = 0.d0
1064 0 : theta_reshaped(1,1:chebfi%neigenpairs) = theta(1,1:chebfi%neigenpairs)
1065 : #ifdef HAVE_OPENMP_OFFLOAD
1066 : !$OMP TARGET ENTER DATA MAP(to:theta_reshaped) IF(chebfi%gpu_option==ABI_GPU_OPENMP)
1067 : #endif
1068 : call xgBlock_map(eigen, theta_reshaped_ptr, space_res, rows=1, &
1069 0 : cols=chebfi%neigenpairs, gpu_option=chebfi%gpu_option)
1070 0 : call xgBlock_copy(eigen, DivResults%self)
1071 : #ifdef HAVE_OPENMP_OFFLOAD
1072 : !$OMP TARGET EXIT DATA MAP(delete:theta_reshaped) IF(chebfi%gpu_option==ABI_GPU_OPENMP)
1073 : #endif
1074 0 : ABI_SFREE(theta_reshaped)
1075 : end if
1076 : else
1077 0 : call xg_init(DivResults, space_res, chebfi%bandpp, 1, gpu_option=chebfi%gpu_option)
1078 0 : if (xmpi_comm_size(chebfi%spacecom) > 1) then
1079 0 : my_rank = xmpi_comm_rank(chebfi%spacecom)
1080 : !shift = my_rank * chebfi%bandpp ! FIXME not working for different bandpp per rank
1081 0 : num_proc = xmpi_comm_size(chebfi%spacecom)
1082 0 : ABI_MALLOC_IFNOT(allbandpp,(num_proc))
1083 0 : allbandpp_ptr => allbandpp
1084 0 : call xmpi_allgather(chebfi%bandpp, allbandpp_ptr, chebfi%spacecom, ierr)
1085 0 : if ( ierr /= xmpi_success ) then
1086 0 : ABI_ERROR("Error while gathering number of bandpp for spacecom")
1087 : end if
1088 0 : if (my_rank==0) then
1089 : shift = 0
1090 : else
1091 0 : shift = sum(allbandpp(1:my_rank)) ! fixed
1092 : end if
1093 0 : ABI_SFREE(allbandpp)
1094 : else
1095 : shift = 0
1096 : end if
1097 : ! Fill DivResults(bandpp,1) with block of eigen(neigenpairs,1) of size bandpp
1098 : ! reshape to access column range
1099 0 : call xgBlock_reshape(DivResults%self, 1, chebfi%bandpp)
1100 0 : call xgBlock_reshape(eigen, 1, chebfi%neigenpairs)
1101 0 : if (space_res==SPACE_R) then
1102 0 : call xgBlock_setBlock(eigen, eigen_block, rows=1, cols=chebfi%bandpp, fcol=1+shift)
1103 0 : call xgBlock_copy(eigen_block, DivResults%self)
1104 : else
1105 : ! workaround to copy from SPACE_R to SPACE_C
1106 0 : ABI_MALLOC_IFNOT(theta_reshaped,(2,chebfi%bandpp))
1107 0 : theta_reshaped_ptr => theta_reshaped
1108 0 : call xgBlock_setBlock(eigen, eigen_block, rows=1, cols=chebfi%bandpp, fcol=1+shift)
1109 0 : call xgBlock_reverseMap(eigen_block, theta, rows=1, cols=chebfi%bandpp)
1110 0 : theta_reshaped = 0.d0
1111 0 : theta_reshaped(1,1:chebfi%bandpp) = theta(1,1:chebfi%bandpp)
1112 : #ifdef HAVE_OPENMP_OFFLOAD
1113 : !$OMP TARGET ENTER DATA MAP(to:theta_reshaped) IF(chebfi%gpu_option==ABI_GPU_OPENMP)
1114 : #endif
1115 : call xgBlock_map(eigen_block, theta_reshaped_ptr, space_res, rows=1, &
1116 0 : cols=chebfi%bandpp, gpu_option=chebfi%gpu_option)
1117 0 : call xgBlock_copy(eigen_block, DivResults%self)
1118 : #ifdef HAVE_OPENMP_OFFLOAD
1119 : !$OMP TARGET EXIT DATA MAP(delete:theta_reshaped) IF(chebfi%gpu_option==ABI_GPU_OPENMP)
1120 : #endif
1121 0 : ABI_SFREE(theta_reshaped)
1122 : end if
1123 : ! restore dimensions
1124 0 : call xgBlock_reshape(eigen, chebfi%neigenpairs, 1)
1125 0 : call xgBlock_reshape(DivResults%self, chebfi%bandpp, 1)
1126 : end if
1127 :
1128 : ! DivResults must be on CPU for ampfactor routine
1129 0 : if (chebfi%gpu_option==ABI_GPU_OPENMP) then
1130 0 : call xgBlock_copy_from_gpu(DivResults%self)
1131 : end if
1132 :
1133 0 : end subroutine chebfi_prepAmpfactor
1134 : !!***
1135 :
1136 : !----------------------------------------------------------------------
1137 :
1138 : !!****f* m_chebfi2/chebfi_ampfactor
1139 : !! NAME
1140 : !! chebfi_ampfactor
1141 : !!
1142 : !! FUNCTION
1143 : !! Compute amplification factor
1144 : !!
1145 : !! INPUTS
1146 : !! eig (:,:)= eigenvalues
1147 : !! lambda_minus,lambda_plus=
1148 : !! ndeg_filter_bands(:)= degree of Chebyshev polynomial filter for each band
1149 : !!
1150 : !! OUTPUT
1151 : !!
1152 : !! SIDE EFFECTS
1153 : !! residu<type(xgBlock_t)>= vector of residuals
1154 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
1155 : !!
1156 : !! SOURCE
1157 :
1158 9024 : subroutine chebfi_ampfactor(chebfi,DivResults,lambda_minus,lambda_plus,ndeg_filter_bands)
1159 :
1160 : ! Arguments ------------------------------------
1161 : integer, intent(in ) :: ndeg_filter_bands(:)
1162 : type(xgBlock_t), intent(in ) :: DivResults
1163 : real(dp), intent(in ) :: lambda_minus
1164 : real(dp), intent(in ) :: lambda_plus
1165 : type(chebfi_t), intent(inout) :: chebfi
1166 :
1167 : ! Local variables-------------------------------
1168 : ! scalars
1169 : integer :: iband,nbands
1170 : real(dp) :: ampfactor
1171 : real(dp) :: eig_per_band
1172 : type(xgBlock_t) :: X_part
1173 : type(xgBlock_t) :: AX_part
1174 : type(xgBlock_t) :: BX_part
1175 9024 : real(dp),pointer :: eig(:,:)
1176 :
1177 : ! *********************************************************************
1178 :
1179 9024 : if (chebfi%paral_kgb == 0) then
1180 2720 : nbands = chebfi%neigenpairs
1181 : else
1182 6304 : nbands = chebfi%bandpp
1183 : end if
1184 :
1185 9024 : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=chebfi%bandpp)
1186 :
1187 53520 : do iband = 1, nbands
1188 :
1189 44496 : eig_per_band = eig(1,iband)
1190 :
1191 : !cheb_poly1(x, n, a, b)
1192 44496 : ampfactor = cheb_poly1(eig_per_band, ndeg_filter_bands(iband), lambda_minus, lambda_plus)
1193 :
1194 44496 : if(abs(ampfactor) < 1e-3) ampfactor = 1e-3 !just in case, avoid amplifying too much
1195 :
1196 44496 : call xgBlock_setBlock(chebfi%xXColsRows, X_part, chebfi%total_spacedim, 1, fcol=iband)
1197 44496 : call xgBlock_setBlock(chebfi%xAXColsRows, AX_part, chebfi%total_spacedim, 1, fcol=iband)
1198 44496 : call xgBlock_setBlock(chebfi%xBXColsRows, BX_part, chebfi%total_spacedim, 1, fcol=iband)
1199 :
1200 44496 : call xgBlock_scale(X_part, 1/ampfactor, 1)
1201 44496 : call xgBlock_scale(AX_part, 1/ampfactor, 1)
1202 53520 : call xgBlock_scale(BX_part, 1/ampfactor, 1)
1203 :
1204 : end do
1205 :
1206 9024 : end subroutine chebfi_ampfactor
1207 : !!***
1208 :
1209 : !----------------------------------------------------------------------
1210 :
1211 : !!****f* m_chebfi2/chebfi_ampfactorBandpass
1212 : !! NAME
1213 : !! chebfi_ampfactorBandpass
1214 : !!
1215 : !! FUNCTION
1216 : !! Compute amplification factor for bandpass polynomial
1217 : !! Numerical zero is 1e-3. Assumes prepAmpfactor prior to this.
1218 : !!
1219 : !! INPUTS
1220 : !! eig (:,:)= eigenvalues
1221 : !! lambda_minus,lambda_plus=
1222 : !! center,radius= used to rescale bandpass filter to interval
1223 : !! ndeg_filter= degree of bandpass polynomial filter
1224 : !!
1225 : !! OUTPUT
1226 : !!
1227 : !! SIDE EFFECTS
1228 : !! residu<type(xgBlock_t)>= vector of residuals
1229 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
1230 : !!
1231 : !! SOURCE
1232 :
1233 : subroutine chebfi_ampfactorBandpass(chebfi,DivResults,lambda_minus,lambda_plus,center,radius,ndeg_filter)
1234 :
1235 : implicit none
1236 :
1237 : ! Arguments ------------------------------------
1238 : integer, intent(in ) :: ndeg_filter
1239 : type(xgBlock_t), intent(in ) :: DivResults
1240 : real(dp), intent(in ) :: lambda_minus
1241 : real(dp), intent(in ) :: lambda_plus
1242 : real(dp), intent(in ) :: center
1243 : real(dp), intent(in ) :: radius
1244 : type(chebfi_t), intent(inout) :: chebfi
1245 :
1246 : ! Local variables-------------------------------
1247 : ! scalars
1248 : integer :: iband,nbands
1249 : real(dp) :: ampfactor
1250 : real(dp) :: xred, ls, us
1251 : real(dp) :: eig_per_band
1252 : type(xgBlock_t) :: X_part
1253 : type(xgBlock_t) :: AX_part
1254 : type(xgBlock_t) :: BX_part
1255 : real(dp),pointer :: eig(:,:)
1256 :
1257 : ! *********************************************************************
1258 :
1259 : if (chebfi%paral_kgb == 0) then
1260 : nbands = chebfi%neigenpairs
1261 : else
1262 : nbands = chebfi%bandpp
1263 : end if
1264 : ls = (lambda_minus-center)/radius
1265 : us = (lambda_plus-center)/radius
1266 :
1267 : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=chebfi%bandpp)
1268 :
1269 : do iband = 1, nbands
1270 :
1271 : eig_per_band = eig(1,iband)
1272 :
1273 : !poly(x, a, b, n), where x,a,b are scaled!!!
1274 : xred = (eig_per_band-center)/radius
1275 : ampfactor = bandpassIndicator_sca(xred, ls, us, ndeg_filter)
1276 :
1277 : if(abs(ampfactor) < 1e-3) ampfactor = 1e-3 !just in case, avoid amplifying too much
1278 :
1279 : call xgBlock_setBlock(chebfi%xXColsRows, X_part, chebfi%total_spacedim, 1, fcol=iband)
1280 : call xgBlock_setBlock(chebfi%xAXColsRows, AX_part, chebfi%total_spacedim, 1, fcol=iband)
1281 : call xgBlock_setBlock(chebfi%xBXColsRows, BX_part, chebfi%total_spacedim, 1, fcol=iband)
1282 :
1283 : !write(std_out,*) 'ampfactor, eig, iband=', eig_per_band, ampfactor, iband
1284 :
1285 : call xgBlock_scale(X_part, 1/ampfactor, 1)
1286 : call xgBlock_scale(AX_part, 1/ampfactor, 1)
1287 : call xgBlock_scale(BX_part, 1/ampfactor, 1)
1288 :
1289 : end do
1290 :
1291 : end subroutine chebfi_ampfactorBandpass
1292 : !!***
1293 :
1294 : !----------------------------------------------------------------------
1295 :
1296 : !!****f* m_chebfi/chebfi_runSlice
1297 : !! NAME
1298 : !! chebfi_runSlice
1299 : !!
1300 : !! FUNCTION
1301 : !! Apply Polynomial filtering to set of vectors.
1302 : !!
1303 : !! NOTES
1304 : !! List of differences with chebfi_run:
1305 : !! - X0 input and output is in ColsRows MPI representation (not Linalg!).
1306 : !! - the filter polynomial can be lowpass or bandpass.
1307 : !! - Rayleigh Quotients are read from eigen.
1308 : !!
1309 : !! INPUT
1310 : !! chebfi <type(chebfi_t)>= all data used to apply Polynomial Filtering algorithm
1311 : !! X0= eigenvector guess distributed in ColsRows representation
1312 : !! getAX_BX= pointer to the function giving A|X> and B|X>
1313 : !! A is typically the Hamiltonian H, and B the overlap operator S
1314 : !! getBm1X= pointer to the function giving B^-1|X>
1315 : !! B is typically the overlap operator S
1316 : !! eigen= Rayleigh quotients associated to X0
1317 : !! residu= empty array
1318 : !! nspinor= number of spinorial components of the wavefunctions
1319 : !! lambda_minus= lower interval to amplify/diminish for bandpass/lowpass
1320 : !! lambda_plus= upper interval to amplify/diminish for bandpass/lowpass
1321 : !! mineig_global= guaranteed lower bound for entire spectrum
1322 : !! maxeig_global= guaranteed upper bound for entire spectrum
1323 : !! is_lowpass= flag. True if Chebyshev otherwise use bandpass Chebyshev-Jackson
1324 : !! nrows_blockrows= number of rows per MPI block in Linalg representation
1325 : !!
1326 : !! SIDE EFFECTS
1327 : !! chebfi= workspaces used
1328 : !! X0= full set of vectors distributed in ColsRows representation
1329 : !! eigen= full eigenvalues
1330 : !! residu= residuals, i.e. norm of (A-lambdaB)|X>
1331 : !!
1332 : !! SOURCE
1333 :
1334 0 : subroutine chebfi_runSlice(chebfi,X0,getAX_BX,getBm1X,eigen,residu,nspinor,&
1335 : mineig_global,maxeig_global,lambda_minus,lambda_plus,is_lowpass,nrows_blockrows)
1336 :
1337 : implicit none
1338 :
1339 : !Arguments ------------------------------------
1340 : type(chebfi_t) , intent(inout) :: chebfi
1341 : type(xgBlock_t), intent(inout) :: X0
1342 : type(xgBlock_t), intent(inout) :: eigen
1343 : type(xgBlock_t), intent(inout) :: residu
1344 : integer , intent(in ) :: nspinor
1345 : integer, pointer, intent(in ) :: nrows_blockrows(:)
1346 : real(dp) , intent(in ) :: mineig_global
1347 : real(dp) , intent(in ) :: maxeig_global
1348 : real(dp) , intent(in ) :: lambda_minus
1349 : real(dp) , intent(in ) :: lambda_plus
1350 : logical , intent(in ) :: is_lowpass
1351 : interface
1352 : subroutine getAX_BX(X,AX,BX)
1353 : use m_xg, only : xgBlock_t
1354 : type(xgBlock_t), intent(inout) :: X
1355 : type(xgBlock_t), intent(inout) :: AX
1356 : type(xgBlock_t), intent(inout) :: BX
1357 : end subroutine getAX_BX
1358 : end interface
1359 : interface
1360 : subroutine getBm1X(X,Bm1X)
1361 : use m_xg, only : xgBlock_t
1362 : type(xgBlock_t), intent(inout) :: X
1363 : type(xgBlock_t), intent(inout) :: Bm1X
1364 : end subroutine getBm1X
1365 : end interface
1366 :
1367 : !Local variables-------------------------------
1368 : integer :: spacedim, neigenpairs, num_proc, ierr
1369 : !type(xg_t) :: X_k
1370 : ! Arrays
1371 0 : integer, target, allocatable :: nrowsLinalg(:)
1372 : integer, pointer :: nrowsLinalg_ptr(:) => null()
1373 : real(dp) :: tsec(2)
1374 :
1375 : ! *********************************************************************
1376 :
1377 0 : if (chebfi%from_linalg) then
1378 0 : ABI_ERROR("chebfi should be from colsrows")
1379 : end if
1380 :
1381 0 : spacedim = chebfi%spacedim
1382 0 : neigenpairs = chebfi%neigenpairs
1383 0 : num_proc = xmpi_comm_size(chebfi%spacecom)
1384 0 : chebfi%eigenvalues = eigen
1385 :
1386 0 : ABI_MALLOC_IFNOT(nrowsLinalg,(num_proc))
1387 0 : nrowsLinalg_ptr => nrowsLinalg
1388 0 : nrowsLinalg = nrows_blockrows
1389 :
1390 : !write(std_out,*) 'getid inside runSlice xXColsRows', xgBlock_getId(chebfi%xXColsRows)
1391 : !write(std_out,*) 'getid inside runSlice xAXColsRows', xgBlock_getId(chebfi%xAXColsRows)
1392 : !flush(std_out)
1393 :
1394 : !A * Psi
1395 0 : call timab(tim_getAX_BX,1,tsec)
1396 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
1397 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
1398 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
1399 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
1400 : ABI_NVTX_END_RANGE()
1401 0 : call timab(tim_getAX_BX,2,tsec)
1402 :
1403 : !write(std_out,*) 'getid inside runSlice xXColsRows (filtered 1)', xgBlock_getId(chebfi%xXColsRows)
1404 : !write(std_out,*) 'getid inside runSlice xAXColsRows (filtered 1)', xgBlock_getId(chebfi%xAXColsRows)
1405 :
1406 0 : write(std_out,*) 'starting filter in proc', xmpi_comm_rank(chebfi%spacecom)
1407 0 : flush(std_out)
1408 :
1409 : ! Apply polynomial filtering to active MPI ColsRows block-column
1410 0 : if (is_lowpass) then
1411 : ! [lambda_minus,lambda_plus) is diminished using Chebyshev
1412 0 : write(std_out,*) 'lambda_minus=', lambda_minus
1413 0 : write(std_out,*) 'lambda_plus=', lambda_plus
1414 0 : flush(std_out)
1415 0 : call chebfi_lowpassFilter(chebfi,eigen,lambda_minus,lambda_plus,getAX_BX,getBm1X)
1416 : else
1417 0 : write(std_out,*) 'lambda_minus=', lambda_minus
1418 0 : write(std_out,*) 'lambda_plus=', lambda_plus
1419 0 : write(std_out,*) 'mineig_global=', mineig_global
1420 0 : write(std_out,*) 'maxeig_global=', maxeig_global
1421 0 : flush(std_out)
1422 : call chebfi_bandpassFilter(chebfi,lambda_minus,lambda_plus,mineig_global,&
1423 0 : maxeig_global,getAX_BX,getBm1X)
1424 : end if
1425 :
1426 : !write(std_out,*) 'getid inside runSlice xXColsRows (filtered N)', xgBlock_getId(chebfi%xXColsRows)
1427 : !write(std_out,*) 'getid inside runSlice xAXColsRows (filtered N)', xgBlock_getId(chebfi%xAXColsRows)
1428 : !flush(std_out)
1429 :
1430 : ! MPI transpose to linalg state
1431 0 : call timab(tim_transpose,1,tsec)
1432 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1433 0 : if (chebfi%paral_kgb==1) then
1434 :
1435 : ! Allocate chebfi%X
1436 : call xgTransposer_constructor(chebfi%xgTransposerX,chebfi%X,chebfi%xXColsRows,nspinor,&
1437 : STATE_COLSROWS,TRANS_ALL2ALL,chebfi%comm_rows,chebfi%comm_cols,0,0,chebfi%me_g0,&
1438 : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit,&
1439 0 : custom_ncolsColsRows=.true.,nrowsLinalg_sub=nrowsLinalg_ptr)
1440 : ! Note: bandpp is custom because it is created from resource allocator
1441 :
1442 : ! Allocate chebfi%AX, chebfi%BX
1443 : call xgTransposer_copyConstructor(chebfi%xgTransposerAX,chebfi%xgTransposerX,&
1444 0 : chebfi%AX%self,chebfi%xAXColsRows,STATE_COLSROWS)
1445 : call xgTransposer_copyConstructor(chebfi%xgTransposerBX,chebfi%xgTransposerX,&
1446 0 : chebfi%BX%self,chebfi%xBXColsRows,STATE_COLSROWS)
1447 : ! Note: at this point chebfi%AX and chebfi%BX are empty. Must transpose
1448 : ! to fill with correct values.
1449 :
1450 : ! todo use copy constructor to create an object by copying an existing object
1451 : ! actually copy constructor *allocates* memory for chebfi%AX. Write a version
1452 : ! that does not allocate memory and only reassigns pointers. Can pointers be reassigned
1453 : ! directly then used in the global constructor? Perform tests.
1454 :
1455 : !write(std_out,*) 'getid before transpose AX', xgBlock_getId(chebfi%AX%self)
1456 : !write(std_out,*) 'getid before transpose xAX', xgBlock_getId(chebfi%xAXColsRows)
1457 : !write(std_out,*) 'getid before transpose xX', xgBlock_getId(chebfi%xXColsRows)
1458 : !flush(std_out)
1459 :
1460 0 : chebfi%xgTransposerX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1461 0 : chebfi%xgTransposerAX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1462 0 : chebfi%xgTransposerBX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1463 :
1464 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_LINALG)
1465 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_LINALG)
1466 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_LINALG)
1467 0 : call xmpi_barrier(chebfi%spacecom)
1468 :
1469 : !write(std_out,*) 'getid after transpose AX', xgBlock_getId(chebfi%AX%self)
1470 : !flush(std_out)
1471 :
1472 : else
1473 0 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
1474 0 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
1475 0 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
1476 : end if
1477 0 : call timab(tim_transpose,2,tsec)
1478 : ABI_NVTX_END_RANGE()
1479 :
1480 0 : if (chebfi%paral_kgb==1 .and. rows(chebfi%X) /= nrowsLinalg(xmpi_comm_rank(chebfi%spacecom)+1)) then
1481 0 : ABI_ERROR("wrong linalg representation")
1482 : end if
1483 0 : write(std_out,'(a,i6,i6)') 'local # proc has # rows ', xmpi_comm_rank(chebfi%spacecom), rows(chebfi%X)
1484 :
1485 : !write(std_out,*) 'chebfi%eigenvalues before RR'
1486 : !call xgBlock_print(chebfi%eigenvalues,std_out)
1487 :
1488 : !write(std_out,*) 'id of X, (before RR) ncols=', xgBlock_getId(chebfi%X), cols(chebfi%X)
1489 : !flush(std_out)
1490 :
1491 : !call xg_Borthonormalize(chebfi%xXColsRows,chebfi%xBxColsRows,ierr,1,chebfi%gpu_option,AX=chebfi%xAXColsRows)
1492 :
1493 : ! Apply Rayleigh-Ritz to active MPI Linalg row-block
1494 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RR)
1495 : call xg_RayleighRitz(chebfi%X,chebfi%AX%self,chebfi%BX%self,eigen,ierr,0,tim_RR,&
1496 0 : chebfi%gpu_option,solve_ax_bx=.true.)
1497 : ABI_NVTX_END_RANGE()
1498 :
1499 : !write(std_out,*) 'id of X, (after RR) ncols=', xgBlock_getId(chebfi%X), cols(chebfi%X)
1500 :
1501 0 : if ( ierr /= 0 ) then
1502 0 : ABI_WARNING("RayleighRitz did not work")
1503 : else
1504 : !write(std_out,*) 'is lowpass=', is_lowpass
1505 : !write(std_out,*) 'chebfi%eigenvalues after RR'
1506 : !call xgBlock_print(chebfi%eigenvalues,std_out)
1507 : !flush(std_out)
1508 : end if
1509 :
1510 : ! Compute residual norm *squared*
1511 0 : if (chebfi%paw) then
1512 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%BX%self,chebfi%AX%self)
1513 : else
1514 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%X,chebfi%AX%self)
1515 : end if
1516 0 : call xgBlock_colwiseNorm2(chebfi%AX%self,residu) ! performs MPI comm
1517 :
1518 : !write(std_out,*) 'max colwise residual norm squared='; call xgBlock_print(residu, std_out)
1519 : !flush(std_out)
1520 :
1521 : ! Copy in Linalg representation (see chebfi_run, kept for reference)
1522 : ! call xgBlock_copy(chebfi%X,X0)
1523 :
1524 : ! MPI Transpose to recover colsrows state (X only)
1525 0 : call timab(tim_transpose,1,tsec)
1526 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1527 0 : if (chebfi%paral_kgb == 1) then
1528 0 : call xmpi_barrier(chebfi%spacecom)
1529 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
1530 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
1531 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
1532 0 : if (xmpi_comm_size(chebfi%spacecom) == 1) then
1533 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
1534 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
1535 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
1536 : end if
1537 : else
1538 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
1539 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
1540 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
1541 : end if
1542 : ABI_NVTX_END_RANGE()
1543 0 : call timab(tim_transpose,2,tsec)
1544 :
1545 : ! Copy in ColsRows representation
1546 0 : call xgBlock_copy(chebfi%xXColsRows, X0)
1547 :
1548 0 : if (cols(X0) /= chebfi%bandpp) then
1549 0 : ABI_ERROR('wrong colsrows representation')
1550 : end if
1551 0 : write(std_out,'(a,i6,i6,i6)') 'local # proc has # rows cols ', xmpi_comm_rank(chebfi%spacecom), rows(X0), cols(X0)
1552 :
1553 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
1554 : if (chebfi%gpu_option==ABI_GPU_KOKKOS) then
1555 : call gpu_device_synchronize()
1556 : end if
1557 : #endif
1558 :
1559 : ! Free transposer objects
1560 0 : if (chebfi%paral_kgb == 1) then
1561 0 : call xgTransposer_free(chebfi%xgTransposerX)
1562 0 : call xgTransposer_free(chebfi%xgTransposerAX)
1563 0 : call xgTransposer_free(chebfi%xgTransposerBX)
1564 : end if
1565 :
1566 : ! Free temporary memory
1567 0 : ABI_SFREE(nrowsLinalg)
1568 :
1569 0 : end subroutine chebfi_runSlice
1570 : !!***
1571 :
1572 : !----------------------------------------------------------------------
1573 :
1574 : !!****f* m_chebfi/chebfi_runSubspaceIteration
1575 : !! NAME
1576 : !! chebfi_runSubspaceIteration
1577 : !!
1578 : !! FUNCTION
1579 : !! Applies subspace iteration by Chebyshev polynomial filtering on slice
1580 : !! then extracts eigenvectors using Rayleigh-Ritz. This execution is completely
1581 : !! local on slice processors. Each slice does not see others.
1582 : !!
1583 : !! NOTES
1584 : !! Restart logic in subspace iteration:
1585 : !! Allocate constructor (state=colsrows)
1586 : !! while (convergence not reached):
1587 : !! | if (state=linalg) Transpose
1588 : !! | Filter
1589 : !! | Transpose (state=linalg)
1590 : !! | Orthogonalize
1591 : !! Rayleigh-Ritz
1592 : !!
1593 : !! INPUTS
1594 : !! X0 input vectors distibuted by colsrows along slice processors
1595 : !!
1596 : !! OUTPUTS
1597 : !! X0 (in-place) converged vectors distributed by linalg along slice processors
1598 : !!
1599 : !! dev=============
1600 : !! - version 1: count number of converged vectors within iteration
1601 : !!
1602 : !! SOURCE
1603 :
1604 0 : subroutine chebfi_runSubspaceIteration(chebfi,X0,getAX_BX,getBm1X,eigen,residu,nspinor,&
1605 : mineig_global,maxeig_global,lambda_minus,lambda_plus,is_lowpass,k_conv,nrows_blockrows)
1606 :
1607 : implicit none
1608 :
1609 : !Arguments ------------------------------------
1610 : type(chebfi_t) , intent(inout) :: chebfi
1611 : type(xgBlock_t), intent(inout) :: X0
1612 : type(xgBlock_t), intent(inout) :: eigen
1613 : type(xgBlock_t), intent(inout) :: residu
1614 : integer , intent(in ) :: nspinor
1615 : integer , intent(in ) :: k_conv
1616 : integer, pointer, intent(in ) :: nrows_blockrows(:)
1617 : real(dp) , intent(in ) :: mineig_global
1618 : real(dp) , intent(in ) :: maxeig_global
1619 : real(dp) , intent(in ) :: lambda_minus
1620 : real(dp) , intent(in ) :: lambda_plus
1621 : logical , intent(in ) :: is_lowpass
1622 : interface
1623 : subroutine getAX_BX(X,AX,BX)
1624 : use m_xg, only : xgBlock_t
1625 : type(xgBlock_t), intent(inout) :: X
1626 : type(xgBlock_t), intent(inout) :: AX
1627 : type(xgBlock_t), intent(inout) :: BX
1628 : end subroutine getAX_BX
1629 : end interface
1630 : interface
1631 : subroutine getBm1X(X,Bm1X)
1632 : use m_xg, only : xgBlock_t
1633 : type(xgBlock_t), intent(inout) :: X
1634 : type(xgBlock_t), intent(inout) :: Bm1X
1635 : end subroutine getBm1X
1636 : end interface
1637 :
1638 : !Local variables-------------------------------
1639 : integer :: iter_subspace, niter_subspace_max, n_locked, n_locked_prev
1640 : integer :: spacedim, neigenpairs, num_proc, nrows, nrows_mpi, ierr
1641 : !type(xg_t) :: X_k
1642 : type(xg_t) :: resid_active
1643 : type(bandPartitionInfo_t) :: mpiInfo
1644 0 : type(bandPartitionData_t) :: X_part
1645 0 : type(bandPartitionData_t) :: AX_part
1646 0 : type(bandPartitionData_t) :: BX_part
1647 : ! Arrays
1648 0 : integer, target, allocatable :: nrowsLinalg(:)
1649 : integer, pointer :: nrowsLinalg_ptr(:) => null()
1650 : real(dp) :: tsec(2)
1651 :
1652 : ! *********************************************************************
1653 :
1654 0 : if (chebfi%from_linalg) then
1655 0 : ABI_ERROR("chebfi should be from colsrows")
1656 : end if
1657 :
1658 0 : spacedim = chebfi%spacedim
1659 0 : neigenpairs = chebfi%neigenpairs
1660 0 : num_proc = xmpi_comm_size(chebfi%spacecom)
1661 0 : chebfi%eigenvalues = eigen
1662 :
1663 0 : ABI_MALLOC_IFNOT(nrowsLinalg,(num_proc))
1664 0 : nrowsLinalg_ptr => nrowsLinalg
1665 0 : nrowsLinalg = nrows_blockrows
1666 :
1667 0 : call xg_init(resid_active, SPACE_R, neigenpairs, 1, gpu_option=chebfi%gpu_option)
1668 :
1669 : ! Initialize values of AX and BX
1670 0 : call timab(tim_getAX_BX,1,tsec)
1671 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
1672 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
1673 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
1674 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
1675 : ABI_NVTX_END_RANGE()
1676 0 : call timab(tim_getAX_BX,2,tsec)
1677 :
1678 : ! Allocate memory for linalg distribution from colsrows
1679 : ! Note: bandpp is custom because it is created from resource allocator
1680 0 : call timab(tim_transpose,1,tsec)
1681 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1682 0 : if (chebfi%paral_kgb==1) then
1683 :
1684 : ! Allocate chebfi%X
1685 : call xgTransposer_constructor(chebfi%xgTransposerX,chebfi%X,chebfi%xXColsRows,nspinor,&
1686 : STATE_COLSROWS,TRANS_ALL2ALL,chebfi%comm_rows,chebfi%comm_cols,0,0,chebfi%me_g0,&
1687 : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit,&
1688 0 : custom_ncolsColsRows=.true.,nrowsLinalg_sub=nrowsLinalg_ptr)
1689 :
1690 : ! Allocate chebfi%AX, chebfi%BX
1691 : call xgTransposer_copyConstructor(chebfi%xgTransposerAX,chebfi%xgTransposerX,&
1692 0 : chebfi%AX%self,chebfi%xAXColsRows,STATE_COLSROWS)
1693 : call xgTransposer_copyConstructor(chebfi%xgTransposerBX,chebfi%xgTransposerX,&
1694 0 : chebfi%BX%self,chebfi%xBXColsRows,STATE_COLSROWS)
1695 : ! Note: at this point chebfi%AX and chebfi%BX are empty. Must transpose
1696 : ! to fill with correct values.
1697 :
1698 0 : chebfi%xgTransposerX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1699 0 : chebfi%xgTransposerAX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1700 0 : chebfi%xgTransposerBX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1701 :
1702 : else
1703 0 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
1704 0 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
1705 0 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
1706 : end if
1707 0 : call timab(tim_transpose,2,tsec)
1708 : ABI_NVTX_END_RANGE()
1709 :
1710 0 : write(std_out,*) 'starting filter in proc', xmpi_comm_rank(chebfi%spacecom)
1711 0 : flush(std_out)
1712 :
1713 0 : niter_subspace_max = 10
1714 0 : n_locked = 0
1715 0 : n_locked_prev = 0
1716 0 : nrows = rows(chebfi%xXColsRows)
1717 0 : nrows_mpi = rows(chebfi%X)
1718 :
1719 : ! Initialize active vectors (in colsrows distribution)
1720 0 : call xgBlock_setBlock(chebfi%xXColsRows , X_part%colsrows_active , nrows, chebfi%bandpp)
1721 0 : call xgBlock_setBlock(chebfi%xAXColsRows, AX_part%colsrows_active, nrows, chebfi%bandpp)
1722 0 : call xgBlock_setBlock(chebfi%xBXColsRows, BX_part%colsrows_active, nrows, chebfi%bandpp)
1723 :
1724 0 : do iter_subspace=1, niter_subspace_max
1725 :
1726 0 : write(std_out,*) 'subspace iteration no=', iter_subspace
1727 0 : flush(std_out)
1728 :
1729 : ! ############################ Manage MPI ##############################
1730 0 : if (iter_subspace>1 .and. chebfi%paral_kgb==1) then
1731 :
1732 0 : if (n_locked > 0 .and. n_locked /= n_locked_prev ) then
1733 :
1734 : ! Assign active bands to processes
1735 0 : call bandPartitionInfo_initActiveSubcomm(mpiInfo, chebfi%comm_cols, chebfi%bandpp)
1736 :
1737 : ! Allocate colsrows distributions for active vectors (from linalg distribution)
1738 0 : call bandPartitionData_allocActiveColsrows(X_part, chebfi, mpiInfo, nspinor)
1739 0 : call bandPartitionData_allocActiveColsrows(AX_part, chebfi, mpiInfo, nspinor)
1740 0 : call bandPartitionData_allocActiveColsrows(BX_part, chebfi, mpiInfo, nspinor)
1741 :
1742 : ! Transpose active vectors to colsrows distribution
1743 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1744 0 : call timab(tim_transpose,1,tsec)
1745 0 : call xgTransposer_transpose(X_part%transposer_active, STATE_COLSROWS)
1746 0 : call xgTransposer_transpose(AX_part%transposer_active, STATE_COLSROWS)
1747 0 : call xgTransposer_transpose(BX_part%transposer_active, STATE_COLSROWS)
1748 0 : call timab(tim_transpose,2,tsec)
1749 : ABI_NVTX_END_RANGE()
1750 :
1751 : else
1752 :
1753 : ! Use global transposer
1754 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1755 0 : call timab(tim_transpose,1,tsec)
1756 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
1757 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
1758 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
1759 0 : call timab(tim_transpose,2,tsec)
1760 : ABI_NVTX_END_RANGE()
1761 : end if
1762 :
1763 : end if
1764 :
1765 : ! ############################ Filter active ##############################
1766 : ! ############################ column vectors ##############################
1767 0 : if (is_lowpass) then
1768 : call chebfi_lowpassFilterActive(chebfi,X_part%colsrows_active,AX_part%colsrows_active,&
1769 0 : BX_part%colsrows_active,eigen,lambda_minus,lambda_plus,getAX_BX,getBm1X)
1770 : else
1771 : call chebfi_bandpassFilterActive(chebfi,X_part%colsrows_active,AX_part%colsrows_active,&
1772 : BX_part%colsrows_active,lambda_minus,lambda_plus,mineig_global,maxeig_global,&
1773 0 : getAX_BX,getBm1X)
1774 : end if
1775 : ! todo prevent bug X_next, X_prev might need resetting
1776 :
1777 : ! ############################ Manage MPI ##############################
1778 0 : if (chebfi%paral_kgb==1) then
1779 :
1780 0 : if (n_locked > 0) then
1781 :
1782 0 : call timab(tim_transpose,1,tsec)
1783 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1784 0 : call xgTransposer_transpose(X_part%transposer_active, STATE_LINALG)
1785 0 : call xgTransposer_transpose(AX_part%transposer_active, STATE_LINALG)
1786 0 : call xgTransposer_transpose(BX_part%transposer_active, STATE_LINALG)
1787 : ABI_NVTX_END_RANGE()
1788 0 : call timab(tim_transpose,2,tsec)
1789 :
1790 : else
1791 :
1792 : ! Use global transposer
1793 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1794 0 : call timab(tim_transpose,1,tsec)
1795 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_LINALG)
1796 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_LINALG)
1797 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_LINALG)
1798 0 : call timab(tim_transpose,2,tsec)
1799 : ABI_NVTX_END_RANGE()
1800 :
1801 : ! Initialize active vectors (in linalg distribution)
1802 0 : call xgBlock_setBlock(chebfi%X , X_part%linalg_active , nrows_mpi, neigenpairs)
1803 0 : call xgBlock_setBlock(chebfi%AX%self, AX_part%linalg_active, nrows_mpi, neigenpairs)
1804 0 : call xgBlock_setBlock(chebfi%BX%self, BX_part%linalg_active, nrows_mpi, neigenpairs)
1805 :
1806 : end if
1807 :
1808 : end if
1809 :
1810 0 : write(std_out,'(a,i6,i6)') 'chebfi%X # rows # cols ', rows(chebfi%X), cols(chebfi%X)
1811 0 : write(std_out,'(a,i6,i6)') 'X_active(linalg) # rows # cols ', nrows_mpi, neigenpairs
1812 0 : flush(std_out)
1813 :
1814 0 : write(std_out,*) 'getid X ', xgBlock_getid(chebfi%X)
1815 0 : write(std_out,*) 'getid AX ', xgBlock_getid(chebfi%AX%self)
1816 0 : write(std_out,*) 'getid BX ', xgBlock_getid(chebfi%BX%self)
1817 0 : write(std_out,*) 'space(X) ', space(chebfi%X)
1818 0 : write(std_out,*) 'comm(X) ', comm(chebfi%X)
1819 0 : write(std_out,*) 'getid Xpart ', xgBlock_getid(X_part%linalg_active)
1820 0 : write(std_out,*) 'getid AXpart', xgBlock_getid(AX_part%linalg_active)
1821 0 : write(std_out,*) 'getid BXpart', xgBlock_getid(BX_part%linalg_active)
1822 0 : write(std_out,*) 'space(Xpart)', space(X_part%linalg_active)
1823 0 : write(std_out,*) 'comm(Xpart) ', comm(X_part%linalg_active)
1824 0 : flush(std_out)
1825 :
1826 : ! ####################### Deflate+Ortho+Residual ##############################
1827 : ! ####################### for active column vectors ##############################
1828 :
1829 0 : if (n_locked > 0) then
1830 0 : call chebfi_deflateActiveWrtLocked(chebfi, n_locked, X_part, AX_part, BX_part)
1831 : end if
1832 :
1833 : call xg_Borthonormalize(X_part%linalg_active,BX_part%linalg_active,ierr,1,&
1834 0 : chebfi%gpu_option,AX=AX_part%linalg_active)
1835 :
1836 0 : call chebfi_getSubspaceResidual(chebfi, resid_active%self)
1837 :
1838 :
1839 :
1840 0 : call chebfi_swapConvergedVectors(chebfi, resid_active%self, 1e-3_dp, n_locked)
1841 : ! todo check if residual of wanted bands has converged if yes skip parts
1842 : ! use k_conv as number of wanted bands
1843 0 : write(std_out,*) 'keeping', k_conv, 'bands'
1844 0 : flush(std_out)
1845 :
1846 0 : if (chebfi%paral_kgb==1 .and. iter_subspace>1 .and. n_locked > 0 .and. n_locked /= n_locked_prev) then
1847 0 : call bandPartitionData_setLinalg(chebfi%X , X_part , n_locked)
1848 0 : call bandPartitionData_setLinalg(chebfi%AX%self, AX_part, n_locked)
1849 0 : call bandPartitionData_setLinalg(chebfi%BX%self, BX_part, n_locked)
1850 :
1851 : ! free colsrows active objects
1852 0 : call xgTransposer_free(X_part%transposer_active)
1853 0 : call xgTransposer_free(AX_part%transposer_active)
1854 0 : call xgTransposer_free(BX_part%transposer_active)
1855 : end if
1856 :
1857 0 : n_locked_prev = n_locked
1858 0 : mpiInfo%n_locked = n_locked
1859 0 : mpiInfo%n_active = chebfi%neigenpairs - n_locked
1860 :
1861 : ! todo debug by recomputing the getSubspaceResidual for locked vectors and verify it is smaller than tol
1862 :
1863 : end do
1864 :
1865 : ! Transfer active data to full workspace
1866 0 : call xgBlock_copy(X_part%linalg_active, chebfi%X)
1867 0 : call xgBlock_copy(AX_part%linalg_active, chebfi%BX%self)
1868 0 : call xgBlock_copy(BX_part%linalg_active, chebfi%BX%self)
1869 :
1870 : ! Apply Rayleigh-Ritz to active MPI Linalg row-block
1871 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RR)
1872 : call xg_RayleighRitz(chebfi%X,chebfi%AX%self,chebfi%BX%self,eigen,ierr,0,tim_RR,&
1873 0 : chebfi%gpu_option,solve_ax_bx=.true.)
1874 : ABI_NVTX_END_RANGE()
1875 :
1876 : !write(std_out,*) 'id of X, (after RR) ncols=', xgBlock_getId(chebfi%X), cols(chebfi%X)
1877 :
1878 0 : if ( ierr /= 0 ) then
1879 0 : ABI_WARNING("RayleighRitz did not work")
1880 : end if
1881 :
1882 : ! Compute residual norm *squared*
1883 0 : if (chebfi%paw) then
1884 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%BX%self,chebfi%AX%self)
1885 : else
1886 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%X,chebfi%AX%self)
1887 : end if
1888 0 : call xgBlock_colwiseNorm2(chebfi%AX%self,residu) ! performs MPI comm
1889 :
1890 : ! MPI Transpose to recover colsrows state
1891 0 : call timab(tim_transpose,1,tsec)
1892 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1893 0 : if (chebfi%paral_kgb == 1) then
1894 0 : call xmpi_barrier(chebfi%spacecom)
1895 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
1896 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
1897 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
1898 0 : if (xmpi_comm_size(chebfi%spacecom) == 1) then
1899 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
1900 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
1901 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
1902 : end if
1903 : else
1904 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
1905 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
1906 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
1907 : end if
1908 : ABI_NVTX_END_RANGE()
1909 0 : call timab(tim_transpose,2,tsec)
1910 :
1911 : ! Copy in ColsRows representation
1912 0 : call xgBlock_copy(chebfi%xXColsRows, X0)
1913 :
1914 0 : if (cols(X0) /= chebfi%bandpp) then
1915 0 : ABI_ERROR('wrong colsrows representation')
1916 : end if
1917 0 : write(std_out,'(a,i6,i6,i6)') 'local # proc has # rows cols ', xmpi_comm_rank(chebfi%spacecom), rows(X0), cols(X0)
1918 :
1919 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
1920 : if (chebfi%gpu_option==ABI_GPU_KOKKOS) then
1921 : call gpu_device_synchronize()
1922 : end if
1923 : #endif
1924 :
1925 : ! Free transposer objects
1926 0 : if (chebfi%paral_kgb == 1) then
1927 0 : call xgTransposer_free(chebfi%xgTransposerX)
1928 0 : call xgTransposer_free(chebfi%xgTransposerAX)
1929 0 : call xgTransposer_free(chebfi%xgTransposerBX)
1930 : end if
1931 :
1932 : ! Free temporary memory
1933 0 : ABI_SFREE(nrowsLinalg)
1934 0 : call xg_free(resid_active)
1935 :
1936 0 : end subroutine chebfi_runSubspaceIteration
1937 : !!***
1938 :
1939 : !----------------------------------------------------------------------
1940 :
1941 : !!****f* m_chebfi/chebfi_runSubspaceIterationDummy
1942 : !! NAME
1943 : !! chebfi_runSubspaceIterationDummy
1944 : !!
1945 : !! FUNCTION
1946 : !! Applies subspace iteration by Chebyshev polynomial filtering on slice
1947 : !! then extracts eigenvectors using Rayleigh-Ritz. This execution is completely
1948 : !! local on slice processors. Each slice does not see others.
1949 : !!
1950 : !! NOTES
1951 : !! Restart logic in subspace iteration:
1952 : !! Allocate constructor (state=colsrows)
1953 : !! while (convergence not reached):
1954 : !! | if (state=linalg) Transpose
1955 : !! | Filter
1956 : !! | Transpose (state=linalg)
1957 : !! | Orthogonalize
1958 : !! Rayleigh-Ritz
1959 : !!
1960 : !! INPUTS
1961 : !! X0 input vectors distibuted by colsrows along slice processors
1962 : !!
1963 : !! OUTPUTS
1964 : !! X0 (in-place) converged vectors distributed by linalg along slice processors
1965 : !!
1966 : !! dev=============
1967 : !! - version 1: count number of converged vectors within iteration
1968 : !!
1969 : !! SOURCE
1970 :
1971 0 : subroutine chebfi_runSubspaceIterationDummy(chebfi,X0,getAX_BX,getBm1X,eigen,residu,nspinor,&
1972 : mineig_global,maxeig_global,lambda_minus,lambda_plus,is_lowpass,k_rank,nrows_blockrows)
1973 :
1974 : implicit none
1975 :
1976 : !Arguments ------------------------------------
1977 : type(chebfi_t) , intent(inout) :: chebfi
1978 : type(xgBlock_t), intent(inout) :: X0
1979 : type(xgBlock_t), intent(inout) :: eigen
1980 : type(xgBlock_t), intent(inout) :: residu
1981 : integer , intent(in ) :: nspinor
1982 : integer , intent(in ) :: k_rank
1983 : integer, pointer, intent(in ) :: nrows_blockrows(:)
1984 : real(dp) , intent(in ) :: mineig_global
1985 : real(dp) , intent(in ) :: maxeig_global
1986 : real(dp) , intent(in ) :: lambda_minus
1987 : real(dp) , intent(in ) :: lambda_plus
1988 : logical , intent(in ) :: is_lowpass
1989 : interface
1990 : subroutine getAX_BX(X,AX,BX)
1991 : use m_xg, only : xgBlock_t
1992 : type(xgBlock_t), intent(inout) :: X
1993 : type(xgBlock_t), intent(inout) :: AX
1994 : type(xgBlock_t), intent(inout) :: BX
1995 : end subroutine getAX_BX
1996 : end interface
1997 : interface
1998 : subroutine getBm1X(X,Bm1X)
1999 : use m_xg, only : xgBlock_t
2000 : type(xgBlock_t), intent(inout) :: X
2001 : type(xgBlock_t), intent(inout) :: Bm1X
2002 : end subroutine getBm1X
2003 : end interface
2004 :
2005 : !Local variables-------------------------------
2006 : integer :: iter_subspace, niter_subspace_max, n_locked, n_active, n_locked_prev
2007 : integer :: spacedim, neigenpairs, num_proc, nrows, ierr
2008 : !type(xg_t) :: X_k
2009 : type(xg_t) :: resid_active
2010 : type(xg_t) :: X_locked, AX_locked, BX_locked ! spaces needed for deflation in linalg distr
2011 : type(xgBlock_t) :: chebfi_X_locked, chebfi_AX_locked, chebfi_BX_locked ! pointers used for copy
2012 : type(xgBlock_t) :: X_active, AX_active, BX_active
2013 : ! Arrays
2014 0 : integer, target, allocatable :: nrowsLinalg(:)
2015 : integer, pointer :: nrowsLinalg_ptr(:) => null()
2016 : real(dp) :: tsec(2)
2017 :
2018 : ! *********************************************************************
2019 :
2020 0 : if (chebfi%from_linalg) then
2021 0 : ABI_ERROR("chebfi should be from colsrows")
2022 : end if
2023 :
2024 0 : spacedim = chebfi%spacedim
2025 0 : neigenpairs = chebfi%neigenpairs
2026 0 : num_proc = xmpi_comm_size(chebfi%spacecom)
2027 0 : chebfi%eigenvalues = eigen
2028 :
2029 0 : ABI_MALLOC_IFNOT(nrowsLinalg,(num_proc))
2030 0 : nrowsLinalg_ptr => nrowsLinalg
2031 0 : nrowsLinalg = nrows_blockrows
2032 :
2033 0 : call xg_init(resid_active, SPACE_R, neigenpairs, 1, gpu_option=chebfi%gpu_option)
2034 :
2035 0 : write(std_out,*) 'wanted rank=', k_rank
2036 :
2037 : !write(std_out,*) 'getid inside runSlice xXColsRows', xgBlock_getId(chebfi%xXColsRows)
2038 : !write(std_out,*) 'getid inside runSlice xAXColsRows', xgBlock_getId(chebfi%xAXColsRows)
2039 : !flush(std_out)
2040 :
2041 : ! Initialize values of AX and BX
2042 0 : call timab(tim_getAX_BX,1,tsec)
2043 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2044 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
2045 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
2046 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
2047 : ABI_NVTX_END_RANGE()
2048 0 : call timab(tim_getAX_BX,2,tsec)
2049 :
2050 : ! Allocate memory for linalg distribution
2051 0 : call timab(tim_transpose,1,tsec)
2052 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
2053 0 : if (chebfi%paral_kgb==1) then
2054 :
2055 : ! Allocate chebfi%X
2056 : call xgTransposer_constructor(chebfi%xgTransposerX,chebfi%X,chebfi%xXColsRows,nspinor,&
2057 : STATE_COLSROWS,TRANS_ALL2ALL,chebfi%comm_rows,chebfi%comm_cols,0,0,chebfi%me_g0,&
2058 : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit,&
2059 0 : custom_ncolsColsRows=.true.,nrowsLinalg_sub=nrowsLinalg_ptr)
2060 : ! Note: bandpp is custom because it is created from resource allocator
2061 :
2062 : ! Allocate chebfi%AX, chebfi%BX
2063 : call xgTransposer_copyConstructor(chebfi%xgTransposerAX,chebfi%xgTransposerX,&
2064 0 : chebfi%AX%self,chebfi%xAXColsRows,STATE_COLSROWS)
2065 : call xgTransposer_copyConstructor(chebfi%xgTransposerBX,chebfi%xgTransposerX,&
2066 0 : chebfi%BX%self,chebfi%xBXColsRows,STATE_COLSROWS)
2067 : ! Note: at this point chebfi%AX and chebfi%BX are empty. Must transpose
2068 : ! to fill with correct values.
2069 :
2070 : ! todo use copy constructor to create an object by copying an existing object
2071 : ! actually copy constructor *allocates* memory for chebfi%AX. Write a version
2072 : ! that does not allocate memory and only reassigns pointers. Can pointers be reassigned
2073 : ! directly then used in the global constructor? Perform tests.
2074 :
2075 0 : chebfi%xgTransposerX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
2076 0 : chebfi%xgTransposerAX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
2077 0 : chebfi%xgTransposerBX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
2078 :
2079 : else
2080 0 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
2081 0 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
2082 0 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
2083 : end if
2084 0 : call timab(tim_transpose,2,tsec)
2085 : ABI_NVTX_END_RANGE()
2086 :
2087 0 : write(std_out,*) '@dummy starting filter in proc', xmpi_comm_rank(chebfi%spacecom)
2088 0 : flush(std_out)
2089 :
2090 0 : niter_subspace_max = 5
2091 0 : n_locked = 0
2092 0 : n_locked_prev = 0
2093 :
2094 : ! chebfi%X contains active vectors
2095 : ! X_lock contains locked vectors
2096 :
2097 : !xX_active = chebfi%xXColsRows
2098 : !xAX_active = chebfi%xAXColsRows
2099 : !xBX_active = chebfi%xBXColsRows
2100 :
2101 0 : do iter_subspace=1, niter_subspace_max
2102 :
2103 0 : write(std_out,*) 'subspace iteration no=', iter_subspace
2104 0 : flush(std_out)
2105 :
2106 0 : if (iter_subspace>1 .and. chebfi%paral_kgb==1) then
2107 :
2108 0 : call timab(tim_transpose,1,tsec)
2109 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
2110 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
2111 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
2112 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
2113 0 : call timab(tim_transpose,2,tsec)
2114 : ABI_NVTX_END_RANGE()
2115 :
2116 : !if (n_locked>1) then
2117 : ! ! n_locked_mpi each rank has different number of locked vectors todo (find this in v2)
2118 : ! call xgBlock_setBlock(chebfi%xXColsRows , xX_active , spacedim, n_locked_mpi)
2119 : ! call xgBlock_setBlock(chebfi%xAXColsRows, xAX_active, spacedim, n_locked_mpi)
2120 : ! call xgBlock_setBlock(chebfi%xBXColsRows, xBX_active, spacedim, n_locked_mpi)
2121 : !end if
2122 : end if
2123 : ! todo this will transform the active+locked in colsrows. Locked is not necessary
2124 : ! therefore try to communicate less data is possible.
2125 :
2126 : ! ############################ Filter active ##############################
2127 : ! ############################ column vectors ##############################
2128 0 : if (is_lowpass) then
2129 0 : call chebfi_lowpassFilter(chebfi,eigen,lambda_minus,lambda_plus,getAX_BX,getBm1X)
2130 : else
2131 : call chebfi_bandpassFilter(chebfi,lambda_minus,lambda_plus,mineig_global,&
2132 0 : maxeig_global,getAX_BX,getBm1X)
2133 : end if
2134 :
2135 0 : write(std_out,'(a,i6,i6)') 'chebfi%xXColsRows # rows # cols ', &
2136 0 : rows(chebfi%xXColsRows), cols(chebfi%xXColsRows)
2137 0 : flush(std_out)
2138 :
2139 0 : if (chebfi%paral_kgb==1) then
2140 0 : call timab(tim_transpose,1,tsec)
2141 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
2142 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_LINALG)
2143 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_LINALG)
2144 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_LINALG)
2145 0 : call timab(tim_transpose,2,tsec)
2146 : ABI_NVTX_END_RANGE()
2147 : end if
2148 :
2149 0 : if (n_locked > 0) then
2150 : ! undo the effect of filtering in locked vectors
2151 0 : call xgBlock_setBlock(chebfi%X, chebfi_X_locked, nrows, n_locked)
2152 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi_AX_locked, nrows, n_locked)
2153 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi_BX_locked, nrows, n_locked)
2154 0 : call xgBlock_copy(X_locked%self, chebfi_X_locked)
2155 0 : call xgBlock_copy(AX_locked%self, chebfi_AX_locked)
2156 0 : call xgBlock_copy(BX_locked%self, chebfi_BX_locked)
2157 : end if
2158 :
2159 0 : write(std_out,'(a,i6,i6)') 'chebfi%X # rows # cols ', rows(chebfi%X), cols(chebfi%X)
2160 0 : write(std_out,*) xgBlock_getid(chebfi%X)
2161 0 : flush(std_out)
2162 :
2163 : ! recover active vectors in linalg
2164 0 : n_active = chebfi%neigenpairs - n_locked
2165 0 : nrows = rows(chebfi%X)
2166 :
2167 0 : write(std_out,*) 'n_active=', n_active; flush(std_out)
2168 0 : call xgBlock_setBlock(chebfi%X, X_active, nrows, n_active, fcol=n_locked+1)
2169 0 : call xgBlock_setBlock(chebfi%AX%self, AX_active, nrows, n_active, fcol=n_locked+1)
2170 0 : call xgBlock_setBlock(chebfi%BX%self, BX_active, nrows, n_active, fcol=n_locked+1)
2171 :
2172 0 : write(std_out,*) 'start deflation'
2173 0 : write(std_out,*) xgBlock_getid(chebfi%X), comm(chebfi%X), me_g0(chebfi%X), cols(chebfi%X), space(chebfi%X)
2174 0 : write(std_out,*) xgBlock_getid(X_active), comm(X_active), me_g0(X_active), cols(X_active), space(X_active)
2175 0 : flush(std_out)
2176 :
2177 0 : if (n_locked>0) then
2178 0 : call chebfi_deflateWrtLocked(chebfi, n_locked, X_locked%self, BX_locked%self)
2179 : !call chebfi_deflateActiveWrtLocked(chebfi, n_locked, X_part, AX_part, BX_part)
2180 : end if
2181 :
2182 0 : write(std_out,*) 'exit deflation'; flush(std_out)
2183 0 : write(std_out,*) 'start ortho'; flush(std_out)
2184 :
2185 0 : write(std_out,*) 'cols=', cols(chebfi%X); flush(std_out)
2186 :
2187 0 : call xg_Borthonormalize(X_active,BX_active,ierr,1,chebfi%gpu_option,AX=AX_active)
2188 : !call xg_Borthonormalize(chebfi%X,chebfi%BX%self,ierr,1,chebfi%gpu_option,AX=chebfi%AX%self)
2189 :
2190 0 : write(std_out,*) 'exit ortho'; flush(std_out)
2191 0 : write(std_out,*) 'start residual'; flush(std_out)
2192 :
2193 0 : call chebfi_getSubspaceResidual(chebfi, resid_active%self)
2194 :
2195 : ! prevent errors in pointer corruption to non-contiguous memory
2196 0 : call xgBlock_setBlock(chebfi%X, X_active, nrows, 1)
2197 0 : call xgBlock_setBlock(chebfi%AX%self, AX_active, nrows, 1)
2198 0 : call xgBlock_setBlock(chebfi%BX%self, BX_active, nrows, 1)
2199 0 : call xgBlock_setBlock(chebfi%X, chebfi_X_locked, nrows, 1)
2200 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi_BX_locked, nrows, 1)
2201 :
2202 0 : call chebfi_swapConvergedVectors(chebfi, resid_active%self, 1e-3_dp, n_locked)
2203 :
2204 0 : if (n_locked>0) then
2205 0 : if (n_locked /= n_locked_prev) then
2206 :
2207 0 : if (n_locked_prev>0) then
2208 0 : call xg_free(X_locked)
2209 0 : call xg_free(AX_locked)
2210 0 : call xg_free(BX_locked)
2211 : end if
2212 :
2213 : ! Allocate locked vectors in linalg representation
2214 : call xg_init(X_locked, chebfi%space, nrows, n_locked, &
2215 0 : chebfi%spacecom, me_g0=chebfi%me_g0, gpu_option=chebfi%gpu_option)
2216 :
2217 : ! same for AX
2218 : call xg_init(AX_locked, chebfi%space, nrows, n_locked, &
2219 0 : chebfi%spacecom, me_g0=chebfi%me_g0, gpu_option=chebfi%gpu_option)
2220 :
2221 : ! same for BX
2222 : call xg_init(BX_locked, chebfi%space, nrows, n_locked, &
2223 0 : chebfi%spacecom, me_g0=chebfi%me_g0, gpu_option=chebfi%gpu_option)
2224 : end if
2225 :
2226 0 : call xgBlock_setBlock(chebfi%X, chebfi_X_locked, nrows, n_locked)
2227 0 : call xgBlock_copy(chebfi_X_locked, X_locked%self)
2228 :
2229 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi_AX_locked, nrows, n_locked)
2230 0 : call xgBlock_copy(chebfi_AX_locked, AX_locked%self)
2231 :
2232 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi_BX_locked, nrows, n_locked)
2233 0 : call xgBlock_copy(chebfi_BX_locked, BX_locked%self)
2234 : end if
2235 :
2236 : ! todo debug by recomputing the getSubspaceResidual for locked vectors and verify it is smaller than tol
2237 :
2238 : ! lock vectors in small residual
2239 : ! actually redefine pointers Xactive and Xlock pointing to X that's all
2240 :
2241 0 : n_locked_prev = n_locked
2242 :
2243 : end do
2244 :
2245 0 : call xg_free(X_locked)
2246 0 : call xg_free(AX_locked)
2247 0 : call xg_free(BX_locked)
2248 :
2249 : ! Apply Rayleigh-Ritz to active MPI Linalg row-block
2250 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RR)
2251 : call xg_RayleighRitz(chebfi%X,chebfi%AX%self,chebfi%BX%self,eigen,ierr,0,tim_RR,&
2252 0 : chebfi%gpu_option,solve_ax_bx=.true.)
2253 : ABI_NVTX_END_RANGE()
2254 :
2255 : !write(std_out,*) 'id of X, (after RR) ncols=', xgBlock_getId(chebfi%X), cols(chebfi%X)
2256 :
2257 0 : if ( ierr /= 0 ) then
2258 0 : ABI_WARNING("RayleighRitz did not work")
2259 : else
2260 : !write(std_out,*) 'is lowpass=', is_lowpass
2261 : !write(std_out,*) 'chebfi%eigenvalues after RR'
2262 : !call xgBlock_print(chebfi%eigenvalues,std_out)
2263 : !flush(std_out)
2264 : end if
2265 :
2266 : ! Compute residual norm *squared*
2267 0 : if (chebfi%paw) then
2268 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%BX%self,chebfi%AX%self)
2269 : else
2270 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%X,chebfi%AX%self)
2271 : end if
2272 0 : call xgBlock_colwiseNorm2(chebfi%AX%self,residu) ! performs MPI comm
2273 :
2274 : !write(std_out,*) 'max colwise residual norm squared='; call xgBlock_print(residu, std_out)
2275 : !flush(std_out)
2276 :
2277 : ! Copy in Linalg representation (see chebfi_run, kept for reference)
2278 : ! call xgBlock_copy(chebfi%X,X0)
2279 :
2280 : ! MPI Transpose to recover colsrows state (X only)
2281 0 : call timab(tim_transpose,1,tsec)
2282 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
2283 0 : if (chebfi%paral_kgb == 1) then
2284 0 : call xmpi_barrier(chebfi%spacecom)
2285 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
2286 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
2287 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
2288 0 : if (xmpi_comm_size(chebfi%spacecom) == 1) then
2289 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
2290 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
2291 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
2292 : end if
2293 : else
2294 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
2295 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
2296 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
2297 : end if
2298 : ABI_NVTX_END_RANGE()
2299 0 : call timab(tim_transpose,2,tsec)
2300 :
2301 : ! Copy in ColsRows representation
2302 0 : call xgBlock_copy(chebfi%xXColsRows, X0)
2303 :
2304 0 : if (cols(X0) /= chebfi%bandpp) then
2305 0 : ABI_ERROR('wrong colsrows representation')
2306 : end if
2307 0 : write(std_out,'(a,i6,i6,i6)') 'local # proc has # rows cols ', xmpi_comm_rank(chebfi%spacecom), rows(X0), cols(X0)
2308 :
2309 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
2310 : if (chebfi%gpu_option==ABI_GPU_KOKKOS) then
2311 : call gpu_device_synchronize()
2312 : end if
2313 : #endif
2314 :
2315 : ! Free transposer objects
2316 0 : if (chebfi%paral_kgb == 1) then
2317 0 : call xgTransposer_free(chebfi%xgTransposerX)
2318 0 : call xgTransposer_free(chebfi%xgTransposerAX)
2319 0 : call xgTransposer_free(chebfi%xgTransposerBX)
2320 : end if
2321 :
2322 : ! Free temporary memory
2323 0 : ABI_SFREE(nrowsLinalg)
2324 0 : call xg_free(resid_active)
2325 :
2326 0 : end subroutine chebfi_runSubspaceIterationDummy
2327 : !!***
2328 :
2329 : !----------------------------------------------------------------------
2330 :
2331 : !!****f* m_chebfi/chebfi_lowpassFilter
2332 : !! NAME
2333 : !! chebfi_lowpassFilter
2334 : !!
2335 : !! FUNCTION
2336 : !! Apply Lowpass filter using Chebyshev polynomial on a set of vectors.
2337 : !! Amplifies interval [-oo, lambda_minus) and diminishes [lambda_minus,lambda_plus).
2338 : !!
2339 : !! INPUTS
2340 : !! chebfi <type(chebfi_t)>=memory workspace used to apply filter
2341 : !! eigen= Rayleigh quotients to use in amplification of chebfi%xXColsRows
2342 : !! lambda_minus= lower bound of interval to diminish
2343 : !! lambda_plus= upper bound of interval to diminish
2344 : !! getAX_BX= pointer to the function giving A|X> and B|X>
2345 : !! A is typically the Hamiltonian H, and B the overlap operator S
2346 : !! getBm1X= pointer to the function giving B^-1|X>
2347 : !! B is typically the overlap operator S
2348 : !!
2349 : !! SIDE EFFECTS
2350 : !! chebfi%xXColsRows= Filtered vectors to use in Subspace iteration
2351 : !!
2352 : !! SOURCE
2353 :
2354 0 : subroutine chebfi_lowpassFilter(chebfi,eigen,lambda_minus,lambda_plus,getAX_BX,getBm1X)
2355 :
2356 : implicit none
2357 :
2358 : ! Arguments ------------------------------------
2359 : type(chebfi_t), intent(inout) :: chebfi
2360 : type(xgBlock_t), intent(inout) :: eigen
2361 : real(dp), intent(in) :: lambda_minus
2362 : real(dp), intent(in) :: lambda_plus
2363 : interface
2364 : subroutine getAX_BX(X,AX,BX)
2365 : use m_xg, only : xgBlock_t
2366 : type(xgBlock_t), intent(inout) :: X
2367 : type(xgBlock_t), intent(inout) :: AX
2368 : type(xgBlock_t), intent(inout) :: BX
2369 : end subroutine getAX_BX
2370 : end interface
2371 : interface
2372 : subroutine getBm1X(X,Bm1X)
2373 : use m_xg, only : xgBlock_t
2374 : type(xgBlock_t), intent(inout) :: X
2375 : type(xgBlock_t), intent(inout) :: Bm1X
2376 : end subroutine getBm1X
2377 : end interface
2378 :
2379 : !Local variables-------------------------------
2380 : integer :: ideg
2381 : real(dp) :: center, radius, one_over_r, two_over_r
2382 : type(xg_t) :: DivResults ! Rayleigh quotients
2383 : ! Arrays
2384 : real(dp) :: tsec(2)
2385 0 : integer, allocatable :: ndeg_filter_bands(:)
2386 :
2387 : ! *********************************************************************
2388 :
2389 0 : if (chebfi%paral_kgb == 0) then
2390 0 : ABI_MALLOC_IFNOT(ndeg_filter_bands,(chebfi%neigenpairs))
2391 : else
2392 0 : ABI_MALLOC_IFNOT(ndeg_filter_bands,(chebfi%bandpp))
2393 : end if
2394 :
2395 : ! [lambda_minus,lambda_plus) is diminished using Chebyshev
2396 0 : write(std_out,*) '@lowpass lambda_minus=', lambda_minus
2397 0 : write(std_out,*) '@lowpass lambda_plus=', lambda_plus
2398 0 : flush(std_out)
2399 :
2400 : ! Filter parameters
2401 0 : ndeg_filter_bands(:) = chebfi%ndeg_filter
2402 0 : center = (lambda_plus + lambda_minus)*0.5
2403 0 : radius = (lambda_plus - lambda_minus)*0.5
2404 0 : one_over_r = 1/radius
2405 0 : two_over_r = 2/radius
2406 :
2407 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
2408 0 : do ideg = 0, chebfi%ndeg_filter - 1
2409 :
2410 : ! X_next=2/r*(AX_next-c*X_next)-X_prev
2411 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
2412 0 : call chebfi_computeNextOrderChebfiPolynom(chebfi, ideg, center, one_over_r, two_over_r, getBm1X)
2413 : ABI_NVTX_END_RANGE()
2414 :
2415 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
2416 0 : if (chebfi%paral_kgb == 0) then
2417 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%spacedim, chebfi%neigenpairs)
2418 : else
2419 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
2420 : end if
2421 : ABI_NVTX_END_RANGE()
2422 :
2423 : !A * Psi (=AX_next=A*X_next)
2424 0 : call timab(tim_getAX_BX,1,tsec)
2425 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2426 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
2427 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
2428 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
2429 : ABI_NVTX_END_RANGE()
2430 0 : call timab(tim_getAX_BX,2,tsec)
2431 :
2432 : end do ! ideg
2433 : ABI_NVTX_END_RANGE()
2434 :
2435 : ! Avoid overflow rescale
2436 0 : call chebfi_prepAmpfactor(chebfi, eigen, DivResults)
2437 0 : call chebfi_ampfactor(chebfi, DivResults%self, lambda_minus, lambda_plus, ndeg_filter_bands)
2438 :
2439 : ! Free temporary memory
2440 0 : call xg_free(DivResults)
2441 0 : ABI_SFREE(ndeg_filter_bands)
2442 :
2443 0 : end subroutine chebfi_lowpassFilter
2444 : !!***
2445 :
2446 : !----------------------------------------------------------------------
2447 :
2448 : !!****f* m_chebfi/chebfi_lowpassFilterActive
2449 : !! NAME
2450 : !! chebfi_lowpassFilterActive
2451 : !!
2452 : !! FUNCTION
2453 : !! Apply Lowpass filter on Active vectors only
2454 : !! TODO this function will replace lowpassFilter
2455 : !!
2456 : !! SOURCE
2457 :
2458 0 : subroutine chebfi_lowpassFilterActive(chebfi,X_active,AX_active,BX_active,eigen,&
2459 : lambda_minus,lambda_plus,getAX_BX,getBm1X)
2460 :
2461 : implicit none
2462 :
2463 : ! Arguments ------------------------------------
2464 : type(chebfi_t), intent(inout) :: chebfi
2465 : type(xgBlock_t), intent(inout) :: X_active
2466 : type(xgBlock_t), intent(inout) :: AX_active
2467 : type(xgBlock_t), intent(inout) :: BX_active
2468 : type(xgBlock_t), intent(inout) :: eigen
2469 : real(dp), intent(in) :: lambda_minus
2470 : real(dp), intent(in) :: lambda_plus
2471 : interface
2472 : subroutine getAX_BX(X,AX,BX)
2473 : use m_xg, only : xgBlock_t
2474 : type(xgBlock_t), intent(inout) :: X
2475 : type(xgBlock_t), intent(inout) :: AX
2476 : type(xgBlock_t), intent(inout) :: BX
2477 : end subroutine getAX_BX
2478 : end interface
2479 : interface
2480 : subroutine getBm1X(X,Bm1X)
2481 : use m_xg, only : xgBlock_t
2482 : type(xgBlock_t), intent(inout) :: X
2483 : type(xgBlock_t), intent(inout) :: Bm1X
2484 : end subroutine getBm1X
2485 : end interface
2486 :
2487 : !Local variables-------------------------------
2488 : integer :: ideg
2489 : real(dp) :: center, radius, one_over_r, two_over_r
2490 : type(xg_t) :: DivResults ! Rayleigh quotients
2491 : ! Arrays
2492 : real(dp) :: tsec(2)
2493 0 : integer, allocatable :: ndeg_filter_bands(:)
2494 :
2495 : ! *********************************************************************
2496 :
2497 0 : if (chebfi%paral_kgb == 0) then
2498 0 : ABI_MALLOC_IFNOT(ndeg_filter_bands,(chebfi%neigenpairs))
2499 : else
2500 0 : ABI_MALLOC_IFNOT(ndeg_filter_bands,(chebfi%bandpp))
2501 : end if
2502 :
2503 : ! [lambda_minus,lambda_plus) is diminished using Chebyshev
2504 0 : write(std_out,*) '@lowpass lambda_minus=', lambda_minus
2505 0 : write(std_out,*) '@lowpass lambda_plus=', lambda_plus
2506 0 : flush(std_out)
2507 :
2508 : ! Filter parameters
2509 0 : ndeg_filter_bands(:) = chebfi%ndeg_filter
2510 0 : center = (lambda_plus + lambda_minus)*0.5
2511 0 : radius = (lambda_plus - lambda_minus)*0.5
2512 0 : one_over_r = 1/radius
2513 0 : two_over_r = 2/radius
2514 :
2515 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
2516 0 : do ideg = 0, chebfi%ndeg_filter - 1
2517 :
2518 : ! X_next=2/r*(AX_next-c*X_next)-X_prev
2519 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
2520 0 : call chebfi_computeNextOrderChebfiPolynom(chebfi, ideg, center, one_over_r, two_over_r, getBm1X)
2521 : ABI_NVTX_END_RANGE()
2522 :
2523 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
2524 0 : if (chebfi%paral_kgb == 0) then
2525 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%spacedim, chebfi%neigenpairs)
2526 : else
2527 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
2528 : end if
2529 : ABI_NVTX_END_RANGE()
2530 :
2531 : !A * Psi (=AX_next=A*X_next)
2532 0 : call timab(tim_getAX_BX,1,tsec)
2533 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2534 0 : call getAX_BX(X_active, AX_active, BX_active)
2535 0 : call xgBlock_zero_im_g0(AX_active)
2536 0 : call xgBlock_zero_im_g0(BX_active)
2537 : ABI_NVTX_END_RANGE()
2538 0 : call timab(tim_getAX_BX,2,tsec)
2539 :
2540 : end do ! ideg
2541 : ABI_NVTX_END_RANGE()
2542 :
2543 : ! Avoid overflow rescale
2544 0 : call chebfi_prepAmpfactor(chebfi, eigen, DivResults)
2545 0 : call chebfi_ampfactor(chebfi, DivResults%self, lambda_minus, lambda_plus, ndeg_filter_bands)
2546 :
2547 : ! Free temporary memory
2548 0 : call xg_free(DivResults)
2549 0 : ABI_SFREE(ndeg_filter_bands)
2550 :
2551 0 : end subroutine chebfi_lowpassFilterActive
2552 : !!***
2553 :
2554 : !----------------------------------------------------------------------
2555 :
2556 : !!****f* m_chebfi/chebfi_bandpassFilter
2557 : !! NAME
2558 : !! chebfi_bandpassFilter
2559 : !!
2560 : !! FUNCTION
2561 : !! Apply Bandpass filter using Chebyshev-Jackson polynomial, that is an
2562 : !! approximation of Heaviside step function by a Chebyshev expansion plus
2563 : !! Jackson damping to reduce oscillations, applied on a set of vectors.
2564 : !! Amplifies interval [lambda_minus, lambda_plus).
2565 : !!
2566 : !! INPUTS
2567 : !! chebfi <type(chebfi_t)>=memory workspace used to apply filter
2568 : !! eigen= Rayleigh quotients to use in amplification of chebfi%xXColsRows
2569 : !! lambda_minus= lower bound of interval to amplify
2570 : !! lambda_plus= upper bound of interval to amplify
2571 : !! mineig_global= used to rescale to [-1,1), will le -1
2572 : !! maxeig_global= used to rescale to [-1,1), will be 1
2573 : !! getAX_BX= pointer to the function giving A|X> and B|X>
2574 : !! A is typically the Hamiltonian H, and B the overlap operator S
2575 : !! getBm1X= pointer to the function giving B^-1|X>
2576 : !! B is typically the overlap operator S
2577 : !!
2578 : !! SIDE EFFECTS
2579 : !! chebfi%xXColsRows= Filtered vectors to use in Subspace iteration
2580 : !!
2581 : !! SOURCE
2582 :
2583 0 : subroutine chebfi_bandpassFilter(chebfi,lambda_minus,lambda_plus,mineig_global,&
2584 : maxeig_global,getAX_BX,getBm1X)
2585 :
2586 : implicit none
2587 :
2588 : ! Arguments ------------------------------------
2589 : type(chebfi_t), intent(inout) :: chebfi
2590 : !type(xgBlock_t), intent(inout) :: eigen
2591 : real(dp), intent(in) :: lambda_minus
2592 : real(dp), intent(in) :: lambda_plus
2593 : real(dp), intent(in) :: mineig_global
2594 : real(dp), intent(in) :: maxeig_global
2595 : interface
2596 : subroutine getAX_BX(X,AX,BX)
2597 : use m_xg, only : xgBlock_t
2598 : type(xgBlock_t), intent(inout) :: X
2599 : type(xgBlock_t), intent(inout) :: AX
2600 : type(xgBlock_t), intent(inout) :: BX
2601 : end subroutine getAX_BX
2602 : end interface
2603 : interface
2604 : subroutine getBm1X(X,Bm1X)
2605 : use m_xg, only : xgBlock_t
2606 : type(xgBlock_t), intent(inout) :: X
2607 : type(xgBlock_t), intent(inout) :: Bm1X
2608 : end subroutine getBm1X
2609 : end interface
2610 :
2611 : ! Local variables-------------------------------
2612 : integer :: ndeg, n
2613 : real(dp) :: center, radius, one_over_r, two_over_r
2614 : real(dp) :: ls, us, cdeg, mu, damp
2615 : type(xg_t) :: Heaviside
2616 : type(xg_t) :: DivResults ! Rayleigh quotients
2617 : real(dp) :: tsec(2)
2618 :
2619 : ! *********************************************************************
2620 :
2621 0 : ndeg = chebfi%ndeg_filter
2622 :
2623 : ! Allocate memory for Chebyshev expansion of Heaviside step function
2624 : call xg_init(Heaviside, chebfi%space, chebfi%total_spacedim, chebfi%bandpp, chebfi%spacecom, &
2625 0 : gpu_option=chebfi%gpu_option)
2626 : !call xgBlock_zero(Heaviside%self)
2627 :
2628 0 : write(std_out,*) '@bandpass lambda_minus=', lambda_minus
2629 0 : write(std_out,*) '@bandpass lambda_plus=', lambda_plus
2630 0 : write(std_out,*) '@bandpass mineig_global=', mineig_global
2631 0 : write(std_out,*) '@bandpass maxeig_global=', maxeig_global
2632 0 : flush(std_out)
2633 :
2634 : ! Filter parameters
2635 0 : center = (maxeig_global + mineig_global)*0.5
2636 0 : radius = (maxeig_global - mineig_global)*0.5
2637 0 : one_over_r = 1/radius
2638 0 : two_over_r = 2/radius
2639 :
2640 : ! Scaled slice bounds to be amplified
2641 0 : ls = (lambda_minus - center) / radius
2642 0 : us = (lambda_plus - center) / radius
2643 :
2644 : ! TODO IL 10/3/2025 Deflate vectors to reduce linear dependence: Y=X-(B-projection)
2645 : !call xg_Borthonormalize(chebfi%xXColsRows,chebfi%xBxColsRows,ierr,1,chebfi%gpu_option,AX=chebfi%xAXColsRows)
2646 :
2647 : ! Initialize expansion: Heaviside = mu(0)*damp(0)*X + Heaviside
2648 0 : cdeg = Pi/(ndeg+2)
2649 0 : mu = 1/Pi*(ACOS(ls)-ACOS(us))
2650 0 : damp = 1.d0 ! Jackson damping
2651 0 : call xgBlock_saxpy(Heaviside%self, mu*damp, chebfi%xXColsRows)
2652 :
2653 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
2654 0 : do n = 0, ndeg - 1
2655 :
2656 : ! X_next=2/r*(AX_next-c*X_next)-X_prev
2657 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
2658 0 : call chebfi_computeNextOrderChebfiPolynom(chebfi, n, center, one_over_r, two_over_r, getBm1X)
2659 : ABI_NVTX_END_RANGE()
2660 :
2661 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
2662 0 : if (chebfi%paral_kgb == 0) then
2663 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%spacedim, chebfi%neigenpairs)
2664 : else
2665 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
2666 : end if
2667 : ABI_NVTX_END_RANGE()
2668 :
2669 : ! Update expansion: Heaviside = damp(i+1)*mu(i+1)*X_next + Heaviside
2670 0 : mu = 2/Pi * (SIN((n+1)*ACOS(ls)) - SIN((n+1)*ACOS(us)))/(n+1)
2671 0 : damp = ((1 - (n+1)/(ndeg+2))*SIN(cdeg)*COS((n+1)*cdeg) + 1/(ndeg+2)*COS(cdeg)*SIN((n+1)*cdeg))/SIN(cdeg)
2672 0 : call xgBlock_saxpy(Heaviside%self, mu*damp, chebfi%xXColsRows)
2673 :
2674 : ! Store final expansion before exit, X_next=Heaviside
2675 0 : if (n==ndeg-1) then
2676 0 : call xgBlock_copy(Heaviside%self, chebfi%xXColsRows)
2677 : end if
2678 :
2679 : !A * Psi (=AX_next=A*X_next)
2680 0 : call timab(tim_getAX_BX,1,tsec)
2681 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2682 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
2683 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
2684 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
2685 : ABI_NVTX_END_RANGE()
2686 0 : call timab(tim_getAX_BX,2,tsec)
2687 :
2688 : end do ! end n
2689 : ABI_NVTX_END_RANGE()
2690 :
2691 : ! Avoid overflow
2692 : !call chebfi_prepAmpfactor(chebfi, eigen, DivResults)
2693 : !call chebfi_ampfactorBandpass(chebfi, DivResults%self, lambda_minus, lambda_plus, center, radius, ndeg)
2694 :
2695 : ! Free temporary memory
2696 0 : call xg_free(DivResults)
2697 0 : call xg_free(Heaviside)
2698 :
2699 0 : end subroutine chebfi_bandpassFilter
2700 : !!***
2701 :
2702 : !----------------------------------------------------------------------
2703 :
2704 : !!****f* m_chebfi/chebfi_bandpassFilterActive
2705 : !! NAME
2706 : !! chebfi_bandpassFilterActive
2707 : !!
2708 : !! FUNCTION
2709 : !! Apply Bandpass filter using Chebyshev-Jackson polynomial, that is an
2710 : !! approximation of Heaviside step function by a Chebyshev expansion plus
2711 : !! Jackson damping to reduce oscillations, applied on a set of vectors.
2712 : !! Amplifies interval [lambda_minus, lambda_plus).
2713 : !!
2714 : !! INPUTS
2715 : !! chebfi <type(chebfi_t)>=memory workspace used to apply filter
2716 : !! eigen= Rayleigh quotients to use in amplification of chebfi%xXColsRows
2717 : !! lambda_minus= lower bound of interval to amplify
2718 : !! lambda_plus= upper bound of interval to amplify
2719 : !! mineig_global= used to rescale to [-1,1), will le -1
2720 : !! maxeig_global= used to rescale to [-1,1), will be 1
2721 : !! getAX_BX= pointer to the function giving A|X> and B|X>
2722 : !! A is typically the Hamiltonian H, and B the overlap operator S
2723 : !! getBm1X= pointer to the function giving B^-1|X>
2724 : !! B is typically the overlap operator S
2725 : !!
2726 : !! SIDE EFFECTS
2727 : !! chebfi%xXColsRows= Filtered vectors to use in Subspace iteration
2728 : !!
2729 : !! SOURCE
2730 :
2731 0 : subroutine chebfi_bandpassFilterActive(chebfi,X_active,AX_active,BX_active,&
2732 : lambda_minus,lambda_plus,mineig_global,maxeig_global,getAX_BX,getBm1X)
2733 :
2734 : implicit none
2735 :
2736 : ! Arguments ------------------------------------
2737 : type(chebfi_t), intent(inout) :: chebfi
2738 : type(xgBlock_t), intent(inout) :: X_active
2739 : type(xgBlock_t), intent(inout) :: AX_active
2740 : type(xgBlock_t), intent(inout) :: BX_active
2741 : !type(xgBlock_t), intent(inout) :: eigen
2742 : real(dp), intent(in) :: lambda_minus
2743 : real(dp), intent(in) :: lambda_plus
2744 : real(dp), intent(in) :: mineig_global
2745 : real(dp), intent(in) :: maxeig_global
2746 : interface
2747 : subroutine getAX_BX(X,AX,BX)
2748 : use m_xg, only : xgBlock_t
2749 : type(xgBlock_t), intent(inout) :: X
2750 : type(xgBlock_t), intent(inout) :: AX
2751 : type(xgBlock_t), intent(inout) :: BX
2752 : end subroutine getAX_BX
2753 : end interface
2754 : interface
2755 : subroutine getBm1X(X,Bm1X)
2756 : use m_xg, only : xgBlock_t
2757 : type(xgBlock_t), intent(inout) :: X
2758 : type(xgBlock_t), intent(inout) :: Bm1X
2759 : end subroutine getBm1X
2760 : end interface
2761 :
2762 : ! Local variables-------------------------------
2763 : integer :: ndeg, n
2764 : real(dp) :: center, radius, one_over_r, two_over_r
2765 : real(dp) :: ls, us, cdeg, mu, damp
2766 : type(xg_t) :: Heaviside
2767 : type(xg_t) :: DivResults ! Rayleigh quotients
2768 : real(dp) :: tsec(2)
2769 :
2770 : ! *********************************************************************
2771 :
2772 0 : ndeg = chebfi%ndeg_filter
2773 :
2774 : ! Allocate memory for Chebyshev expansion of Heaviside step function
2775 : call xg_init(Heaviside, chebfi%space, chebfi%total_spacedim, chebfi%bandpp, chebfi%spacecom, &
2776 0 : gpu_option=chebfi%gpu_option)
2777 : !call xgBlock_zero(Heaviside%self)
2778 :
2779 0 : write(std_out,*) '@bandpass lambda_minus=', lambda_minus
2780 0 : write(std_out,*) '@bandpass lambda_plus=', lambda_plus
2781 0 : write(std_out,*) '@bandpass mineig_global=', mineig_global
2782 0 : write(std_out,*) '@bandpass maxeig_global=', maxeig_global
2783 0 : flush(std_out)
2784 :
2785 : ! Filter parameters
2786 0 : center = (maxeig_global + mineig_global)*0.5
2787 0 : radius = (maxeig_global - mineig_global)*0.5
2788 0 : one_over_r = 1/radius
2789 0 : two_over_r = 2/radius
2790 :
2791 : ! Scaled slice bounds to be amplified
2792 0 : ls = (lambda_minus - center) / radius
2793 0 : us = (lambda_plus - center) / radius
2794 :
2795 : ! TODO IL 10/3/2025 Deflate vectors to reduce linear dependence: Y=X-(B-projection)
2796 : !call xg_Borthonormalize(chebfi%xXCols,chebfi%xBxColsRows,ierr,1,chebfi%gpu_option,AX=chebfi%xAXColsRows)
2797 :
2798 : ! Initialize expansion: Heaviside = mu(0)*damp(0)*X + Heaviside
2799 0 : cdeg = Pi/(ndeg+2)
2800 0 : mu = 1/Pi*(ACOS(ls)-ACOS(us))
2801 0 : damp = 1.d0 ! Jackson damping
2802 0 : call xgBlock_saxpy(Heaviside%self, mu*damp, X_active)
2803 :
2804 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
2805 0 : do n = 0, ndeg - 1
2806 :
2807 : ! X_next=2/r*(AX_next-c*X_next)-X_prev
2808 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
2809 0 : call chebfi_computeNextOrderChebfiPolynom(chebfi, n, center, one_over_r, two_over_r, getBm1X)
2810 : ABI_NVTX_END_RANGE()
2811 :
2812 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
2813 0 : if (chebfi%paral_kgb == 0) then
2814 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%spacedim, chebfi%neigenpairs)
2815 : else
2816 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
2817 : end if
2818 : ABI_NVTX_END_RANGE()
2819 :
2820 : ! Update expansion: Heaviside = damp(i+1)*mu(i+1)*X_next + Heaviside
2821 0 : mu = 2/Pi * (SIN((n+1)*ACOS(ls)) - SIN((n+1)*ACOS(us)))/(n+1)
2822 0 : damp = ((1 - (n+1)/(ndeg+2))*SIN(cdeg)*COS((n+1)*cdeg) + 1/(ndeg+2)*COS(cdeg)*SIN((n+1)*cdeg))/SIN(cdeg)
2823 0 : call xgBlock_saxpy(Heaviside%self, mu*damp, X_active)
2824 :
2825 : ! Store final expansion before exit, X_next=Heaviside
2826 0 : if (n==ndeg-1) then
2827 0 : call xgBlock_copy(Heaviside%self, X_active)
2828 : end if
2829 :
2830 : !A * Psi (=AX_next=A*X_next)
2831 0 : call timab(tim_getAX_BX,1,tsec)
2832 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2833 0 : call getAX_BX(X_active, AX_active, BX_active)
2834 0 : call xgBlock_zero_im_g0(AX_active)
2835 0 : call xgBlock_zero_im_g0(BX_active)
2836 : ABI_NVTX_END_RANGE()
2837 0 : call timab(tim_getAX_BX,2,tsec)
2838 :
2839 : end do ! end n
2840 : ABI_NVTX_END_RANGE()
2841 :
2842 : ! Avoid overflow
2843 : !call chebfi_prepAmpfactor(chebfi, eigen, DivResults)
2844 : !call chebfi_ampfactorBandpass(chebfi, DivResults%self, lambda_minus, lambda_plus, center, radius, ndeg)
2845 :
2846 : ! Free temporary memory
2847 0 : call xg_free(DivResults)
2848 0 : call xg_free(Heaviside)
2849 :
2850 0 : end subroutine chebfi_bandpassFilterActive
2851 : !!***
2852 :
2853 : !----------------------------------------------------------------------
2854 :
2855 : !!****f* m_chebfi2/chebfi_oracle1
2856 : !! NAME
2857 : !! chebfi_oracle1
2858 : !!
2859 : !! FUNCTION
2860 : !! Compute order of Chebyshev polynom necessary to converge to a given tol
2861 : !!
2862 : !! INPUTS
2863 : !! xx= input variable
2864 : !! aa= left bound of the interval
2865 : !! bb= right bound of the interval
2866 : !! tol= needed precision
2867 : !! nmax= max number of iterations
2868 : !!
2869 : !! OUTPUT
2870 : !!
2871 : !! SIDE EFFECTS
2872 : !!
2873 : !! SOURCE
2874 :
2875 16276 : function cheb_oracle1(xx,aa,bb,tol,nmax) result(nn)
2876 :
2877 : ! Arguments ------------------------------------
2878 : integer :: nn
2879 : integer, intent(in) :: nmax
2880 : real(dp), intent(in) :: xx,aa,bb
2881 : real(dp), intent(in) :: tol
2882 :
2883 : ! Local variables-------------------------------
2884 : integer :: ii
2885 : real(dp) :: yy,yim1,xred,temp
2886 :
2887 : ! *************************************************************************
2888 :
2889 16276 : xred = (xx-(aa+bb)/2)/(bb-aa)*2
2890 16276 : yy = xred
2891 16276 : yim1 = 1 !ONE
2892 :
2893 16276 : nn = nmax
2894 16276 : if(1/(yy**2) < tol) then
2895 : nn = 1
2896 : else
2897 294686 : do ii=2, nmax-1
2898 294362 : temp = yy
2899 294362 : yy = 2*xred*yy - yim1
2900 294362 : yim1 = temp
2901 294686 : if(1/(yy**2) < tol) then
2902 : nn = ii
2903 : exit
2904 : end if
2905 : end do
2906 : end if
2907 :
2908 16276 : end function cheb_oracle1
2909 : !!***
2910 :
2911 : !----------------------------------------------------------------------
2912 :
2913 : !!****f* m_chebfi2/chebfi_poly1
2914 : !! NAME
2915 : !! chebfi_poly1
2916 : !!
2917 : !! FUNCTION
2918 : !! Compute Chebyshev polynomial???
2919 : !!
2920 : !! INPUTS
2921 : !! xx= input variable
2922 : !! aa= left bound of the interval
2923 : !! bb= right bound of the interval
2924 : !! nn=
2925 : !!
2926 : !! OUTPUT
2927 : !!
2928 : !! SIDE EFFECTS
2929 : !!
2930 : !! SOURCE
2931 :
2932 44496 : function cheb_poly1(xx,nn,aa,bb) result(yy)
2933 :
2934 : ! Arguments ------------------------------------
2935 : integer, intent(in) :: nn
2936 : real(dp), intent(in) :: xx, aa, bb
2937 : real(dp) :: yy
2938 :
2939 : ! Local variables-------------------------------
2940 : integer :: ii
2941 : real(dp) :: xred,yim1,temp
2942 :
2943 : ! *************************************************************************
2944 :
2945 44496 : xred = (xx-(aa+bb)/2)/(bb-aa)*2
2946 44496 : yy = xred
2947 44496 : yim1 = 1
2948 317168 : do ii= 2, nn
2949 272672 : temp = yy
2950 272672 : yy = 2*xred*yy - yim1
2951 317168 : yim1 = temp
2952 : end do
2953 :
2954 44496 : end function cheb_poly1
2955 : !!***
2956 :
2957 : !----------------------------------------------------------------------
2958 :
2959 : !!****f* m_chebfi2/chebfi_set_ndeg_from_residu
2960 : !! NAME
2961 : !! chebfi_set_ndeg_from_residu
2962 : !!
2963 : !! FUNCTION
2964 : !! Compute ndeg_filter using the oracle and residuals.
2965 : !!
2966 : !! INPUTS
2967 : !!
2968 : !! OUTPUT
2969 : !!
2970 : !! SIDE EFFECTS
2971 : !!
2972 : !! SOURCE
2973 :
2974 2312 : subroutine chebfi_set_ndeg_from_residu(chebfi,lambda_minus,lambda_plus,occ,DivResults,ndeg_filter_max,ndeg_filter)
2975 :
2976 : integer,intent(in) :: ndeg_filter_max
2977 : integer,intent(out) :: ndeg_filter
2978 : type(chebfi_t), intent(inout) :: chebfi
2979 : type(xgBlock_t), intent(in) :: occ
2980 : type(xgBlock_t), intent(in) :: DivResults
2981 : real(dp), intent(in) :: lambda_minus, lambda_plus
2982 :
2983 : logical :: test1,test2,test3
2984 : integer :: iband_tot,iband
2985 : integer :: bandpp,ierr,ndeg_filter_tolwfr,ndeg_filter_decrease,nbdbuf,ndeg_filter_all,shift
2986 2312 : integer,allocatable :: ndeg_filter_bands(:)
2987 : type(xgBlock_t) :: occBlock,occ_reshaped
2988 : type(xg_t) :: residu
2989 2312 : real(dp),pointer :: residu_(:,:),occ_(:,:)
2990 : real(dp) :: eig_iband,res_iband,occ_iband
2991 2312 : real(dp),pointer :: eig(:,:)
2992 : ! character(len=500) :: msg
2993 :
2994 2312 : bandpp = chebfi%bandpp
2995 :
2996 : !Compute residu here for oracle, use X_next as a work space
2997 : ! X_next = S|Psi>
2998 2312 : call xgBlock_copy(chebfi%xBXColsRows,chebfi%X_next)
2999 : ! X_next = - eig * S|Psi>
3000 2312 : call xgBlock_ymax(chebfi%X_next,DivResults,0,1)
3001 : ! X_next = H|Psi> - eig * S|Psi>
3002 2312 : call xgBlock_add(chebfi%X_next,chebfi%xAXColsRows)
3003 : ! resid = |X_next|^2
3004 2312 : call xg_init(residu,SPACE_R,bandpp,1)
3005 2312 : call xgBlock_colwiseNorm2(chebfi%X_next, residu%self,comm_loc=xmpi_comm_null)
3006 :
3007 2312 : occ_reshaped = occ
3008 2312 : shift=xmpi_comm_rank(chebfi%comm_cols)*bandpp
3009 2312 : call xgBlock_reshape(occ_reshaped,1,chebfi%neigenpairs)
3010 2312 : call xgBlock_setBlock(occ_reshaped,occBlock,1,bandpp,fcol=1+shift)
3011 2312 : call xgBlock_reshape(occBlock,bandpp,1)
3012 2312 : if (chebfi%nbdbuf==-101) then
3013 952 : call xgBlock_apply_diag(residu%self,occBlock,1)
3014 : end if
3015 :
3016 6936 : ABI_MALLOC(ndeg_filter_bands,(bandpp))
3017 :
3018 : ! DivResults could be complex (with null imaginary part), so bandpp has to be in cols, not rows
3019 2312 : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=bandpp)
3020 2312 : call xgBlock_reverseMap(residu%self,residu_,rows=1,cols=bandpp)
3021 2312 : call xgBlock_reverseMap(occBlock,occ_,rows=1,cols=bandpp)
3022 :
3023 2312 : if (chebfi%nbdbuf>0) then
3024 : nbdbuf = chebfi%nbdbuf
3025 952 : else if (chebfi%nbdbuf==-101) then
3026 952 : nbdbuf = 0
3027 : end if
3028 :
3029 12104 : do iband=1, bandpp
3030 9792 : eig_iband = eig(1,iband)
3031 9792 : res_iband = residu_(1,iband)
3032 9792 : occ_iband = occ_(1,iband)
3033 9792 : iband_tot = iband + shift
3034 9792 : test1 = res_iband<chebfi%tolerance ! band already converged
3035 9792 : test2 = iband_tot>chebfi%neigenpairs-nbdbuf ! band in the buffer
3036 9792 : test3 = chebfi%nbdbuf==-101.and.occ_iband<chebfi%oracle_min_occ ! occupancy is too low
3037 12104 : if (test1.or.test2.or.test3) then
3038 5750 : ndeg_filter_bands(iband) = 0
3039 : else
3040 : !ndeg_filter necessary to converge to tolerance
3041 4042 : ndeg_filter_tolwfr = cheb_oracle1(eig_iband, lambda_minus, lambda_plus, chebfi%tolerance / res_iband, 1000)
3042 4042 : if (chebfi%oracle==1) then
3043 832 : ndeg_filter_bands(iband) = MIN(ndeg_filter_max, ndeg_filter_tolwfr, chebfi%ndeg_filter)
3044 3210 : else if (chebfi%oracle==2) then
3045 : !ndeg_filter necessary to decrease residual by a constant factor
3046 3210 : ndeg_filter_decrease = cheb_oracle1(eig_iband, lambda_minus, lambda_plus, chebfi%oracle_factor, 15)
3047 3210 : ndeg_filter_bands(iband) = MIN(ndeg_filter_max, ndeg_filter_tolwfr, ndeg_filter_decrease)
3048 : else
3049 0 : ABI_ERROR('Wrong value for chebfi%oracle')
3050 : end if
3051 : end if
3052 : end do
3053 12104 : ndeg_filter = MAXVAL(ndeg_filter_bands)
3054 2312 : call xmpi_max(ndeg_filter,ndeg_filter_all,chebfi%comm_cols,ierr)
3055 2312 : ndeg_filter=ndeg_filter_all
3056 :
3057 2312 : call xg_free(residu)
3058 2312 : ABI_SFREE(ndeg_filter_bands)
3059 :
3060 4624 : end subroutine chebfi_set_ndeg_from_residu
3061 : !!***
3062 :
3063 : !----------------------------------------------------------------------
3064 :
3065 : !!****f* m_chebfi2/chebfi_getSubspaceResidual
3066 : !! NAME
3067 : !! chebfi_getSubspaceResidual
3068 : !!
3069 : !! FUNCTION
3070 : !! || (I-QQ'B)Aq_j ||^2 for j=1,..,bandpp
3071 : !! stored to resid
3072 :
3073 0 : subroutine chebfi_getSubspaceResidual(chebfi, resid)
3074 :
3075 : implicit none
3076 :
3077 : type(chebfi_t), intent(inout) :: chebfi
3078 : type(xgBlock_t), intent(inout) :: resid
3079 :
3080 : type(xg_t) :: M, R
3081 : integer :: nrows, ncols
3082 : integer :: space_buf
3083 :
3084 : ! *********************************************************************
3085 :
3086 0 : nrows = rows(chebfi%X) ! todo this should be chebfi%spacedim
3087 0 : ncols = chebfi%neigenpairs ! todo active part only
3088 :
3089 0 : if (chebfi%space==SPACE_C) then
3090 0 : space_buf = SPACE_C
3091 0 : else if (chebfi%space==SPACE_CR) then
3092 0 : space_buf = SPACE_R
3093 : else
3094 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
3095 : end if
3096 :
3097 0 : call xg_init(M, space_buf, ncols, ncols, comm=chebfi%spacecom, gpu_option=chebfi%gpu_option)
3098 0 : call xg_init(R, space_buf, nrows, ncols, comm=chebfi%spacecom, gpu_option=chebfi%gpu_option)
3099 :
3100 : ! Compute M = X^T*AX
3101 : ! sum all process contribution
3102 0 : call xgBlock_gemm('t','n',1.0d0,chebfi%X,chebfi%AX%self,0.d0,M%self,comm=chebfi%spacecom)
3103 :
3104 : ! Compute R = AX - BX*M
3105 0 : call xgBlock_copy(chebfi%AX%self, R%self)
3106 0 : call xgBlock_gemm('n','n',-1.0d0,chebfi%BX%self,M%self,1.0d0,R%self)
3107 :
3108 : ! sum along rows communicator
3109 0 : call xgBlock_colwiseNorm2(R%self, resid)
3110 :
3111 0 : call xg_free(M)
3112 0 : call xg_free(R)
3113 :
3114 0 : end subroutine chebfi_getSubspaceResidual
3115 : !!***
3116 :
3117 : !----------------------------------------------------------------------
3118 :
3119 : !!****f* m_chebfi2/chebfi_swapConvergedVectors
3120 : !! NAME
3121 : !! chebfi_swapConvergedVectors
3122 : !!
3123 : !! FUNCTION
3124 : !! In-place contiguous-in-memory locking of column vectors in row distribution.
3125 : !! Algorithm is based on swapping: at the end first k are locked last are active.
3126 : !! | locked columns | active columns |
3127 : !! | 1 ... k | k+1 ... m |
3128 : !! TODO logic could also be applied to divide Xext to slices.
3129 : !!
3130 :
3131 0 : subroutine chebfi_swapConvergedVectors(chebfi, resid, tol, n_locked)
3132 :
3133 : implicit none
3134 :
3135 : type(chebfi_t), intent(inout) :: chebfi
3136 : type(xgBlock_t), intent(in) :: resid
3137 : real(dp), intent(in) :: tol
3138 : integer, intent(out) :: n_locked
3139 :
3140 : integer :: i, nrows
3141 : integer :: j, left
3142 0 : real(dp), pointer :: resid_vals(:,:)
3143 0 : logical, allocatable :: mask(:)
3144 0 : logical, allocatable :: is_locked(:)
3145 0 : integer, allocatable :: idx(:)
3146 : integer :: space_buf
3147 : type(xg_t) :: T_swap ! temporary column buffer
3148 :
3149 : ! *********************************************************************
3150 :
3151 0 : nrows = rows(resid)
3152 :
3153 0 : if (chebfi%space==SPACE_C) then
3154 0 : space_buf = SPACE_C
3155 0 : else if (chebfi%space==SPACE_CR) then
3156 0 : space_buf = SPACE_R
3157 : else
3158 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
3159 : end if
3160 :
3161 0 : ABI_MALLOC(mask, (nrows))
3162 0 : ABI_MALLOC(is_locked, (nrows))
3163 :
3164 0 : call xg_init(T_swap, space_buf, nrows, 1, comm=chebfi%spacecom, gpu_option=chebfi%gpu_option)
3165 :
3166 0 : call xgBlock_reverseMap(resid, resid_vals, rows=nrows, cols=1)
3167 :
3168 0 : mask = resid_vals(:,1) < tol
3169 0 : idx = pack([(i, i=1,nrows)], mask)
3170 0 : n_locked = size(idx)
3171 :
3172 0 : write(std_out,*) 'n_locked=', n_locked
3173 0 : write(std_out,*) 'locked indices=', idx(:)
3174 0 : write(std_out,*) 'residuals=', resid_vals(:,1)
3175 0 : flush(std_out)
3176 :
3177 0 : is_locked = .false.
3178 0 : is_locked(idx) = .true.
3179 :
3180 : ! partition via swapping (in-place)
3181 0 : left = 1
3182 0 : do j = 1, nrows
3183 0 : if (is_locked(j)) then
3184 0 : if (j /= left) then
3185 0 : call xgBlock_colwiseSwap(chebfi%X, j, left, T_swap%self)
3186 0 : call xgBlock_colwiseSwap(chebfi%AX%self, j, left, T_swap%self)
3187 0 : call xgBlock_colwiseSwap(chebfi%BX%self, j, left, T_swap%self)
3188 :
3189 : ! keep mask consistent after swap
3190 0 : is_locked(j) = is_locked(left)
3191 0 : is_locked(left) = .true.
3192 : end if
3193 0 : left = left + 1
3194 : end if
3195 : end do
3196 :
3197 0 : ABI_FREE(mask)
3198 0 : ABI_FREE(is_locked)
3199 0 : call xg_free(T_swap)
3200 :
3201 0 : end subroutine chebfi_swapConvergedVectors
3202 : !!***
3203 :
3204 : !----------------------------------------------------------------------
3205 :
3206 : !!****f* m_chebfi2/chebfi_deflateWrtLocked
3207 : !! NAME
3208 : !! chebfi_deflateWrtLocked
3209 : !!
3210 :
3211 0 : subroutine chebfi_deflateWrtLocked(chebfi, n_locked, X_locked, BX_locked)
3212 :
3213 : type(chebfi_t), intent(inout) :: chebfi
3214 : integer, intent(in) :: n_locked
3215 : type(xgBlock_t), intent(inout) :: X_locked
3216 : type(xgBlock_t), intent(inout) :: BX_locked
3217 :
3218 : integer :: nrows
3219 : integer :: space_buf
3220 : integer :: n_active
3221 : type(xg_t) :: buffer
3222 : type(xgBlock_t) :: X_active, AX_active, BX_active
3223 :
3224 : ! *********************************************************************
3225 :
3226 0 : if (chebfi%space==SPACE_C) then
3227 0 : space_buf = SPACE_C
3228 0 : else if (chebfi%space==SPACE_CR) then
3229 0 : space_buf = SPACE_R
3230 : else
3231 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
3232 : end if
3233 :
3234 0 : nrows = rows(chebfi%X)
3235 0 : n_active = chebfi%neigenpairs - n_locked
3236 0 : call xg_init(buffer,space_buf,n_locked,n_active,comm=chebfi%spacecom,gpu_option=chebfi%gpu_option)
3237 :
3238 0 : write(std_out,*) 'start gemm with sum over comm'; flush(std_out)
3239 0 : write(std_out,*) 'getid', xgBlock_getid(buffer%self); flush(std_out)
3240 0 : write(std_out,*) 'getid', xgBlock_getid(buffer%self,comm=xmpi_comm_null); flush(std_out)
3241 0 : write(std_out,*) 'getid', xgBlock_getid(buffer%self,comm=chebfi%spacecom); flush(std_out)
3242 :
3243 : ! buffer = BX0^T*X
3244 : ! sum all process contribution
3245 0 : call xgBlock_setBlock(chebfi%X, X_active, nrows, n_active, fcol=n_locked+1)
3246 0 : call xgBlock_gemm('t','n',1.0d0,BX_locked,X_active,0.d0,buffer%self, comm=chebfi%spacecom)
3247 :
3248 0 : write(std_out,*) 'exit gemm with sum over comm'; flush(std_out)
3249 :
3250 0 : write(std_out,*) '@deflate W', rows(buffer%self), cols(buffer%self)
3251 0 : write(std_out,*) '@deflate A', rows(BX_locked), cols(BX_locked)
3252 0 : write(std_out,*) '@deflate B', rows(X_active), cols(X_active)
3253 0 : flush(std_out)
3254 :
3255 : ! X = X - X0*buffer
3256 0 : call xgBlock_gemm('n','n',-1.0d0,X_locked,buffer%self,1.0d0,X_active)
3257 : ! AX = AX - X0*buffer
3258 0 : call xgBlock_setBlock(chebfi%AX%self, AX_active, nrows, n_active, fcol=n_locked+1)
3259 0 : call xgBlock_gemm('n','n',-1.0d0,X_locked,buffer%self,1.0d0,AX_active)
3260 : ! BX = BX - X0*buffer
3261 0 : call xgBlock_setBlock(chebfi%BX%self, BX_active, nrows, n_active, fcol=n_locked+1)
3262 0 : call xgBlock_gemm('n','n',-1.0d0,X_locked,buffer%self,1.0d0,BX_active)
3263 :
3264 0 : call xg_free(buffer)
3265 :
3266 0 : end subroutine chebfi_deflateWrtLocked
3267 : !!***
3268 :
3269 : !----------------------------------------------------------------------
3270 :
3271 : !!****f* m_chebfi2/chebfi_deflateActiveWrtLocked
3272 : !! NAME
3273 : !! chebfi_deflateActiveWrtLocked
3274 : !!
3275 : !! FUNCTION
3276 : !! Deflate active vectors with respect to locked vectors in linalg representation
3277 :
3278 0 : subroutine chebfi_deflateActiveWrtLocked(chebfi, n_locked, X_part, AX_part, BX_part)
3279 :
3280 : type(chebfi_t), intent(inout) :: chebfi
3281 : integer, intent(in) :: n_locked
3282 : type(bandPartitionData_t), intent(inout) :: X_part
3283 : type(bandPartitionData_t), intent(inout) :: AX_part
3284 : type(bandPartitionData_t), intent(inout) :: BX_part
3285 :
3286 : integer :: space_buf
3287 : integer :: n_active
3288 : type(xg_t) :: buffer
3289 :
3290 : ! *********************************************************************
3291 :
3292 0 : if (chebfi%space==SPACE_C) then
3293 0 : space_buf = SPACE_C
3294 0 : else if (chebfi%space==SPACE_CR) then
3295 0 : space_buf = SPACE_R
3296 : else
3297 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
3298 : end if
3299 :
3300 0 : n_active = chebfi%neigenpairs - n_locked
3301 0 : call xg_init(buffer,space_buf,n_locked,n_active,comm=chebfi%spacecom,gpu_option=chebfi%gpu_option)
3302 :
3303 0 : write(std_out,*) '@deflate W', rows(buffer%self), cols(buffer%self)
3304 0 : write(std_out,*) '@deflate A', rows(BX_part%linalg_locked), cols(BX_part%linalg_locked)
3305 0 : write(std_out,*) '@deflate B', rows(X_part%linalg_active), cols(X_part%linalg_active)
3306 0 : flush(std_out)
3307 :
3308 : ! buffer = BX0^T*X
3309 : call xgBlock_gemm('t','n',1.0d0,BX_part%linalg_locked,X_part%linalg_active,0.d0,&
3310 0 : buffer%self,comm=chebfi%spacecom)
3311 :
3312 : ! sum all process contribution
3313 : ! X = X - X0*buffer
3314 0 : call xgBlock_gemm('n','n',-1.0d0,X_part%linalg_locked,buffer%self,1.0d0,X_part%linalg_active)
3315 : ! AX = AX - X0*buffer
3316 0 : call xgBlock_gemm('n','n',-1.0d0,X_part%linalg_locked,buffer%self,1.0d0,AX_part%linalg_active)
3317 : ! BX = BX - X0*buffer
3318 0 : call xgBlock_gemm('n','n',-1.0d0,X_part%linalg_locked,buffer%self,1.0d0,BX_part%linalg_active)
3319 :
3320 0 : call xg_free(buffer)
3321 :
3322 0 : end subroutine chebfi_deflateActiveWrtLocked
3323 : !!***
3324 :
3325 : !----------------------------------------------------------------------
3326 :
3327 : !!****f* m_chebfi2/bandPartitionInfo_initActiveSubcomm
3328 : !! NAME
3329 : !! bandPartitionInfo_initActiveSubcomm
3330 : !!
3331 : !! FUNCTION
3332 : !! Every process can treat between 1 and bandpp bands
3333 : !! with p processes we cover p x bandpp bands at maximum.
3334 : !! Solve problem find minimal p to cover n_active bands.
3335 : !! Constraints: use at least 1 process and at max nproc.
3336 : !!
3337 : !! INPUT
3338 : !! comm global communicator
3339 : !!
3340 : !! OUTPUT
3341 : !! mpiInfo%n_active_bandpp stores the result per mpi rank
3342 : !! mpiInfo%rank_active true if process treats active bands
3343 :
3344 0 : subroutine bandPartitionInfo_initActiveSubcomm(mpiInfo, comm, bandpp_max)
3345 :
3346 : implicit none
3347 :
3348 : type(bandPartitionInfo_t), intent(inout) :: mpiInfo
3349 : integer, intent(in) :: comm
3350 : integer, intent(in) :: bandpp_max
3351 :
3352 0 : integer, allocatable :: ranks(:)
3353 : integer :: nproc
3354 : integer :: nproc_active, min_p, ii
3355 :
3356 : ! *********************************************************************
3357 :
3358 0 : nproc = xmpi_comm_size(comm)
3359 0 : min_p = ceiling(1.d0 * mpiInfo%n_active / bandpp_max)
3360 0 : nproc_active = min(nproc, max(1, min_p))
3361 :
3362 0 : ABI_MALLOC(ranks,(0:nproc_active-1))
3363 0 : ranks(0:nproc_active-1)=(/((ii),ii=0,nproc_active-1)/)
3364 0 : mpiInfo%comm_active=xmpi_subcomm(comm,nproc_active,ranks)
3365 0 : ABI_FREE(ranks)
3366 :
3367 : ! todo @debug all processes of slice must create the same active subcomm
3368 0 : write(std_out,*) '@initActiveSubcomm comm_active=', mpiInfo%comm_active
3369 0 : flush(std_out)
3370 :
3371 :
3372 0 : end subroutine bandPartitionInfo_initActiveSubcomm
3373 : !!***
3374 :
3375 : !----------------------------------------------------------------------
3376 :
3377 : !!****f* m_chebfi2/bandPartitionData_allocActiveColsrows
3378 : !! NAME
3379 : !! bandPartitionData_allocActiveColsrows
3380 :
3381 0 : subroutine bandPartitionData_allocActiveColsrows(bpart, chebfi, mpiInfo, nspinor)
3382 :
3383 : implicit none
3384 :
3385 : type(bandPartitionData_t), intent(inout) :: bpart
3386 : type(chebfi_t), intent(in) :: chebfi
3387 : type(bandPartitionInfo_t), intent(inout) :: mpiInfo
3388 : integer, intent(in) :: nspinor
3389 :
3390 : integer :: n_active
3391 : integer :: comm_active
3392 : real(dp) :: tsec(2)
3393 :
3394 : ! *********************************************************************
3395 :
3396 0 : n_active = mpiInfo%n_active
3397 0 : comm_active = mpiInfo%comm_active
3398 :
3399 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
3400 0 : call timab(tim_transpose,1,tsec)
3401 :
3402 : call xgTransposer_constructor(bpart%transposer_active,bpart%linalg_active,bpart%colsrows_active,&
3403 : nspinor,STATE_LINALG,TRANS_ALL2ALL,chebfi%comm_rows,comm_active,0,0,chebfi%me_g0_fft,&
3404 0 : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit)
3405 :
3406 0 : bpart%transposer_active%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
3407 :
3408 : ABI_NVTX_END_RANGE()
3409 0 : call timab(tim_transpose,2,tsec)
3410 :
3411 0 : end subroutine bandPartitionData_allocActiveColsrows
3412 : !!***
3413 :
3414 : !----------------------------------------------------------------------
3415 :
3416 : !!****f* m_chebfi2/bandPartitionData_setLinalg
3417 : !! NAME
3418 : !! bandPartitionData_setLinalg
3419 : !!
3420 : !! FUNCTION
3421 : !! Set active and locked vectors in linalg representation
3422 :
3423 0 : subroutine bandPartitionData_setLinalg(X_in, bpart, n_locked)
3424 :
3425 : implicit none
3426 :
3427 : type(xgBlock_t), intent(in) :: X_in
3428 : type(bandPartitionData_t), intent(inout) :: bpart
3429 : integer, intent(in) :: n_locked
3430 :
3431 : integer :: nrows
3432 : integer :: n_active
3433 :
3434 : ! *********************************************************************
3435 :
3436 0 : nrows = rows(X_in)
3437 0 : n_active = cols(X_in) - n_locked
3438 :
3439 0 : call xgBlock_setBlock(X_in, bpart%linalg_locked, nrows, n_locked)
3440 0 : call xgBlock_setBlock(X_in, bpart%linalg_active, nrows, n_active, fcol=n_locked+1)
3441 :
3442 0 : end subroutine bandPartitionData_setLinalg
3443 : !!***
3444 :
3445 0 : end module m_chebfi2
3446 : !!***
|
