Branch data 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 : 11274 : 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 : 11274 : call timab(tim_init,1,tsec)
260 : :
261 : 11274 : chebfi%space = space
262 : 11274 : chebfi%neigenpairs = neigenpairs
263 : 11274 : chebfi%spacedim = spacedim
264 [ + + ]: 11274 : if (tolerance > 0.0) then
265 : 5242 : chebfi%tolerance = tolerance
266 : : else
267 : 6032 : chebfi%tolerance = 1.0e-20
268 : : end if
269 : 11274 : chebfi%ecut = ecut
270 : 11274 : chebfi%paral_kgb = paral_kgb
271 : 11274 : chebfi%comm_cols = comm_cols
272 : 11274 : chebfi%bandpp = bandpp
273 : 11274 : chebfi%comm_rows = comm_rows
274 : 11274 : chebfi%ndeg_filter = ndeg_filter
275 : 11274 : chebfi%nbdbuf = nbdbuf
276 : 11274 : chebfi%spacecom = spacecom
277 : 11274 : chebfi%eigenProblem = eigenProblem
278 : 11274 : chebfi%me_g0 = me_g0
279 : 11274 : chebfi%me_g0_fft = me_g0_fft
280 : 11274 : chebfi%paw = paw
281 : 11274 : chebfi%gpu_option = gpu_option
282 : 11274 : chebfi%oracle = oracle
283 : 11274 : chebfi%oracle_factor = oracle_factor
284 : 11274 : chebfi%oracle_min_occ = oracle_min_occ
285 : :
286 : 11274 : chebfi%gpu_kokkos_nthrd = 1
287 [ + - ]: 11274 : if (present(gpu_kokkos_nthrd)) chebfi%gpu_kokkos_nthrd = gpu_kokkos_nthrd
288 : 11274 : chebfi%gpu_thread_limit = 0
289 [ + - ]: 11274 : if (present(gpu_thread_limit)) chebfi%gpu_thread_limit = gpu_thread_limit
290 : 11274 : chebfi%from_linalg = .true.
291 [ + + ]: 11274 : if (present(from_linalg)) chebfi%from_linalg = from_linalg
292 : :
293 : 11274 : call chebfi_allocateAll(chebfi)
294 : :
295 : 11274 : call timab(tim_init,2,tsec)
296 : :
297 : 11274 : 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 : 11274 : 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 : 11274 : space = chebfi%space
333 : 11274 : spacedim = chebfi%spacedim
334 : 11274 : neigenpairs = chebfi%neigenpairs
335 : :
336 : 11274 : call chebfi_free(chebfi)
337 : :
338 [ + + ]: 11274 : if (chebfi%paral_kgb == 0) then
339 : 4954 : chebfi%total_spacedim = spacedim
340 : 4954 : 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 : 4954 : call xg_setBlock(chebfi%X_NP, chebfi%X_next,spacedim, neigenpairs)
342 : 4954 : 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 [ + + ]: 11274 : if (chebfi%from_linalg) then
359 : 11256 : call xg_init(chebfi%AX,space,spacedim,neigenpairs,chebfi%spacecom,me_g0=chebfi%me_g0,gpu_option=chebfi%gpu_option)
360 : 11256 : 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 : 11274 : 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 : 22548 : subroutine chebfi_free(chebfi)
393 : :
394 : : !Arguments ------------------------------------
395 : : type(chebfi_t) , intent(inout) :: chebfi
396 : :
397 : : ! *********************************************************************
398 : :
399 : 22548 : call xg_free(chebfi%X_NP)
400 : :
401 [ + + ]: 22548 : if (chebfi%from_linalg) then
402 : 22512 : call xg_free(chebfi%AX)
403 : 22512 : 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 : 22548 : 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 : 11256 : 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 : 11256 : integer,allocatable :: ndeg_filter_bands(:) !Oracle variable
582 : :
583 : : ! *********************************************************************
584 : :
585 : : ! call timab(tim_run,1,tsec)
586 : :
587 : 11256 : spacedim = chebfi%spacedim
588 : 11256 : neigenpairs = chebfi%neigenpairs
589 : 11256 : bandpp = chebfi%bandpp
590 : 11256 : ndeg_filter = chebfi%ndeg_filter
591 : 11256 : chebfi%eigenvalues = eigen
592 : :
593 : 11256 : if (chebfi%space==SPACE_C) then
594 : 7589 : space_res = SPACE_C
595 [ + - ]: 3667 : else if (chebfi%space==SPACE_CR) then
596 : 3667 : space_res = SPACE_R
597 : : else
598 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
599 : : end if
600 : :
601 [ + + ]: 11256 : if (chebfi%paral_kgb == 0) then
602 [ - + + - : 14856 : ABI_MALLOC(ndeg_filter_bands,(neigenpairs))
+ - + - -
+ - + ]
603 : 4952 : 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 : 11256 : tolerance = chebfi%tolerance
610 : 11256 : lambda_plus = chebfi%ecut
611 : 11256 : chebfi%X = X0
612 : :
613 : : ! Transpose
614 [ + + ]: 11256 : 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 : 4952 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs) !use xXColsRows instead of X notion
636 : 4952 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs) !use xAXColsRows instead of AX notion
637 : 4952 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
638 : : end if
639 : :
640 : 11256 : call timab(tim_getAX_BX,1,tsec)
641 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
642 : 11256 : call getAX_BX(chebfi%xXColsRows,chebfi%xAXColsRows,chebfi%xBXColsRows)
643 : 11256 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
644 : 11256 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
645 : : ABI_NVTX_END_RANGE()
646 : 11256 : call timab(tim_getAX_BX,2,tsec)
647 : :
648 [ + + ]: 11256 : 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 : 11256 : call timab(tim_RR_q, 1, tsec)
670 : 11256 : call chebfi_rayleighRitzQuotients(chebfi, maxeig, mineig, DivResults%self)
671 : :
672 [ + + ]: 11256 : 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 : 4952 : maxeig_global = maxeig
677 : 4952 : mineig_global = mineig
678 : : end if
679 : 11256 : call timab(tim_RR_q, 2, tsec)
680 : : ABI_NVTX_END_RANGE()
681 : :
682 : 11256 : lambda_minus = maxeig_global
683 : :
684 : 11256 : 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 : 11256 : ndeg_filter_max = cheb_oracle1(mineig_global, lambda_minus, lambda_plus, 1D-16, 40)
689 : 11256 : ndeg_filter = MIN(ndeg_filter_max,chebfi%ndeg_filter)
690 [ + + ]: 11256 : 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 [ + + ]: 74696 : ndeg_filter_bands(:) = ndeg_filter
694 : :
695 : 11256 : call timab(tim_oracle,2,tsec)
696 : :
697 : 11256 : center = (lambda_plus + lambda_minus)*0.5
698 : 11256 : radius = (lambda_plus - lambda_minus)*0.5
699 : :
700 : 11256 : one_over_r = 1/radius
701 : 11256 : two_over_r = 2/radius
702 : :
703 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
704 [ + + ]: 92674 : do ideg = 0, ndeg_filter - 1
705 : :
706 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
707 : 81418 : call chebfi_computeNextOrderChebfiPolynom(chebfi, ideg, center, one_over_r, two_over_r, getBm1X)
708 : : ABI_NVTX_END_RANGE()
709 : :
710 : 81418 : call timab(tim_swap,1,tsec)
711 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
712 [ + + ]: 81418 : if (chebfi%paral_kgb == 0) then
713 : 36906 : 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 : 81418 : call timab(tim_swap,2,tsec)
719 : :
720 : : !A * Psi
721 : 81418 : call timab(tim_getAX_BX,1,tsec)
722 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
723 : 81418 : call getAX_BX(chebfi%xXColsRows,chebfi%xAXColsRows,chebfi%xBXColsRows)
724 : 81418 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
725 : 81418 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
726 : : ABI_NVTX_END_RANGE()
727 : 92674 : call timab(tim_getAX_BX,2,tsec)
728 : :
729 : : end do ! ideg
730 : : ABI_NVTX_END_RANGE()
731 : :
732 [ + + ]: 11256 : 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 : 11256 : call timab(tim_amp_f,1,tsec)
739 : 11256 : call chebfi_ampfactor(chebfi, DivResults%self, lambda_minus, lambda_plus, ndeg_filter_bands)
740 : 11256 : call timab(tim_amp_f,2,tsec)
741 : :
742 : 11256 : call xg_free(DivResults)
743 [ + - ]: 11256 : ABI_SFREE(ndeg_filter_bands)
744 : :
745 : 11256 : call timab(tim_transpose,1,tsec)
746 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
747 [ + + ]: 11256 : 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 : 4952 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
762 : 4952 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
763 : 4952 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
764 : : end if
765 : : ABI_NVTX_END_RANGE()
766 : 11256 : 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 : 11256 : & chebfi%gpu_option,solve_ax_bx=.true.)
772 : : ABI_NVTX_END_RANGE()
773 : :
774 : : ! Compute residual
775 : 11256 : call timab(tim_residu, 1, tsec)
776 [ + + ]: 11256 : if (chebfi%paw) then
777 : 9760 : 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 : 11256 : call xgBlock_colwiseNorm2(chebfi%AX%self, residu)
783 : 11256 : call timab(tim_residu, 2, tsec)
784 : :
785 : 11256 : call timab(tim_copy, 1, tsec)
786 : 11256 : call xgBlock_copy(chebfi%X,X0)
787 : 11256 : 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 [ + + ]: 11256 : 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 [ + + - + ]: 22512 : 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 : 11256 : 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 : 11256 : if (chebfi%space==SPACE_C) then
847 : 7589 : space_res = SPACE_C
848 [ + - ]: 3667 : else if (chebfi%space==SPACE_CR) then
849 : 3667 : 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 [ + + ]: 11256 : if (chebfi%paral_kgb == 0) then
856 : 4952 : call xg_init(Results1, space_res, chebfi%neigenpairs, 1, gpu_option=chebfi%gpu_option)
857 : 4952 : 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 : 11256 : call xgBlock_colwiseDotProduct(chebfi%xXColsRows, chebfi%xAXColsRows, Results1%self, comm_loc=xmpi_comm_null)
865 : :
866 : : ! <Psi|S|Psi>
867 : 11256 : 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 : 11256 : & maxeig, maxeig_pos, mineig, mineig_pos)
872 : :
873 : 11256 : call xg_free(Results1)
874 : 11256 : call xg_free(Results2)
875 : :
876 [ + + ]: 11256 : 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 : 81598 : 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 [ + + ]: 81598 : if (chebfi%paw) then
925 : 69754 : call timab(tim_invovl, 1, tsec)
926 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_BM1X)
927 : 69754 : call getBm1X(chebfi%xAXColsRows, chebfi%X_next)
928 : : ABI_NVTX_END_RANGE()
929 : 69754 : 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 : 81598 : call timab(tim_postinvovl, 1, tsec)
937 : : ABI_NVTX_START_RANGE(NVTX_INVOVL_POST3)
938 : 81598 : call xgBlock_scale(chebfi%xXColsRows, center, 1) !scale by center
939 : :
940 : : !(B-1 * A * Psi^i-1 - c * Psi^i-1)
941 : 81598 : call xgBlock_saxpy(chebfi%X_next, dble(-1.0), chebfi%xXColsRows)
942 : :
943 : : !Psi^i-1 = 1/c * Psi^i-1
944 : 81598 : call xgBlock_scale(chebfi%xXColsRows, dble(1.0)/center, 1) !counter scale by 1/center
945 : :
946 [ + + ]: 81598 : if (ideg == 0) then
947 : 10890 : call xgBlock_scale(chebfi%X_next, one_over_r, 1)
948 : : else
949 : 70708 : call xgBlock_scale(chebfi%X_next, two_over_r, 1)
950 : :
951 : 70708 : 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 : 81598 : call timab(tim_postinvovl, 2, tsec)
961 : :
962 : 81598 : 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 : 81598 : 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 : 81598 : call xgBlock_setBlock(chebfi%X_prev, chebfi%X_swap, spacedim, neigenpairs) !X_swap = X_prev
995 : 81598 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X_prev, spacedim, neigenpairs) !X_prev = xXColsRows
996 : 81598 : call xgBlock_setBlock(chebfi%X_next, chebfi%xXColsRows, spacedim, neigenpairs) !xXColsRows = X_next
997 : 81598 : call xgBlock_setBlock(chebfi%X_swap, chebfi%X_next, spacedim, neigenpairs) !X_next = X_swap
998 : :
999 : 81598 : 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 [ - + ]: 11256 : 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 : 11256 : real(dp),pointer :: eig(:,:)
1176 : :
1177 : : ! *********************************************************************
1178 : :
1179 : 11256 : if (chebfi%paral_kgb == 0) then
1180 : 4952 : nbands = chebfi%neigenpairs
1181 : : else
1182 : 6304 : nbands = chebfi%bandpp
1183 : : end if
1184 : :
1185 : 11256 : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=chebfi%bandpp)
1186 : :
1187 [ + + ]: 74696 : do iband = 1, nbands
1188 : :
1189 : 63440 : eig_per_band = eig(1,iband)
1190 : :
1191 : : !cheb_poly1(x, n, a, b)
1192 : 63440 : ampfactor = cheb_poly1(eig_per_band, ndeg_filter_bands(iband), lambda_minus, lambda_plus)
1193 : :
1194 [ - + ]: 63440 : if(abs(ampfactor) < 1e-3) ampfactor = 1e-3 !just in case, avoid amplifying too much
1195 : :
1196 : 63440 : call xgBlock_setBlock(chebfi%xXColsRows, X_part, chebfi%total_spacedim, 1, fcol=iband)
1197 : 63440 : call xgBlock_setBlock(chebfi%xAXColsRows, AX_part, chebfi%total_spacedim, 1, fcol=iband)
1198 : 63440 : call xgBlock_setBlock(chebfi%xBXColsRows, BX_part, chebfi%total_spacedim, 1, fcol=iband)
1199 : :
1200 : 63440 : call xgBlock_scale(X_part, 1/ampfactor, 1)
1201 : 63440 : call xgBlock_scale(AX_part, 1/ampfactor, 1)
1202 : 74696 : call xgBlock_scale(BX_part, 1/ampfactor, 1)
1203 : :
1204 : : end do
1205 : :
1206 [ + + ]: 11256 : 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 : : if (.not. allocated(nrowsLinalg)) then
1387 [ # # # # : 0 : ABI_MALLOC(nrowsLinalg,(num_proc))
# # # # #
# # # ]
1388 : : endif
1389 : 0 : nrowsLinalg_ptr => nrowsLinalg
1390 [ # # # # : 0 : nrowsLinalg = nrows_blockrows
# # # # #
# # # #
# ]
1391 : :
1392 : : !write(std_out,*) 'getid inside runSlice xXColsRows', xgBlock_getId(chebfi%xXColsRows)
1393 : : !write(std_out,*) 'getid inside runSlice xAXColsRows', xgBlock_getId(chebfi%xAXColsRows)
1394 : : !flush(std_out)
1395 : :
1396 : : !A * Psi
1397 : 0 : call timab(tim_getAX_BX,1,tsec)
1398 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
1399 : 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
1400 : 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
1401 : 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
1402 : : ABI_NVTX_END_RANGE()
1403 : 0 : call timab(tim_getAX_BX,2,tsec)
1404 : :
1405 : : !write(std_out,*) 'getid inside runSlice xXColsRows (filtered 1)', xgBlock_getId(chebfi%xXColsRows)
1406 : : !write(std_out,*) 'getid inside runSlice xAXColsRows (filtered 1)', xgBlock_getId(chebfi%xAXColsRows)
1407 : :
1408 : 0 : write(std_out,*) 'starting filter in proc', xmpi_comm_rank(chebfi%spacecom)
1409 : 0 : flush(std_out)
1410 : :
1411 : : ! Apply polynomial filtering to active MPI ColsRows block-column
1412 [ # # ]: 0 : if (is_lowpass) then
1413 : : ! [lambda_minus,lambda_plus) is diminished using Chebyshev
1414 : 0 : write(std_out,*) 'lambda_minus=', lambda_minus
1415 : 0 : write(std_out,*) 'lambda_plus=', lambda_plus
1416 : 0 : flush(std_out)
1417 : 0 : call chebfi_lowpassFilter(chebfi,eigen,lambda_minus,lambda_plus,getAX_BX,getBm1X)
1418 : : else
1419 : 0 : write(std_out,*) 'lambda_minus=', lambda_minus
1420 : 0 : write(std_out,*) 'lambda_plus=', lambda_plus
1421 : 0 : write(std_out,*) 'mineig_global=', mineig_global
1422 : 0 : write(std_out,*) 'maxeig_global=', maxeig_global
1423 : 0 : flush(std_out)
1424 : : call chebfi_bandpassFilter(chebfi,lambda_minus,lambda_plus,mineig_global,&
1425 : 0 : maxeig_global,getAX_BX,getBm1X)
1426 : : end if
1427 : :
1428 : : !write(std_out,*) 'getid inside runSlice xXColsRows (filtered N)', xgBlock_getId(chebfi%xXColsRows)
1429 : : !write(std_out,*) 'getid inside runSlice xAXColsRows (filtered N)', xgBlock_getId(chebfi%xAXColsRows)
1430 : : !flush(std_out)
1431 : :
1432 : : ! MPI transpose to linalg state
1433 : 0 : call timab(tim_transpose,1,tsec)
1434 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1435 [ # # ]: 0 : if (chebfi%paral_kgb==1) then
1436 : :
1437 : : ! Allocate chebfi%X
1438 : : call xgTransposer_constructor(chebfi%xgTransposerX,chebfi%X,chebfi%xXColsRows,nspinor,&
1439 : : STATE_COLSROWS,TRANS_ALL2ALL,chebfi%comm_rows,chebfi%comm_cols,0,0,chebfi%me_g0,&
1440 : : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit,&
1441 : 0 : custom_ncolsColsRows=.true.,nrowsLinalg_sub=nrowsLinalg_ptr)
1442 : : ! Note: bandpp is custom because it is created from resource allocator
1443 : :
1444 : : ! Allocate chebfi%AX, chebfi%BX
1445 : : call xgTransposer_copyConstructor(chebfi%xgTransposerAX,chebfi%xgTransposerX,&
1446 : 0 : chebfi%AX%self,chebfi%xAXColsRows,STATE_COLSROWS)
1447 : : call xgTransposer_copyConstructor(chebfi%xgTransposerBX,chebfi%xgTransposerX,&
1448 : 0 : chebfi%BX%self,chebfi%xBXColsRows,STATE_COLSROWS)
1449 : : ! Note: at this point chebfi%AX and chebfi%BX are empty. Must transpose
1450 : : ! to fill with correct values.
1451 : :
1452 : : ! todo use copy constructor to create an object by copying an existing object
1453 : : ! actually copy constructor *allocates* memory for chebfi%AX. Write a version
1454 : : ! that does not allocate memory and only reassigns pointers. Can pointers be reassigned
1455 : : ! directly then used in the global constructor? Perform tests.
1456 : :
1457 : : !write(std_out,*) 'getid before transpose AX', xgBlock_getId(chebfi%AX%self)
1458 : : !write(std_out,*) 'getid before transpose xAX', xgBlock_getId(chebfi%xAXColsRows)
1459 : : !write(std_out,*) 'getid before transpose xX', xgBlock_getId(chebfi%xXColsRows)
1460 : : !flush(std_out)
1461 : :
1462 : 0 : chebfi%xgTransposerX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1463 : 0 : chebfi%xgTransposerAX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1464 : 0 : chebfi%xgTransposerBX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1465 : :
1466 : 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_LINALG)
1467 : 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_LINALG)
1468 : 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_LINALG)
1469 : 0 : call xmpi_barrier(chebfi%spacecom)
1470 : :
1471 : : !write(std_out,*) 'getid after transpose AX', xgBlock_getId(chebfi%AX%self)
1472 : : !flush(std_out)
1473 : :
1474 : : else
1475 : 0 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
1476 : 0 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
1477 : 0 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
1478 : : end if
1479 : 0 : call timab(tim_transpose,2,tsec)
1480 : : ABI_NVTX_END_RANGE()
1481 : :
1482 [ # # # # ]: 0 : if (chebfi%paral_kgb==1 .and. rows(chebfi%X) /= nrowsLinalg(xmpi_comm_rank(chebfi%spacecom)+1)) then
1483 : 0 : ABI_ERROR("wrong linalg representation")
1484 : : end if
1485 : 0 : write(std_out,'(a,i6,i6)') 'local # proc has # rows ', xmpi_comm_rank(chebfi%spacecom), rows(chebfi%X)
1486 : :
1487 : : !write(std_out,*) 'chebfi%eigenvalues before RR'
1488 : : !call xgBlock_print(chebfi%eigenvalues,std_out)
1489 : :
1490 : : !write(std_out,*) 'id of X, (before RR) ncols=', xgBlock_getId(chebfi%X), cols(chebfi%X)
1491 : : !flush(std_out)
1492 : :
1493 : : !call xg_Borthonormalize(chebfi%xXColsRows,chebfi%xBxColsRows,ierr,1,chebfi%gpu_option,AX=chebfi%xAXColsRows)
1494 : :
1495 : : ! Apply Rayleigh-Ritz to active MPI Linalg row-block
1496 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RR)
1497 : : call xg_RayleighRitz(chebfi%X,chebfi%AX%self,chebfi%BX%self,eigen,ierr,0,tim_RR,&
1498 : 0 : chebfi%gpu_option,solve_ax_bx=.true.)
1499 : : ABI_NVTX_END_RANGE()
1500 : :
1501 : : !write(std_out,*) 'id of X, (after RR) ncols=', xgBlock_getId(chebfi%X), cols(chebfi%X)
1502 : :
1503 [ # # ]: 0 : if ( ierr /= 0 ) then
1504 : 0 : ABI_WARNING("RayleighRitz did not work")
1505 : : else
1506 : : !write(std_out,*) 'is lowpass=', is_lowpass
1507 : : !write(std_out,*) 'chebfi%eigenvalues after RR'
1508 : : !call xgBlock_print(chebfi%eigenvalues,std_out)
1509 : : !flush(std_out)
1510 : : end if
1511 : :
1512 : : ! Compute residual norm *squared*
1513 [ # # ]: 0 : if (chebfi%paw) then
1514 : 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%BX%self,chebfi%AX%self)
1515 : : else
1516 : 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%X,chebfi%AX%self)
1517 : : end if
1518 : 0 : call xgBlock_colwiseNorm2(chebfi%AX%self,residu) ! performs MPI comm
1519 : :
1520 : : !write(std_out,*) 'max colwise residual norm squared='; call xgBlock_print(residu, std_out)
1521 : : !flush(std_out)
1522 : :
1523 : : ! Copy in Linalg representation (see chebfi_run, kept for reference)
1524 : : ! call xgBlock_copy(chebfi%X,X0)
1525 : :
1526 : : ! MPI Transpose to recover colsrows state (X only)
1527 : 0 : call timab(tim_transpose,1,tsec)
1528 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1529 [ # # ]: 0 : if (chebfi%paral_kgb == 1) then
1530 : 0 : call xmpi_barrier(chebfi%spacecom)
1531 : 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
1532 : 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
1533 : 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
1534 [ # # ]: 0 : if (xmpi_comm_size(chebfi%spacecom) == 1) then
1535 : 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
1536 : 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
1537 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
1538 : : end if
1539 : : else
1540 : 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
1541 : 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
1542 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
1543 : : end if
1544 : : ABI_NVTX_END_RANGE()
1545 : 0 : call timab(tim_transpose,2,tsec)
1546 : :
1547 : : ! Copy in ColsRows representation
1548 : 0 : call xgBlock_copy(chebfi%xXColsRows, X0)
1549 : :
1550 [ # # ]: 0 : if (cols(X0) /= chebfi%bandpp) then
1551 : 0 : ABI_ERROR('wrong colsrows representation')
1552 : : end if
1553 : 0 : write(std_out,'(a,i6,i6,i6)') 'local # proc has # rows cols ', xmpi_comm_rank(chebfi%spacecom), rows(X0), cols(X0)
1554 : :
1555 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
1556 : : if (chebfi%gpu_option==ABI_GPU_KOKKOS) then
1557 : : call gpu_device_synchronize()
1558 : : end if
1559 : : #endif
1560 : :
1561 : : ! Free transposer objects
1562 [ # # ]: 0 : if (chebfi%paral_kgb == 1) then
1563 : 0 : call xgTransposer_free(chebfi%xgTransposerX)
1564 : 0 : call xgTransposer_free(chebfi%xgTransposerAX)
1565 : 0 : call xgTransposer_free(chebfi%xgTransposerBX)
1566 : : end if
1567 : :
1568 : : ! Free temporary memory
1569 [ # # ]: 0 : ABI_SFREE(nrowsLinalg)
1570 : :
1571 : 0 : end subroutine chebfi_runSlice
1572 : : !!***
1573 : :
1574 : : !----------------------------------------------------------------------
1575 : :
1576 : : !!****f* m_chebfi/chebfi_runSubspaceIteration
1577 : : !! NAME
1578 : : !! chebfi_runSubspaceIteration
1579 : : !!
1580 : : !! FUNCTION
1581 : : !! Applies subspace iteration by Chebyshev polynomial filtering on slice
1582 : : !! then extracts eigenvectors using Rayleigh-Ritz. This execution is completely
1583 : : !! local on slice processors. Each slice does not see others.
1584 : : !!
1585 : : !! NOTES
1586 : : !! Restart logic in subspace iteration:
1587 : : !! Allocate constructor (state=colsrows)
1588 : : !! while (convergence not reached):
1589 : : !! | if (state=linalg) Transpose
1590 : : !! | Filter
1591 : : !! | Transpose (state=linalg)
1592 : : !! | Orthogonalize
1593 : : !! Rayleigh-Ritz
1594 : : !!
1595 : : !! INPUTS
1596 : : !! X0 input vectors distibuted by colsrows along slice processors
1597 : : !!
1598 : : !! OUTPUTS
1599 : : !! X0 (in-place) converged vectors distributed by linalg along slice processors
1600 : : !!
1601 : : !! dev=============
1602 : : !! - version 1: count number of converged vectors within iteration
1603 : : !!
1604 : : !! SOURCE
1605 : :
1606 : 0 : subroutine chebfi_runSubspaceIteration(chebfi,X0,getAX_BX,getBm1X,eigen,residu,nspinor,&
1607 : : mineig_global,maxeig_global,lambda_minus,lambda_plus,is_lowpass,k_conv,nrows_blockrows)
1608 : :
1609 : : implicit none
1610 : :
1611 : : !Arguments ------------------------------------
1612 : : type(chebfi_t) , intent(inout) :: chebfi
1613 : : type(xgBlock_t), intent(inout) :: X0
1614 : : type(xgBlock_t), intent(inout) :: eigen
1615 : : type(xgBlock_t), intent(inout) :: residu
1616 : : integer , intent(in ) :: nspinor
1617 : : integer , intent(in ) :: k_conv
1618 : : integer, pointer, intent(in ) :: nrows_blockrows(:)
1619 : : real(dp) , intent(in ) :: mineig_global
1620 : : real(dp) , intent(in ) :: maxeig_global
1621 : : real(dp) , intent(in ) :: lambda_minus
1622 : : real(dp) , intent(in ) :: lambda_plus
1623 : : logical , intent(in ) :: is_lowpass
1624 : : interface
1625 : : subroutine getAX_BX(X,AX,BX)
1626 : : use m_xg, only : xgBlock_t
1627 : : type(xgBlock_t), intent(inout) :: X
1628 : : type(xgBlock_t), intent(inout) :: AX
1629 : : type(xgBlock_t), intent(inout) :: BX
1630 : : end subroutine getAX_BX
1631 : : end interface
1632 : : interface
1633 : : subroutine getBm1X(X,Bm1X)
1634 : : use m_xg, only : xgBlock_t
1635 : : type(xgBlock_t), intent(inout) :: X
1636 : : type(xgBlock_t), intent(inout) :: Bm1X
1637 : : end subroutine getBm1X
1638 : : end interface
1639 : :
1640 : : !Local variables-------------------------------
1641 : : integer :: iter_subspace, niter_subspace_max, n_locked, n_locked_prev
1642 : : integer :: spacedim, neigenpairs, num_proc, nrows, nrows_mpi, ierr
1643 : : !type(xg_t) :: X_k
1644 : : type(xg_t) :: resid_active
1645 : : type(bandPartitionInfo_t) :: mpiInfo
1646 : 0 : type(bandPartitionData_t) :: X_part
1647 : 0 : type(bandPartitionData_t) :: AX_part
1648 : 0 : type(bandPartitionData_t) :: BX_part
1649 : : ! Arrays
1650 : 0 : integer, target, allocatable :: nrowsLinalg(:)
1651 : : integer, pointer :: nrowsLinalg_ptr(:) => null()
1652 : : real(dp) :: tsec(2)
1653 : :
1654 : : ! *********************************************************************
1655 : :
1656 : 0 : if (chebfi%from_linalg) then
1657 : 0 : ABI_ERROR("chebfi should be from colsrows")
1658 : : end if
1659 : :
1660 : 0 : spacedim = chebfi%spacedim
1661 : 0 : neigenpairs = chebfi%neigenpairs
1662 : 0 : num_proc = xmpi_comm_size(chebfi%spacecom)
1663 : 0 : chebfi%eigenvalues = eigen
1664 : :
1665 : : if (.not. allocated(nrowsLinalg)) then
1666 [ # # # # : 0 : ABI_MALLOC(nrowsLinalg,(num_proc))
# # # # #
# # # ]
1667 : : endif
1668 : 0 : nrowsLinalg_ptr => nrowsLinalg
1669 [ # # # # : 0 : nrowsLinalg = nrows_blockrows
# # # # #
# # # #
# ]
1670 : :
1671 : 0 : call xg_init(resid_active, SPACE_R, neigenpairs, 1, gpu_option=chebfi%gpu_option)
1672 : :
1673 : : ! Initialize values of AX and BX
1674 : 0 : call timab(tim_getAX_BX,1,tsec)
1675 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
1676 : 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
1677 : 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
1678 : 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
1679 : : ABI_NVTX_END_RANGE()
1680 : 0 : call timab(tim_getAX_BX,2,tsec)
1681 : :
1682 : : ! Allocate memory for linalg distribution from colsrows
1683 : : ! Note: bandpp is custom because it is created from resource allocator
1684 : 0 : call timab(tim_transpose,1,tsec)
1685 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1686 [ # # ]: 0 : if (chebfi%paral_kgb==1) then
1687 : :
1688 : : ! Allocate chebfi%X
1689 : : call xgTransposer_constructor(chebfi%xgTransposerX,chebfi%X,chebfi%xXColsRows,nspinor,&
1690 : : STATE_COLSROWS,TRANS_ALL2ALL,chebfi%comm_rows,chebfi%comm_cols,0,0,chebfi%me_g0,&
1691 : : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit,&
1692 : 0 : custom_ncolsColsRows=.true.,nrowsLinalg_sub=nrowsLinalg_ptr)
1693 : :
1694 : : ! Allocate chebfi%AX, chebfi%BX
1695 : : call xgTransposer_copyConstructor(chebfi%xgTransposerAX,chebfi%xgTransposerX,&
1696 : 0 : chebfi%AX%self,chebfi%xAXColsRows,STATE_COLSROWS)
1697 : : call xgTransposer_copyConstructor(chebfi%xgTransposerBX,chebfi%xgTransposerX,&
1698 : 0 : chebfi%BX%self,chebfi%xBXColsRows,STATE_COLSROWS)
1699 : : ! Note: at this point chebfi%AX and chebfi%BX are empty. Must transpose
1700 : : ! to fill with correct values.
1701 : :
1702 : 0 : chebfi%xgTransposerX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1703 : 0 : chebfi%xgTransposerAX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1704 : 0 : chebfi%xgTransposerBX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
1705 : :
1706 : : else
1707 : 0 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
1708 : 0 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
1709 : 0 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
1710 : : end if
1711 : 0 : call timab(tim_transpose,2,tsec)
1712 : : ABI_NVTX_END_RANGE()
1713 : :
1714 : 0 : write(std_out,*) 'starting filter in proc', xmpi_comm_rank(chebfi%spacecom)
1715 : 0 : flush(std_out)
1716 : :
1717 : 0 : niter_subspace_max = 10
1718 : 0 : n_locked = 0
1719 : 0 : n_locked_prev = 0
1720 : 0 : nrows = rows(chebfi%xXColsRows)
1721 : 0 : nrows_mpi = rows(chebfi%X)
1722 : :
1723 : : ! Initialize active vectors (in colsrows distribution)
1724 : 0 : call xgBlock_setBlock(chebfi%xXColsRows , X_part%colsrows_active , nrows, chebfi%bandpp)
1725 : 0 : call xgBlock_setBlock(chebfi%xAXColsRows, AX_part%colsrows_active, nrows, chebfi%bandpp)
1726 : 0 : call xgBlock_setBlock(chebfi%xBXColsRows, BX_part%colsrows_active, nrows, chebfi%bandpp)
1727 : :
1728 [ # # ]: 0 : do iter_subspace=1, niter_subspace_max
1729 : :
1730 : 0 : write(std_out,*) 'subspace iteration no=', iter_subspace
1731 : 0 : flush(std_out)
1732 : :
1733 : : ! ############################ Manage MPI ##############################
1734 [ # # # # ]: 0 : if (iter_subspace>1 .and. chebfi%paral_kgb==1) then
1735 : :
1736 [ # # ]: 0 : if (n_locked > 0 .and. n_locked /= n_locked_prev ) then
1737 : :
1738 : : ! Assign active bands to processes
1739 : 0 : call bandPartitionInfo_initActiveSubcomm(mpiInfo, chebfi%comm_cols, chebfi%bandpp)
1740 : :
1741 : : ! Allocate colsrows distributions for active vectors (from linalg distribution)
1742 : 0 : call bandPartitionData_allocActiveColsrows(X_part, chebfi, mpiInfo, nspinor)
1743 : 0 : call bandPartitionData_allocActiveColsrows(AX_part, chebfi, mpiInfo, nspinor)
1744 : 0 : call bandPartitionData_allocActiveColsrows(BX_part, chebfi, mpiInfo, nspinor)
1745 : :
1746 : : ! Transpose active vectors to colsrows distribution
1747 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1748 : 0 : call timab(tim_transpose,1,tsec)
1749 : 0 : call xgTransposer_transpose(X_part%transposer_active, STATE_COLSROWS)
1750 : 0 : call xgTransposer_transpose(AX_part%transposer_active, STATE_COLSROWS)
1751 : 0 : call xgTransposer_transpose(BX_part%transposer_active, STATE_COLSROWS)
1752 : 0 : call timab(tim_transpose,2,tsec)
1753 : : ABI_NVTX_END_RANGE()
1754 : :
1755 : : else
1756 : :
1757 : : ! Use global transposer
1758 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1759 : 0 : call timab(tim_transpose,1,tsec)
1760 : 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
1761 : 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
1762 : 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
1763 : 0 : call timab(tim_transpose,2,tsec)
1764 : : ABI_NVTX_END_RANGE()
1765 : : end if
1766 : :
1767 : : end if
1768 : :
1769 : : ! ############################ Filter active ##############################
1770 : : ! ############################ column vectors ##############################
1771 [ # # ]: 0 : if (is_lowpass) then
1772 : : call chebfi_lowpassFilterActive(chebfi,X_part%colsrows_active,AX_part%colsrows_active,&
1773 : 0 : BX_part%colsrows_active,eigen,lambda_minus,lambda_plus,getAX_BX,getBm1X)
1774 : : else
1775 : : call chebfi_bandpassFilterActive(chebfi,X_part%colsrows_active,AX_part%colsrows_active,&
1776 : : BX_part%colsrows_active,lambda_minus,lambda_plus,mineig_global,maxeig_global,&
1777 : 0 : getAX_BX,getBm1X)
1778 : : end if
1779 : : ! todo prevent bug X_next, X_prev might need resetting
1780 : :
1781 : : ! ############################ Manage MPI ##############################
1782 [ # # ]: 0 : if (chebfi%paral_kgb==1) then
1783 : :
1784 [ # # ]: 0 : if (n_locked > 0) then
1785 : :
1786 : 0 : call timab(tim_transpose,1,tsec)
1787 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1788 : 0 : call xgTransposer_transpose(X_part%transposer_active, STATE_LINALG)
1789 : 0 : call xgTransposer_transpose(AX_part%transposer_active, STATE_LINALG)
1790 : 0 : call xgTransposer_transpose(BX_part%transposer_active, STATE_LINALG)
1791 : : ABI_NVTX_END_RANGE()
1792 : 0 : call timab(tim_transpose,2,tsec)
1793 : :
1794 : : else
1795 : :
1796 : : ! Use global transposer
1797 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1798 : 0 : call timab(tim_transpose,1,tsec)
1799 : 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_LINALG)
1800 : 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_LINALG)
1801 : 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_LINALG)
1802 : 0 : call timab(tim_transpose,2,tsec)
1803 : : ABI_NVTX_END_RANGE()
1804 : :
1805 : : ! Initialize active vectors (in linalg distribution)
1806 : 0 : call xgBlock_setBlock(chebfi%X , X_part%linalg_active , nrows_mpi, neigenpairs)
1807 : 0 : call xgBlock_setBlock(chebfi%AX%self, AX_part%linalg_active, nrows_mpi, neigenpairs)
1808 : 0 : call xgBlock_setBlock(chebfi%BX%self, BX_part%linalg_active, nrows_mpi, neigenpairs)
1809 : :
1810 : : end if
1811 : :
1812 : : end if
1813 : :
1814 : 0 : write(std_out,'(a,i6,i6)') 'chebfi%X # rows # cols ', rows(chebfi%X), cols(chebfi%X)
1815 : 0 : write(std_out,'(a,i6,i6)') 'X_active(linalg) # rows # cols ', nrows_mpi, neigenpairs
1816 : 0 : flush(std_out)
1817 : :
1818 : 0 : write(std_out,*) 'getid X ', xgBlock_getid(chebfi%X)
1819 : 0 : write(std_out,*) 'getid AX ', xgBlock_getid(chebfi%AX%self)
1820 : 0 : write(std_out,*) 'getid BX ', xgBlock_getid(chebfi%BX%self)
1821 : 0 : write(std_out,*) 'space(X) ', space(chebfi%X)
1822 : 0 : write(std_out,*) 'comm(X) ', comm(chebfi%X)
1823 : 0 : write(std_out,*) 'getid Xpart ', xgBlock_getid(X_part%linalg_active)
1824 : 0 : write(std_out,*) 'getid AXpart', xgBlock_getid(AX_part%linalg_active)
1825 : 0 : write(std_out,*) 'getid BXpart', xgBlock_getid(BX_part%linalg_active)
1826 : 0 : write(std_out,*) 'space(Xpart)', space(X_part%linalg_active)
1827 : 0 : write(std_out,*) 'comm(Xpart) ', comm(X_part%linalg_active)
1828 : 0 : flush(std_out)
1829 : :
1830 : : ! ####################### Deflate+Ortho+Residual ##############################
1831 : : ! ####################### for active column vectors ##############################
1832 : :
1833 [ # # ]: 0 : if (n_locked > 0) then
1834 : 0 : call chebfi_deflateActiveWrtLocked(chebfi, n_locked, X_part, AX_part, BX_part)
1835 : : end if
1836 : :
1837 : : call xg_Borthonormalize(X_part%linalg_active,BX_part%linalg_active,ierr,1,&
1838 : 0 : chebfi%gpu_option,AX=AX_part%linalg_active)
1839 : :
1840 : 0 : call chebfi_getSubspaceResidual(chebfi, resid_active%self)
1841 : :
1842 : :
1843 : :
1844 : 0 : call chebfi_swapConvergedVectors(chebfi, resid_active%self, 1e-3_dp, n_locked)
1845 : : ! todo check if residual of wanted bands has converged if yes skip parts
1846 : : ! use k_conv as number of wanted bands
1847 : 0 : write(std_out,*) 'keeping', k_conv, 'bands'
1848 : 0 : flush(std_out)
1849 : :
1850 [ # # # # : 0 : if (chebfi%paral_kgb==1 .and. iter_subspace>1 .and. n_locked > 0 .and. n_locked /= n_locked_prev) then
# # # # ]
1851 : 0 : call bandPartitionData_setLinalg(chebfi%X , X_part , n_locked)
1852 : 0 : call bandPartitionData_setLinalg(chebfi%AX%self, AX_part, n_locked)
1853 : 0 : call bandPartitionData_setLinalg(chebfi%BX%self, BX_part, n_locked)
1854 : :
1855 : : ! free colsrows active objects
1856 : 0 : call xgTransposer_free(X_part%transposer_active)
1857 : 0 : call xgTransposer_free(AX_part%transposer_active)
1858 : 0 : call xgTransposer_free(BX_part%transposer_active)
1859 : : end if
1860 : :
1861 : 0 : n_locked_prev = n_locked
1862 : 0 : mpiInfo%n_locked = n_locked
1863 : 0 : mpiInfo%n_active = chebfi%neigenpairs - n_locked
1864 : :
1865 : : ! todo debug by recomputing the getSubspaceResidual for locked vectors and verify it is smaller than tol
1866 : :
1867 : : end do
1868 : :
1869 : : ! Transfer active data to full workspace
1870 : 0 : call xgBlock_copy(X_part%linalg_active, chebfi%X)
1871 : 0 : call xgBlock_copy(AX_part%linalg_active, chebfi%BX%self)
1872 : 0 : call xgBlock_copy(BX_part%linalg_active, chebfi%BX%self)
1873 : :
1874 : : ! Apply Rayleigh-Ritz to active MPI Linalg row-block
1875 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RR)
1876 : : call xg_RayleighRitz(chebfi%X,chebfi%AX%self,chebfi%BX%self,eigen,ierr,0,tim_RR,&
1877 : 0 : chebfi%gpu_option,solve_ax_bx=.true.)
1878 : : ABI_NVTX_END_RANGE()
1879 : :
1880 : : !write(std_out,*) 'id of X, (after RR) ncols=', xgBlock_getId(chebfi%X), cols(chebfi%X)
1881 : :
1882 [ # # ]: 0 : if ( ierr /= 0 ) then
1883 : 0 : ABI_WARNING("RayleighRitz did not work")
1884 : : end if
1885 : :
1886 : : ! Compute residual norm *squared*
1887 [ # # ]: 0 : if (chebfi%paw) then
1888 : 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%BX%self,chebfi%AX%self)
1889 : : else
1890 : 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%X,chebfi%AX%self)
1891 : : end if
1892 : 0 : call xgBlock_colwiseNorm2(chebfi%AX%self,residu) ! performs MPI comm
1893 : :
1894 : : ! MPI Transpose to recover colsrows state
1895 : 0 : call timab(tim_transpose,1,tsec)
1896 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
1897 [ # # ]: 0 : if (chebfi%paral_kgb == 1) then
1898 : 0 : call xmpi_barrier(chebfi%spacecom)
1899 : 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
1900 : 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
1901 : 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
1902 [ # # ]: 0 : if (xmpi_comm_size(chebfi%spacecom) == 1) then
1903 : 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
1904 : 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
1905 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
1906 : : end if
1907 : : else
1908 : 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
1909 : 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
1910 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
1911 : : end if
1912 : : ABI_NVTX_END_RANGE()
1913 : 0 : call timab(tim_transpose,2,tsec)
1914 : :
1915 : : ! Copy in ColsRows representation
1916 : 0 : call xgBlock_copy(chebfi%xXColsRows, X0)
1917 : :
1918 [ # # ]: 0 : if (cols(X0) /= chebfi%bandpp) then
1919 : 0 : ABI_ERROR('wrong colsrows representation')
1920 : : end if
1921 : 0 : write(std_out,'(a,i6,i6,i6)') 'local # proc has # rows cols ', xmpi_comm_rank(chebfi%spacecom), rows(X0), cols(X0)
1922 : :
1923 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
1924 : : if (chebfi%gpu_option==ABI_GPU_KOKKOS) then
1925 : : call gpu_device_synchronize()
1926 : : end if
1927 : : #endif
1928 : :
1929 : : ! Free transposer objects
1930 [ # # ]: 0 : if (chebfi%paral_kgb == 1) then
1931 : 0 : call xgTransposer_free(chebfi%xgTransposerX)
1932 : 0 : call xgTransposer_free(chebfi%xgTransposerAX)
1933 : 0 : call xgTransposer_free(chebfi%xgTransposerBX)
1934 : : end if
1935 : :
1936 : : ! Free temporary memory
1937 [ # # ]: 0 : ABI_SFREE(nrowsLinalg)
1938 : 0 : call xg_free(resid_active)
1939 : :
1940 [ # # # # : 0 : end subroutine chebfi_runSubspaceIteration
# # # # ]
1941 : : !!***
1942 : :
1943 : : !----------------------------------------------------------------------
1944 : :
1945 : : !!****f* m_chebfi/chebfi_runSubspaceIterationDummy
1946 : : !! NAME
1947 : : !! chebfi_runSubspaceIterationDummy
1948 : : !!
1949 : : !! FUNCTION
1950 : : !! Applies subspace iteration by Chebyshev polynomial filtering on slice
1951 : : !! then extracts eigenvectors using Rayleigh-Ritz. This execution is completely
1952 : : !! local on slice processors. Each slice does not see others.
1953 : : !!
1954 : : !! NOTES
1955 : : !! Restart logic in subspace iteration:
1956 : : !! Allocate constructor (state=colsrows)
1957 : : !! while (convergence not reached):
1958 : : !! | if (state=linalg) Transpose
1959 : : !! | Filter
1960 : : !! | Transpose (state=linalg)
1961 : : !! | Orthogonalize
1962 : : !! Rayleigh-Ritz
1963 : : !!
1964 : : !! INPUTS
1965 : : !! X0 input vectors distibuted by colsrows along slice processors
1966 : : !!
1967 : : !! OUTPUTS
1968 : : !! X0 (in-place) converged vectors distributed by linalg along slice processors
1969 : : !!
1970 : : !! dev=============
1971 : : !! - version 1: count number of converged vectors within iteration
1972 : : !!
1973 : : !! SOURCE
1974 : :
1975 : 0 : subroutine chebfi_runSubspaceIterationDummy(chebfi,X0,getAX_BX,getBm1X,eigen,residu,nspinor,&
1976 : : mineig_global,maxeig_global,lambda_minus,lambda_plus,is_lowpass,k_rank,nrows_blockrows)
1977 : :
1978 : : implicit none
1979 : :
1980 : : !Arguments ------------------------------------
1981 : : type(chebfi_t) , intent(inout) :: chebfi
1982 : : type(xgBlock_t), intent(inout) :: X0
1983 : : type(xgBlock_t), intent(inout) :: eigen
1984 : : type(xgBlock_t), intent(inout) :: residu
1985 : : integer , intent(in ) :: nspinor
1986 : : integer , intent(in ) :: k_rank
1987 : : integer, pointer, intent(in ) :: nrows_blockrows(:)
1988 : : real(dp) , intent(in ) :: mineig_global
1989 : : real(dp) , intent(in ) :: maxeig_global
1990 : : real(dp) , intent(in ) :: lambda_minus
1991 : : real(dp) , intent(in ) :: lambda_plus
1992 : : logical , intent(in ) :: is_lowpass
1993 : : interface
1994 : : subroutine getAX_BX(X,AX,BX)
1995 : : use m_xg, only : xgBlock_t
1996 : : type(xgBlock_t), intent(inout) :: X
1997 : : type(xgBlock_t), intent(inout) :: AX
1998 : : type(xgBlock_t), intent(inout) :: BX
1999 : : end subroutine getAX_BX
2000 : : end interface
2001 : : interface
2002 : : subroutine getBm1X(X,Bm1X)
2003 : : use m_xg, only : xgBlock_t
2004 : : type(xgBlock_t), intent(inout) :: X
2005 : : type(xgBlock_t), intent(inout) :: Bm1X
2006 : : end subroutine getBm1X
2007 : : end interface
2008 : :
2009 : : !Local variables-------------------------------
2010 : : integer :: iter_subspace, niter_subspace_max, n_locked, n_active, n_locked_prev
2011 : : integer :: spacedim, neigenpairs, num_proc, nrows, ierr
2012 : : !type(xg_t) :: X_k
2013 : : type(xg_t) :: resid_active
2014 : : type(xg_t) :: X_locked, AX_locked, BX_locked ! spaces needed for deflation in linalg distr
2015 : : type(xgBlock_t) :: chebfi_X_locked, chebfi_AX_locked, chebfi_BX_locked ! pointers used for copy
2016 : : type(xgBlock_t) :: X_active, AX_active, BX_active
2017 : : ! Arrays
2018 : 0 : integer, target, allocatable :: nrowsLinalg(:)
2019 : : integer, pointer :: nrowsLinalg_ptr(:) => null()
2020 : : real(dp) :: tsec(2)
2021 : :
2022 : : ! *********************************************************************
2023 : :
2024 : 0 : if (chebfi%from_linalg) then
2025 : 0 : ABI_ERROR("chebfi should be from colsrows")
2026 : : end if
2027 : :
2028 : 0 : spacedim = chebfi%spacedim
2029 : 0 : neigenpairs = chebfi%neigenpairs
2030 : 0 : num_proc = xmpi_comm_size(chebfi%spacecom)
2031 : 0 : chebfi%eigenvalues = eigen
2032 : :
2033 : : if (.not. allocated(nrowsLinalg)) then
2034 [ # # # # : 0 : ABI_MALLOC(nrowsLinalg,(num_proc))
# # # # #
# # # ]
2035 : : endif
2036 : 0 : nrowsLinalg_ptr => nrowsLinalg
2037 [ # # # # : 0 : nrowsLinalg = nrows_blockrows
# # # # #
# # # #
# ]
2038 : :
2039 : 0 : call xg_init(resid_active, SPACE_R, neigenpairs, 1, gpu_option=chebfi%gpu_option)
2040 : :
2041 : 0 : write(std_out,*) 'wanted rank=', k_rank
2042 : :
2043 : : !write(std_out,*) 'getid inside runSlice xXColsRows', xgBlock_getId(chebfi%xXColsRows)
2044 : : !write(std_out,*) 'getid inside runSlice xAXColsRows', xgBlock_getId(chebfi%xAXColsRows)
2045 : : !flush(std_out)
2046 : :
2047 : : ! Initialize values of AX and BX
2048 : 0 : call timab(tim_getAX_BX,1,tsec)
2049 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2050 : 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
2051 : 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
2052 : 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
2053 : : ABI_NVTX_END_RANGE()
2054 : 0 : call timab(tim_getAX_BX,2,tsec)
2055 : :
2056 : : ! Allocate memory for linalg distribution
2057 : 0 : call timab(tim_transpose,1,tsec)
2058 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
2059 [ # # ]: 0 : if (chebfi%paral_kgb==1) then
2060 : :
2061 : : ! Allocate chebfi%X
2062 : : call xgTransposer_constructor(chebfi%xgTransposerX,chebfi%X,chebfi%xXColsRows,nspinor,&
2063 : : STATE_COLSROWS,TRANS_ALL2ALL,chebfi%comm_rows,chebfi%comm_cols,0,0,chebfi%me_g0,&
2064 : : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit,&
2065 : 0 : custom_ncolsColsRows=.true.,nrowsLinalg_sub=nrowsLinalg_ptr)
2066 : : ! Note: bandpp is custom because it is created from resource allocator
2067 : :
2068 : : ! Allocate chebfi%AX, chebfi%BX
2069 : : call xgTransposer_copyConstructor(chebfi%xgTransposerAX,chebfi%xgTransposerX,&
2070 : 0 : chebfi%AX%self,chebfi%xAXColsRows,STATE_COLSROWS)
2071 : : call xgTransposer_copyConstructor(chebfi%xgTransposerBX,chebfi%xgTransposerX,&
2072 : 0 : chebfi%BX%self,chebfi%xBXColsRows,STATE_COLSROWS)
2073 : : ! Note: at this point chebfi%AX and chebfi%BX are empty. Must transpose
2074 : : ! to fill with correct values.
2075 : :
2076 : : ! todo use copy constructor to create an object by copying an existing object
2077 : : ! actually copy constructor *allocates* memory for chebfi%AX. Write a version
2078 : : ! that does not allocate memory and only reassigns pointers. Can pointers be reassigned
2079 : : ! directly then used in the global constructor? Perform tests.
2080 : :
2081 : 0 : chebfi%xgTransposerX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
2082 : 0 : chebfi%xgTransposerAX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
2083 : 0 : chebfi%xgTransposerBX%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
2084 : :
2085 : : else
2086 : 0 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X, spacedim, neigenpairs)
2087 : 0 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
2088 : 0 : call xgBlock_setBlock(chebfi%xBXColsRows, chebfi%BX%self, spacedim, neigenpairs)
2089 : : end if
2090 : 0 : call timab(tim_transpose,2,tsec)
2091 : : ABI_NVTX_END_RANGE()
2092 : :
2093 : 0 : write(std_out,*) '@dummy starting filter in proc', xmpi_comm_rank(chebfi%spacecom)
2094 : 0 : flush(std_out)
2095 : :
2096 : 0 : niter_subspace_max = 5
2097 : 0 : n_locked = 0
2098 : 0 : n_locked_prev = 0
2099 : :
2100 : : ! chebfi%X contains active vectors
2101 : : ! X_lock contains locked vectors
2102 : :
2103 : : !xX_active = chebfi%xXColsRows
2104 : : !xAX_active = chebfi%xAXColsRows
2105 : : !xBX_active = chebfi%xBXColsRows
2106 : :
2107 [ # # ]: 0 : do iter_subspace=1, niter_subspace_max
2108 : :
2109 : 0 : write(std_out,*) 'subspace iteration no=', iter_subspace
2110 : 0 : flush(std_out)
2111 : :
2112 [ # # # # ]: 0 : if (iter_subspace>1 .and. chebfi%paral_kgb==1) then
2113 : :
2114 : 0 : call timab(tim_transpose,1,tsec)
2115 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
2116 : 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
2117 : 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
2118 : 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
2119 : 0 : call timab(tim_transpose,2,tsec)
2120 : : ABI_NVTX_END_RANGE()
2121 : :
2122 : : !if (n_locked>1) then
2123 : : ! ! n_locked_mpi each rank has different number of locked vectors todo (find this in v2)
2124 : : ! call xgBlock_setBlock(chebfi%xXColsRows , xX_active , spacedim, n_locked_mpi)
2125 : : ! call xgBlock_setBlock(chebfi%xAXColsRows, xAX_active, spacedim, n_locked_mpi)
2126 : : ! call xgBlock_setBlock(chebfi%xBXColsRows, xBX_active, spacedim, n_locked_mpi)
2127 : : !end if
2128 : : end if
2129 : : ! todo this will transform the active+locked in colsrows. Locked is not necessary
2130 : : ! therefore try to communicate less data is possible.
2131 : :
2132 : : ! ############################ Filter active ##############################
2133 : : ! ############################ column vectors ##############################
2134 [ # # ]: 0 : if (is_lowpass) then
2135 : 0 : call chebfi_lowpassFilter(chebfi,eigen,lambda_minus,lambda_plus,getAX_BX,getBm1X)
2136 : : else
2137 : : call chebfi_bandpassFilter(chebfi,lambda_minus,lambda_plus,mineig_global,&
2138 : 0 : maxeig_global,getAX_BX,getBm1X)
2139 : : end if
2140 : :
2141 : 0 : write(std_out,'(a,i6,i6)') 'chebfi%xXColsRows # rows # cols ', &
2142 : 0 : rows(chebfi%xXColsRows), cols(chebfi%xXColsRows)
2143 : 0 : flush(std_out)
2144 : :
2145 [ # # ]: 0 : if (chebfi%paral_kgb==1) then
2146 : 0 : call timab(tim_transpose,1,tsec)
2147 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
2148 : 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_LINALG)
2149 : 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_LINALG)
2150 : 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_LINALG)
2151 : 0 : call timab(tim_transpose,2,tsec)
2152 : : ABI_NVTX_END_RANGE()
2153 : : end if
2154 : :
2155 [ # # ]: 0 : if (n_locked > 0) then
2156 : : ! undo the effect of filtering in locked vectors
2157 : 0 : call xgBlock_setBlock(chebfi%X, chebfi_X_locked, nrows, n_locked)
2158 : 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi_AX_locked, nrows, n_locked)
2159 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi_BX_locked, nrows, n_locked)
2160 : 0 : call xgBlock_copy(X_locked%self, chebfi_X_locked)
2161 : 0 : call xgBlock_copy(AX_locked%self, chebfi_AX_locked)
2162 : 0 : call xgBlock_copy(BX_locked%self, chebfi_BX_locked)
2163 : : end if
2164 : :
2165 : 0 : write(std_out,'(a,i6,i6)') 'chebfi%X # rows # cols ', rows(chebfi%X), cols(chebfi%X)
2166 : 0 : write(std_out,*) xgBlock_getid(chebfi%X)
2167 : 0 : flush(std_out)
2168 : :
2169 : : ! recover active vectors in linalg
2170 : 0 : n_active = chebfi%neigenpairs - n_locked
2171 : 0 : nrows = rows(chebfi%X)
2172 : :
2173 : 0 : write(std_out,*) 'n_active=', n_active; flush(std_out)
2174 : 0 : call xgBlock_setBlock(chebfi%X, X_active, nrows, n_active, fcol=n_locked+1)
2175 : 0 : call xgBlock_setBlock(chebfi%AX%self, AX_active, nrows, n_active, fcol=n_locked+1)
2176 : 0 : call xgBlock_setBlock(chebfi%BX%self, BX_active, nrows, n_active, fcol=n_locked+1)
2177 : :
2178 : 0 : write(std_out,*) 'start deflation'
2179 : 0 : write(std_out,*) xgBlock_getid(chebfi%X), comm(chebfi%X), me_g0(chebfi%X), cols(chebfi%X), space(chebfi%X)
2180 : 0 : write(std_out,*) xgBlock_getid(X_active), comm(X_active), me_g0(X_active), cols(X_active), space(X_active)
2181 : 0 : flush(std_out)
2182 : :
2183 [ # # ]: 0 : if (n_locked>0) then
2184 : 0 : call chebfi_deflateWrtLocked(chebfi, n_locked, X_locked%self, BX_locked%self)
2185 : : !call chebfi_deflateActiveWrtLocked(chebfi, n_locked, X_part, AX_part, BX_part)
2186 : : end if
2187 : :
2188 : 0 : write(std_out,*) 'exit deflation'; flush(std_out)
2189 : 0 : write(std_out,*) 'start ortho'; flush(std_out)
2190 : :
2191 : 0 : write(std_out,*) 'cols=', cols(chebfi%X); flush(std_out)
2192 : :
2193 : 0 : call xg_Borthonormalize(X_active,BX_active,ierr,1,chebfi%gpu_option,AX=AX_active)
2194 : : !call xg_Borthonormalize(chebfi%X,chebfi%BX%self,ierr,1,chebfi%gpu_option,AX=chebfi%AX%self)
2195 : :
2196 : 0 : write(std_out,*) 'exit ortho'; flush(std_out)
2197 : 0 : write(std_out,*) 'start residual'; flush(std_out)
2198 : :
2199 : 0 : call chebfi_getSubspaceResidual(chebfi, resid_active%self)
2200 : :
2201 : : ! prevent errors in pointer corruption to non-contiguous memory
2202 : 0 : call xgBlock_setBlock(chebfi%X, X_active, nrows, 1)
2203 : 0 : call xgBlock_setBlock(chebfi%AX%self, AX_active, nrows, 1)
2204 : 0 : call xgBlock_setBlock(chebfi%BX%self, BX_active, nrows, 1)
2205 : 0 : call xgBlock_setBlock(chebfi%X, chebfi_X_locked, nrows, 1)
2206 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi_BX_locked, nrows, 1)
2207 : :
2208 : 0 : call chebfi_swapConvergedVectors(chebfi, resid_active%self, 1e-3_dp, n_locked)
2209 : :
2210 [ # # ]: 0 : if (n_locked>0) then
2211 [ # # ]: 0 : if (n_locked /= n_locked_prev) then
2212 : :
2213 [ # # ]: 0 : if (n_locked_prev>0) then
2214 : 0 : call xg_free(X_locked)
2215 : 0 : call xg_free(AX_locked)
2216 : 0 : call xg_free(BX_locked)
2217 : : end if
2218 : :
2219 : : ! Allocate locked vectors in linalg representation
2220 : : call xg_init(X_locked, chebfi%space, nrows, n_locked, &
2221 : 0 : chebfi%spacecom, me_g0=chebfi%me_g0, gpu_option=chebfi%gpu_option)
2222 : :
2223 : : ! same for AX
2224 : : call xg_init(AX_locked, chebfi%space, nrows, n_locked, &
2225 : 0 : chebfi%spacecom, me_g0=chebfi%me_g0, gpu_option=chebfi%gpu_option)
2226 : :
2227 : : ! same for BX
2228 : : call xg_init(BX_locked, chebfi%space, nrows, n_locked, &
2229 : 0 : chebfi%spacecom, me_g0=chebfi%me_g0, gpu_option=chebfi%gpu_option)
2230 : : end if
2231 : :
2232 : 0 : call xgBlock_setBlock(chebfi%X, chebfi_X_locked, nrows, n_locked)
2233 : 0 : call xgBlock_copy(chebfi_X_locked, X_locked%self)
2234 : :
2235 : 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi_AX_locked, nrows, n_locked)
2236 : 0 : call xgBlock_copy(chebfi_AX_locked, AX_locked%self)
2237 : :
2238 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi_BX_locked, nrows, n_locked)
2239 : 0 : call xgBlock_copy(chebfi_BX_locked, BX_locked%self)
2240 : : end if
2241 : :
2242 : : ! todo debug by recomputing the getSubspaceResidual for locked vectors and verify it is smaller than tol
2243 : :
2244 : : ! lock vectors in small residual
2245 : : ! actually redefine pointers Xactive and Xlock pointing to X that's all
2246 : :
2247 : 0 : n_locked_prev = n_locked
2248 : :
2249 : : end do
2250 : :
2251 : 0 : call xg_free(X_locked)
2252 : 0 : call xg_free(AX_locked)
2253 : 0 : call xg_free(BX_locked)
2254 : :
2255 : : ! Apply Rayleigh-Ritz to active MPI Linalg row-block
2256 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RR)
2257 : : call xg_RayleighRitz(chebfi%X,chebfi%AX%self,chebfi%BX%self,eigen,ierr,0,tim_RR,&
2258 : 0 : chebfi%gpu_option,solve_ax_bx=.true.)
2259 : : ABI_NVTX_END_RANGE()
2260 : :
2261 : : !write(std_out,*) 'id of X, (after RR) ncols=', xgBlock_getId(chebfi%X), cols(chebfi%X)
2262 : :
2263 [ # # ]: 0 : if ( ierr /= 0 ) then
2264 : 0 : ABI_WARNING("RayleighRitz did not work")
2265 : : else
2266 : : !write(std_out,*) 'is lowpass=', is_lowpass
2267 : : !write(std_out,*) 'chebfi%eigenvalues after RR'
2268 : : !call xgBlock_print(chebfi%eigenvalues,std_out)
2269 : : !flush(std_out)
2270 : : end if
2271 : :
2272 : : ! Compute residual norm *squared*
2273 [ # # ]: 0 : if (chebfi%paw) then
2274 : 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%BX%self,chebfi%AX%self)
2275 : : else
2276 : 0 : call xgBlock_colwiseCymax(chebfi%AX%self,chebfi%eigenvalues,chebfi%X,chebfi%AX%self)
2277 : : end if
2278 : 0 : call xgBlock_colwiseNorm2(chebfi%AX%self,residu) ! performs MPI comm
2279 : :
2280 : : !write(std_out,*) 'max colwise residual norm squared='; call xgBlock_print(residu, std_out)
2281 : : !flush(std_out)
2282 : :
2283 : : ! Copy in Linalg representation (see chebfi_run, kept for reference)
2284 : : ! call xgBlock_copy(chebfi%X,X0)
2285 : :
2286 : : ! MPI Transpose to recover colsrows state (X only)
2287 : 0 : call timab(tim_transpose,1,tsec)
2288 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
2289 [ # # ]: 0 : if (chebfi%paral_kgb == 1) then
2290 : 0 : call xmpi_barrier(chebfi%spacecom)
2291 : 0 : call xgTransposer_transpose(chebfi%xgTransposerX, STATE_COLSROWS)
2292 : 0 : call xgTransposer_transpose(chebfi%xgTransposerAX, STATE_COLSROWS)
2293 : 0 : call xgTransposer_transpose(chebfi%xgTransposerBX, STATE_COLSROWS)
2294 [ # # ]: 0 : if (xmpi_comm_size(chebfi%spacecom) == 1) then
2295 : 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
2296 : 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
2297 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
2298 : : end if
2299 : : else
2300 : 0 : call xgBlock_setBlock(chebfi%X, chebfi%xXColsRows, spacedim, neigenpairs)
2301 : 0 : call xgBlock_setBlock(chebfi%AX%self, chebfi%xAXColsRows, spacedim, neigenpairs)
2302 : 0 : call xgBlock_setBlock(chebfi%BX%self, chebfi%xBXColsRows, spacedim, neigenpairs)
2303 : : end if
2304 : : ABI_NVTX_END_RANGE()
2305 : 0 : call timab(tim_transpose,2,tsec)
2306 : :
2307 : : ! Copy in ColsRows representation
2308 : 0 : call xgBlock_copy(chebfi%xXColsRows, X0)
2309 : :
2310 [ # # ]: 0 : if (cols(X0) /= chebfi%bandpp) then
2311 : 0 : ABI_ERROR('wrong colsrows representation')
2312 : : end if
2313 : 0 : write(std_out,'(a,i6,i6,i6)') 'local # proc has # rows cols ', xmpi_comm_rank(chebfi%spacecom), rows(X0), cols(X0)
2314 : :
2315 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
2316 : : if (chebfi%gpu_option==ABI_GPU_KOKKOS) then
2317 : : call gpu_device_synchronize()
2318 : : end if
2319 : : #endif
2320 : :
2321 : : ! Free transposer objects
2322 [ # # ]: 0 : if (chebfi%paral_kgb == 1) then
2323 : 0 : call xgTransposer_free(chebfi%xgTransposerX)
2324 : 0 : call xgTransposer_free(chebfi%xgTransposerAX)
2325 : 0 : call xgTransposer_free(chebfi%xgTransposerBX)
2326 : : end if
2327 : :
2328 : : ! Free temporary memory
2329 [ # # ]: 0 : ABI_SFREE(nrowsLinalg)
2330 : 0 : call xg_free(resid_active)
2331 : :
2332 [ # # ]: 0 : end subroutine chebfi_runSubspaceIterationDummy
2333 : : !!***
2334 : :
2335 : : !----------------------------------------------------------------------
2336 : :
2337 : : !!****f* m_chebfi/chebfi_lowpassFilter
2338 : : !! NAME
2339 : : !! chebfi_lowpassFilter
2340 : : !!
2341 : : !! FUNCTION
2342 : : !! Apply Lowpass filter using Chebyshev polynomial on a set of vectors.
2343 : : !! Amplifies interval [-oo, lambda_minus) and diminishes [lambda_minus,lambda_plus).
2344 : : !!
2345 : : !! INPUTS
2346 : : !! chebfi <type(chebfi_t)>=memory workspace used to apply filter
2347 : : !! eigen= Rayleigh quotients to use in amplification of chebfi%xXColsRows
2348 : : !! lambda_minus= lower bound of interval to diminish
2349 : : !! lambda_plus= upper bound of interval to diminish
2350 : : !! getAX_BX= pointer to the function giving A|X> and B|X>
2351 : : !! A is typically the Hamiltonian H, and B the overlap operator S
2352 : : !! getBm1X= pointer to the function giving B^-1|X>
2353 : : !! B is typically the overlap operator S
2354 : : !!
2355 : : !! SIDE EFFECTS
2356 : : !! chebfi%xXColsRows= Filtered vectors to use in Subspace iteration
2357 : : !!
2358 : : !! SOURCE
2359 : :
2360 : 0 : subroutine chebfi_lowpassFilter(chebfi,eigen,lambda_minus,lambda_plus,getAX_BX,getBm1X)
2361 : :
2362 : : implicit none
2363 : :
2364 : : ! Arguments ------------------------------------
2365 : : type(chebfi_t), intent(inout) :: chebfi
2366 : : type(xgBlock_t), intent(inout) :: eigen
2367 : : real(dp), intent(in) :: lambda_minus
2368 : : real(dp), intent(in) :: lambda_plus
2369 : : interface
2370 : : subroutine getAX_BX(X,AX,BX)
2371 : : use m_xg, only : xgBlock_t
2372 : : type(xgBlock_t), intent(inout) :: X
2373 : : type(xgBlock_t), intent(inout) :: AX
2374 : : type(xgBlock_t), intent(inout) :: BX
2375 : : end subroutine getAX_BX
2376 : : end interface
2377 : : interface
2378 : : subroutine getBm1X(X,Bm1X)
2379 : : use m_xg, only : xgBlock_t
2380 : : type(xgBlock_t), intent(inout) :: X
2381 : : type(xgBlock_t), intent(inout) :: Bm1X
2382 : : end subroutine getBm1X
2383 : : end interface
2384 : :
2385 : : !Local variables-------------------------------
2386 : : integer :: ideg
2387 : : real(dp) :: center, radius, one_over_r, two_over_r
2388 : : type(xg_t) :: DivResults ! Rayleigh quotients
2389 : : ! Arrays
2390 : : real(dp) :: tsec(2)
2391 : 0 : integer, allocatable :: ndeg_filter_bands(:)
2392 : :
2393 : : ! *********************************************************************
2394 : :
2395 : 0 : if (chebfi%paral_kgb == 0) then
2396 [ # # # # : 0 : ABI_MALLOC_IFNOT(ndeg_filter_bands,(chebfi%neigenpairs))
# # # # #
# # # ]
2397 : : else
2398 [ # # # # : 0 : ABI_MALLOC_IFNOT(ndeg_filter_bands,(chebfi%bandpp))
# # # # #
# # # ]
2399 : : end if
2400 : :
2401 : : ! [lambda_minus,lambda_plus) is diminished using Chebyshev
2402 : 0 : write(std_out,*) '@lowpass lambda_minus=', lambda_minus
2403 : 0 : write(std_out,*) '@lowpass lambda_plus=', lambda_plus
2404 : 0 : flush(std_out)
2405 : :
2406 : : ! Filter parameters
2407 [ # # ]: 0 : ndeg_filter_bands(:) = chebfi%ndeg_filter
2408 : 0 : center = (lambda_plus + lambda_minus)*0.5
2409 : 0 : radius = (lambda_plus - lambda_minus)*0.5
2410 : 0 : one_over_r = 1/radius
2411 : 0 : two_over_r = 2/radius
2412 : :
2413 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
2414 [ # # ]: 0 : do ideg = 0, chebfi%ndeg_filter - 1
2415 : :
2416 : : ! X_next=2/r*(AX_next-c*X_next)-X_prev
2417 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
2418 : 0 : call chebfi_computeNextOrderChebfiPolynom(chebfi, ideg, center, one_over_r, two_over_r, getBm1X)
2419 : : ABI_NVTX_END_RANGE()
2420 : :
2421 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
2422 [ # # ]: 0 : if (chebfi%paral_kgb == 0) then
2423 : 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%spacedim, chebfi%neigenpairs)
2424 : : else
2425 : 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
2426 : : end if
2427 : : ABI_NVTX_END_RANGE()
2428 : :
2429 : : !A * Psi (=AX_next=A*X_next)
2430 : 0 : call timab(tim_getAX_BX,1,tsec)
2431 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2432 : 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
2433 : 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
2434 : 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
2435 : : ABI_NVTX_END_RANGE()
2436 : 0 : call timab(tim_getAX_BX,2,tsec)
2437 : :
2438 : : end do ! ideg
2439 : : ABI_NVTX_END_RANGE()
2440 : :
2441 : : ! Avoid overflow rescale
2442 : 0 : call chebfi_prepAmpfactor(chebfi, eigen, DivResults)
2443 : 0 : call chebfi_ampfactor(chebfi, DivResults%self, lambda_minus, lambda_plus, ndeg_filter_bands)
2444 : :
2445 : : ! Free temporary memory
2446 : 0 : call xg_free(DivResults)
2447 [ # # ]: 0 : ABI_SFREE(ndeg_filter_bands)
2448 : :
2449 [ # # # # ]: 0 : end subroutine chebfi_lowpassFilter
2450 : : !!***
2451 : :
2452 : : !----------------------------------------------------------------------
2453 : :
2454 : : !!****f* m_chebfi/chebfi_lowpassFilterActive
2455 : : !! NAME
2456 : : !! chebfi_lowpassFilterActive
2457 : : !!
2458 : : !! FUNCTION
2459 : : !! Apply Lowpass filter on Active vectors only
2460 : : !! TODO this function will replace lowpassFilter
2461 : : !!
2462 : : !! SOURCE
2463 : :
2464 : 0 : subroutine chebfi_lowpassFilterActive(chebfi,X_active,AX_active,BX_active,eigen,&
2465 : : lambda_minus,lambda_plus,getAX_BX,getBm1X)
2466 : :
2467 : : implicit none
2468 : :
2469 : : ! Arguments ------------------------------------
2470 : : type(chebfi_t), intent(inout) :: chebfi
2471 : : type(xgBlock_t), intent(inout) :: X_active
2472 : : type(xgBlock_t), intent(inout) :: AX_active
2473 : : type(xgBlock_t), intent(inout) :: BX_active
2474 : : type(xgBlock_t), intent(inout) :: eigen
2475 : : real(dp), intent(in) :: lambda_minus
2476 : : real(dp), intent(in) :: lambda_plus
2477 : : interface
2478 : : subroutine getAX_BX(X,AX,BX)
2479 : : use m_xg, only : xgBlock_t
2480 : : type(xgBlock_t), intent(inout) :: X
2481 : : type(xgBlock_t), intent(inout) :: AX
2482 : : type(xgBlock_t), intent(inout) :: BX
2483 : : end subroutine getAX_BX
2484 : : end interface
2485 : : interface
2486 : : subroutine getBm1X(X,Bm1X)
2487 : : use m_xg, only : xgBlock_t
2488 : : type(xgBlock_t), intent(inout) :: X
2489 : : type(xgBlock_t), intent(inout) :: Bm1X
2490 : : end subroutine getBm1X
2491 : : end interface
2492 : :
2493 : : !Local variables-------------------------------
2494 : : integer :: ideg
2495 : : real(dp) :: center, radius, one_over_r, two_over_r
2496 : : type(xg_t) :: DivResults ! Rayleigh quotients
2497 : : ! Arrays
2498 : : real(dp) :: tsec(2)
2499 : 0 : integer, allocatable :: ndeg_filter_bands(:)
2500 : :
2501 : : ! *********************************************************************
2502 : :
2503 : 0 : if (chebfi%paral_kgb == 0) then
2504 [ # # # # : 0 : ABI_MALLOC_IFNOT(ndeg_filter_bands,(chebfi%neigenpairs))
# # # # #
# # # ]
2505 : : else
2506 [ # # # # : 0 : ABI_MALLOC_IFNOT(ndeg_filter_bands,(chebfi%bandpp))
# # # # #
# # # ]
2507 : : end if
2508 : :
2509 : : ! [lambda_minus,lambda_plus) is diminished using Chebyshev
2510 : 0 : write(std_out,*) '@lowpass lambda_minus=', lambda_minus
2511 : 0 : write(std_out,*) '@lowpass lambda_plus=', lambda_plus
2512 : 0 : flush(std_out)
2513 : :
2514 : : ! Filter parameters
2515 [ # # ]: 0 : ndeg_filter_bands(:) = chebfi%ndeg_filter
2516 : 0 : center = (lambda_plus + lambda_minus)*0.5
2517 : 0 : radius = (lambda_plus - lambda_minus)*0.5
2518 : 0 : one_over_r = 1/radius
2519 : 0 : two_over_r = 2/radius
2520 : :
2521 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
2522 [ # # ]: 0 : do ideg = 0, chebfi%ndeg_filter - 1
2523 : :
2524 : : ! X_next=2/r*(AX_next-c*X_next)-X_prev
2525 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
2526 : 0 : call chebfi_computeNextOrderChebfiPolynom(chebfi, ideg, center, one_over_r, two_over_r, getBm1X)
2527 : : ABI_NVTX_END_RANGE()
2528 : :
2529 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
2530 [ # # ]: 0 : if (chebfi%paral_kgb == 0) then
2531 : 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%spacedim, chebfi%neigenpairs)
2532 : : else
2533 : 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
2534 : : end if
2535 : : ABI_NVTX_END_RANGE()
2536 : :
2537 : : !A * Psi (=AX_next=A*X_next)
2538 : 0 : call timab(tim_getAX_BX,1,tsec)
2539 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2540 : 0 : call getAX_BX(X_active, AX_active, BX_active)
2541 : 0 : call xgBlock_zero_im_g0(AX_active)
2542 : 0 : call xgBlock_zero_im_g0(BX_active)
2543 : : ABI_NVTX_END_RANGE()
2544 : 0 : call timab(tim_getAX_BX,2,tsec)
2545 : :
2546 : : end do ! ideg
2547 : : ABI_NVTX_END_RANGE()
2548 : :
2549 : : ! Avoid overflow rescale
2550 : 0 : call chebfi_prepAmpfactor(chebfi, eigen, DivResults)
2551 : 0 : call chebfi_ampfactor(chebfi, DivResults%self, lambda_minus, lambda_plus, ndeg_filter_bands)
2552 : :
2553 : : ! Free temporary memory
2554 : 0 : call xg_free(DivResults)
2555 [ # # ]: 0 : ABI_SFREE(ndeg_filter_bands)
2556 : :
2557 [ # # # # ]: 0 : end subroutine chebfi_lowpassFilterActive
2558 : : !!***
2559 : :
2560 : : !----------------------------------------------------------------------
2561 : :
2562 : : !!****f* m_chebfi/chebfi_bandpassFilter
2563 : : !! NAME
2564 : : !! chebfi_bandpassFilter
2565 : : !!
2566 : : !! FUNCTION
2567 : : !! Apply Bandpass filter using Chebyshev-Jackson polynomial, that is an
2568 : : !! approximation of Heaviside step function by a Chebyshev expansion plus
2569 : : !! Jackson damping to reduce oscillations, applied on a set of vectors.
2570 : : !! Amplifies interval [lambda_minus, lambda_plus).
2571 : : !!
2572 : : !! INPUTS
2573 : : !! chebfi <type(chebfi_t)>=memory workspace used to apply filter
2574 : : !! eigen= Rayleigh quotients to use in amplification of chebfi%xXColsRows
2575 : : !! lambda_minus= lower bound of interval to amplify
2576 : : !! lambda_plus= upper bound of interval to amplify
2577 : : !! mineig_global= used to rescale to [-1,1), will le -1
2578 : : !! maxeig_global= used to rescale to [-1,1), will be 1
2579 : : !! getAX_BX= pointer to the function giving A|X> and B|X>
2580 : : !! A is typically the Hamiltonian H, and B the overlap operator S
2581 : : !! getBm1X= pointer to the function giving B^-1|X>
2582 : : !! B is typically the overlap operator S
2583 : : !!
2584 : : !! SIDE EFFECTS
2585 : : !! chebfi%xXColsRows= Filtered vectors to use in Subspace iteration
2586 : : !!
2587 : : !! SOURCE
2588 : :
2589 : 0 : subroutine chebfi_bandpassFilter(chebfi,lambda_minus,lambda_plus,mineig_global,&
2590 : : maxeig_global,getAX_BX,getBm1X)
2591 : :
2592 : : implicit none
2593 : :
2594 : : ! Arguments ------------------------------------
2595 : : type(chebfi_t), intent(inout) :: chebfi
2596 : : !type(xgBlock_t), intent(inout) :: eigen
2597 : : real(dp), intent(in) :: lambda_minus
2598 : : real(dp), intent(in) :: lambda_plus
2599 : : real(dp), intent(in) :: mineig_global
2600 : : real(dp), intent(in) :: maxeig_global
2601 : : interface
2602 : : subroutine getAX_BX(X,AX,BX)
2603 : : use m_xg, only : xgBlock_t
2604 : : type(xgBlock_t), intent(inout) :: X
2605 : : type(xgBlock_t), intent(inout) :: AX
2606 : : type(xgBlock_t), intent(inout) :: BX
2607 : : end subroutine getAX_BX
2608 : : end interface
2609 : : interface
2610 : : subroutine getBm1X(X,Bm1X)
2611 : : use m_xg, only : xgBlock_t
2612 : : type(xgBlock_t), intent(inout) :: X
2613 : : type(xgBlock_t), intent(inout) :: Bm1X
2614 : : end subroutine getBm1X
2615 : : end interface
2616 : :
2617 : : ! Local variables-------------------------------
2618 : : integer :: ndeg, n
2619 : : real(dp) :: center, radius, one_over_r, two_over_r
2620 : : real(dp) :: ls, us, cdeg, mu, damp
2621 : : type(xg_t) :: Heaviside
2622 : : type(xg_t) :: DivResults ! Rayleigh quotients
2623 : : real(dp) :: tsec(2)
2624 : :
2625 : : ! *********************************************************************
2626 : :
2627 : 0 : ndeg = chebfi%ndeg_filter
2628 : :
2629 : : ! Allocate memory for Chebyshev expansion of Heaviside step function
2630 : : call xg_init(Heaviside, chebfi%space, chebfi%total_spacedim, chebfi%bandpp, chebfi%spacecom, &
2631 : 0 : gpu_option=chebfi%gpu_option)
2632 : : !call xgBlock_zero(Heaviside%self)
2633 : :
2634 : 0 : write(std_out,*) '@bandpass lambda_minus=', lambda_minus
2635 : 0 : write(std_out,*) '@bandpass lambda_plus=', lambda_plus
2636 : 0 : write(std_out,*) '@bandpass mineig_global=', mineig_global
2637 : 0 : write(std_out,*) '@bandpass maxeig_global=', maxeig_global
2638 : 0 : flush(std_out)
2639 : :
2640 : : ! Filter parameters
2641 : 0 : center = (maxeig_global + mineig_global)*0.5
2642 : 0 : radius = (maxeig_global - mineig_global)*0.5
2643 : 0 : one_over_r = 1/radius
2644 : 0 : two_over_r = 2/radius
2645 : :
2646 : : ! Scaled slice bounds to be amplified
2647 : 0 : ls = (lambda_minus - center) / radius
2648 : 0 : us = (lambda_plus - center) / radius
2649 : :
2650 : : ! TODO IL 10/3/2025 Deflate vectors to reduce linear dependence: Y=X-(B-projection)
2651 : : !call xg_Borthonormalize(chebfi%xXColsRows,chebfi%xBxColsRows,ierr,1,chebfi%gpu_option,AX=chebfi%xAXColsRows)
2652 : :
2653 : : ! Initialize expansion: Heaviside = mu(0)*damp(0)*X + Heaviside
2654 : 0 : cdeg = Pi/(ndeg+2)
2655 : 0 : mu = 1/Pi*(ACOS(ls)-ACOS(us))
2656 : 0 : damp = 1.d0 ! Jackson damping
2657 : 0 : call xgBlock_saxpy(Heaviside%self, mu*damp, chebfi%xXColsRows)
2658 : :
2659 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
2660 [ # # ]: 0 : do n = 0, ndeg - 1
2661 : :
2662 : : ! X_next=2/r*(AX_next-c*X_next)-X_prev
2663 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
2664 : 0 : call chebfi_computeNextOrderChebfiPolynom(chebfi, n, center, one_over_r, two_over_r, getBm1X)
2665 : : ABI_NVTX_END_RANGE()
2666 : :
2667 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
2668 [ # # ]: 0 : if (chebfi%paral_kgb == 0) then
2669 : 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%spacedim, chebfi%neigenpairs)
2670 : : else
2671 : 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
2672 : : end if
2673 : : ABI_NVTX_END_RANGE()
2674 : :
2675 : : ! Update expansion: Heaviside = damp(i+1)*mu(i+1)*X_next + Heaviside
2676 : 0 : mu = 2/Pi * (SIN((n+1)*ACOS(ls)) - SIN((n+1)*ACOS(us)))/(n+1)
2677 : 0 : damp = ((1 - (n+1)/(ndeg+2))*SIN(cdeg)*COS((n+1)*cdeg) + 1/(ndeg+2)*COS(cdeg)*SIN((n+1)*cdeg))/SIN(cdeg)
2678 : 0 : call xgBlock_saxpy(Heaviside%self, mu*damp, chebfi%xXColsRows)
2679 : :
2680 : : ! Store final expansion before exit, X_next=Heaviside
2681 [ # # ]: 0 : if (n==ndeg-1) then
2682 : 0 : call xgBlock_copy(Heaviside%self, chebfi%xXColsRows)
2683 : : end if
2684 : :
2685 : : !A * Psi (=AX_next=A*X_next)
2686 : 0 : call timab(tim_getAX_BX,1,tsec)
2687 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2688 : 0 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
2689 : 0 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
2690 : 0 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
2691 : : ABI_NVTX_END_RANGE()
2692 : 0 : call timab(tim_getAX_BX,2,tsec)
2693 : :
2694 : : end do ! end n
2695 : : ABI_NVTX_END_RANGE()
2696 : :
2697 : : ! Avoid overflow
2698 : : !call chebfi_prepAmpfactor(chebfi, eigen, DivResults)
2699 : : !call chebfi_ampfactorBandpass(chebfi, DivResults%self, lambda_minus, lambda_plus, center, radius, ndeg)
2700 : :
2701 : : ! Free temporary memory
2702 : 0 : call xg_free(DivResults)
2703 : 0 : call xg_free(Heaviside)
2704 : :
2705 : 0 : end subroutine chebfi_bandpassFilter
2706 : : !!***
2707 : :
2708 : : !----------------------------------------------------------------------
2709 : :
2710 : : !!****f* m_chebfi/chebfi_bandpassFilterActive
2711 : : !! NAME
2712 : : !! chebfi_bandpassFilterActive
2713 : : !!
2714 : : !! FUNCTION
2715 : : !! Apply Bandpass filter using Chebyshev-Jackson polynomial, that is an
2716 : : !! approximation of Heaviside step function by a Chebyshev expansion plus
2717 : : !! Jackson damping to reduce oscillations, applied on a set of vectors.
2718 : : !! Amplifies interval [lambda_minus, lambda_plus).
2719 : : !!
2720 : : !! INPUTS
2721 : : !! chebfi <type(chebfi_t)>=memory workspace used to apply filter
2722 : : !! eigen= Rayleigh quotients to use in amplification of chebfi%xXColsRows
2723 : : !! lambda_minus= lower bound of interval to amplify
2724 : : !! lambda_plus= upper bound of interval to amplify
2725 : : !! mineig_global= used to rescale to [-1,1), will le -1
2726 : : !! maxeig_global= used to rescale to [-1,1), will be 1
2727 : : !! getAX_BX= pointer to the function giving A|X> and B|X>
2728 : : !! A is typically the Hamiltonian H, and B the overlap operator S
2729 : : !! getBm1X= pointer to the function giving B^-1|X>
2730 : : !! B is typically the overlap operator S
2731 : : !!
2732 : : !! SIDE EFFECTS
2733 : : !! chebfi%xXColsRows= Filtered vectors to use in Subspace iteration
2734 : : !!
2735 : : !! SOURCE
2736 : :
2737 : 0 : subroutine chebfi_bandpassFilterActive(chebfi,X_active,AX_active,BX_active,&
2738 : : lambda_minus,lambda_plus,mineig_global,maxeig_global,getAX_BX,getBm1X)
2739 : :
2740 : : implicit none
2741 : :
2742 : : ! Arguments ------------------------------------
2743 : : type(chebfi_t), intent(inout) :: chebfi
2744 : : type(xgBlock_t), intent(inout) :: X_active
2745 : : type(xgBlock_t), intent(inout) :: AX_active
2746 : : type(xgBlock_t), intent(inout) :: BX_active
2747 : : !type(xgBlock_t), intent(inout) :: eigen
2748 : : real(dp), intent(in) :: lambda_minus
2749 : : real(dp), intent(in) :: lambda_plus
2750 : : real(dp), intent(in) :: mineig_global
2751 : : real(dp), intent(in) :: maxeig_global
2752 : : interface
2753 : : subroutine getAX_BX(X,AX,BX)
2754 : : use m_xg, only : xgBlock_t
2755 : : type(xgBlock_t), intent(inout) :: X
2756 : : type(xgBlock_t), intent(inout) :: AX
2757 : : type(xgBlock_t), intent(inout) :: BX
2758 : : end subroutine getAX_BX
2759 : : end interface
2760 : : interface
2761 : : subroutine getBm1X(X,Bm1X)
2762 : : use m_xg, only : xgBlock_t
2763 : : type(xgBlock_t), intent(inout) :: X
2764 : : type(xgBlock_t), intent(inout) :: Bm1X
2765 : : end subroutine getBm1X
2766 : : end interface
2767 : :
2768 : : ! Local variables-------------------------------
2769 : : integer :: ndeg, n
2770 : : real(dp) :: center, radius, one_over_r, two_over_r
2771 : : real(dp) :: ls, us, cdeg, mu, damp
2772 : : type(xg_t) :: Heaviside
2773 : : type(xg_t) :: DivResults ! Rayleigh quotients
2774 : : real(dp) :: tsec(2)
2775 : :
2776 : : ! *********************************************************************
2777 : :
2778 : 0 : ndeg = chebfi%ndeg_filter
2779 : :
2780 : : ! Allocate memory for Chebyshev expansion of Heaviside step function
2781 : : call xg_init(Heaviside, chebfi%space, chebfi%total_spacedim, chebfi%bandpp, chebfi%spacecom, &
2782 : 0 : gpu_option=chebfi%gpu_option)
2783 : : !call xgBlock_zero(Heaviside%self)
2784 : :
2785 : 0 : write(std_out,*) '@bandpass lambda_minus=', lambda_minus
2786 : 0 : write(std_out,*) '@bandpass lambda_plus=', lambda_plus
2787 : 0 : write(std_out,*) '@bandpass mineig_global=', mineig_global
2788 : 0 : write(std_out,*) '@bandpass maxeig_global=', maxeig_global
2789 : 0 : flush(std_out)
2790 : :
2791 : : ! Filter parameters
2792 : 0 : center = (maxeig_global + mineig_global)*0.5
2793 : 0 : radius = (maxeig_global - mineig_global)*0.5
2794 : 0 : one_over_r = 1/radius
2795 : 0 : two_over_r = 2/radius
2796 : :
2797 : : ! Scaled slice bounds to be amplified
2798 : 0 : ls = (lambda_minus - center) / radius
2799 : 0 : us = (lambda_plus - center) / radius
2800 : :
2801 : : ! TODO IL 10/3/2025 Deflate vectors to reduce linear dependence: Y=X-(B-projection)
2802 : : !call xg_Borthonormalize(chebfi%xXCols,chebfi%xBxColsRows,ierr,1,chebfi%gpu_option,AX=chebfi%xAXColsRows)
2803 : :
2804 : : ! Initialize expansion: Heaviside = mu(0)*damp(0)*X + Heaviside
2805 : 0 : cdeg = Pi/(ndeg+2)
2806 : 0 : mu = 1/Pi*(ACOS(ls)-ACOS(us))
2807 : 0 : damp = 1.d0 ! Jackson damping
2808 : 0 : call xgBlock_saxpy(Heaviside%self, mu*damp, X_active)
2809 : :
2810 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_CORE)
2811 [ # # ]: 0 : do n = 0, ndeg - 1
2812 : :
2813 : : ! X_next=2/r*(AX_next-c*X_next)-X_prev
2814 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
2815 : 0 : call chebfi_computeNextOrderChebfiPolynom(chebfi, n, center, one_over_r, two_over_r, getBm1X)
2816 : : ABI_NVTX_END_RANGE()
2817 : :
2818 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
2819 [ # # ]: 0 : if (chebfi%paral_kgb == 0) then
2820 : 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%spacedim, chebfi%neigenpairs)
2821 : : else
2822 : 0 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
2823 : : end if
2824 : : ABI_NVTX_END_RANGE()
2825 : :
2826 : : ! Update expansion: Heaviside = damp(i+1)*mu(i+1)*X_next + Heaviside
2827 : 0 : mu = 2/Pi * (SIN((n+1)*ACOS(ls)) - SIN((n+1)*ACOS(us)))/(n+1)
2828 : 0 : damp = ((1 - (n+1)/(ndeg+2))*SIN(cdeg)*COS((n+1)*cdeg) + 1/(ndeg+2)*COS(cdeg)*SIN((n+1)*cdeg))/SIN(cdeg)
2829 : 0 : call xgBlock_saxpy(Heaviside%self, mu*damp, X_active)
2830 : :
2831 : : ! Store final expansion before exit, X_next=Heaviside
2832 [ # # ]: 0 : if (n==ndeg-1) then
2833 : 0 : call xgBlock_copy(Heaviside%self, X_active)
2834 : : end if
2835 : :
2836 : : !A * Psi (=AX_next=A*X_next)
2837 : 0 : call timab(tim_getAX_BX,1,tsec)
2838 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
2839 : 0 : call getAX_BX(X_active, AX_active, BX_active)
2840 : 0 : call xgBlock_zero_im_g0(AX_active)
2841 : 0 : call xgBlock_zero_im_g0(BX_active)
2842 : : ABI_NVTX_END_RANGE()
2843 : 0 : call timab(tim_getAX_BX,2,tsec)
2844 : :
2845 : : end do ! end n
2846 : : ABI_NVTX_END_RANGE()
2847 : :
2848 : : ! Avoid overflow
2849 : : !call chebfi_prepAmpfactor(chebfi, eigen, DivResults)
2850 : : !call chebfi_ampfactorBandpass(chebfi, DivResults%self, lambda_minus, lambda_plus, center, radius, ndeg)
2851 : :
2852 : : ! Free temporary memory
2853 : 0 : call xg_free(DivResults)
2854 : 0 : call xg_free(Heaviside)
2855 : :
2856 : 0 : end subroutine chebfi_bandpassFilterActive
2857 : : !!***
2858 : :
2859 : : !----------------------------------------------------------------------
2860 : :
2861 : : !!****f* m_chebfi2/chebfi_oracle1
2862 : : !! NAME
2863 : : !! chebfi_oracle1
2864 : : !!
2865 : : !! FUNCTION
2866 : : !! Compute order of Chebyshev polynom necessary to converge to a given tol
2867 : : !!
2868 : : !! INPUTS
2869 : : !! xx= input variable
2870 : : !! aa= left bound of the interval
2871 : : !! bb= right bound of the interval
2872 : : !! tol= needed precision
2873 : : !! nmax= max number of iterations
2874 : : !!
2875 : : !! OUTPUT
2876 : : !!
2877 : : !! SIDE EFFECTS
2878 : : !!
2879 : : !! SOURCE
2880 : :
2881 : 18508 : function cheb_oracle1(xx,aa,bb,tol,nmax) result(nn)
2882 : :
2883 : : ! Arguments ------------------------------------
2884 : : integer :: nn
2885 : : integer, intent(in) :: nmax
2886 : : real(dp), intent(in) :: xx,aa,bb
2887 : : real(dp), intent(in) :: tol
2888 : :
2889 : : ! Local variables-------------------------------
2890 : : integer :: ii
2891 : : real(dp) :: yy,yim1,xred,temp
2892 : :
2893 : : ! *************************************************************************
2894 : :
2895 : 18508 : xred = (xx-(aa+bb)/2)/(bb-aa)*2
2896 : 18508 : yy = xred
2897 : 18508 : yim1 = 1 !ONE
2898 : :
2899 : 18508 : nn = nmax
2900 [ + + ]: 18508 : if(1/(yy**2) < tol) then
2901 : : nn = 1
2902 : : else
2903 [ + + ]: 348284 : do ii=2, nmax-1
2904 : 347960 : temp = yy
2905 : 347960 : yy = 2*xred*yy - yim1
2906 : 347960 : yim1 = temp
2907 [ + + ]: 348284 : if(1/(yy**2) < tol) then
2908 : : nn = ii
2909 : : exit
2910 : : end if
2911 : : end do
2912 : : end if
2913 : :
2914 : 18508 : end function cheb_oracle1
2915 : : !!***
2916 : :
2917 : : !----------------------------------------------------------------------
2918 : :
2919 : : !!****f* m_chebfi2/chebfi_poly1
2920 : : !! NAME
2921 : : !! chebfi_poly1
2922 : : !!
2923 : : !! FUNCTION
2924 : : !! Compute Chebyshev polynomial???
2925 : : !!
2926 : : !! INPUTS
2927 : : !! xx= input variable
2928 : : !! aa= left bound of the interval
2929 : : !! bb= right bound of the interval
2930 : : !! nn=
2931 : : !!
2932 : : !! OUTPUT
2933 : : !!
2934 : : !! SIDE EFFECTS
2935 : : !!
2936 : : !! SOURCE
2937 : :
2938 : 63440 : function cheb_poly1(xx,nn,aa,bb) result(yy)
2939 : :
2940 : : ! Arguments ------------------------------------
2941 : : integer, intent(in) :: nn
2942 : : real(dp), intent(in) :: xx, aa, bb
2943 : : real(dp) :: yy
2944 : :
2945 : : ! Local variables-------------------------------
2946 : : integer :: ii
2947 : : real(dp) :: xred,yim1,temp
2948 : :
2949 : : ! *************************************************************************
2950 : :
2951 : 63440 : xred = (xx-(aa+bb)/2)/(bb-aa)*2
2952 : 63440 : yy = xred
2953 : 63440 : yim1 = 1
2954 [ + + ]: 468720 : do ii= 2, nn
2955 : 405280 : temp = yy
2956 : 405280 : yy = 2*xred*yy - yim1
2957 : 468720 : yim1 = temp
2958 : : end do
2959 : :
2960 : 63440 : end function cheb_poly1
2961 : : !!***
2962 : :
2963 : : !----------------------------------------------------------------------
2964 : :
2965 : : !!****f* m_chebfi2/chebfi_set_ndeg_from_residu
2966 : : !! NAME
2967 : : !! chebfi_set_ndeg_from_residu
2968 : : !!
2969 : : !! FUNCTION
2970 : : !! Compute ndeg_filter using the oracle and residuals.
2971 : : !!
2972 : : !! INPUTS
2973 : : !!
2974 : : !! OUTPUT
2975 : : !!
2976 : : !! SIDE EFFECTS
2977 : : !!
2978 : : !! SOURCE
2979 : :
2980 : 2312 : subroutine chebfi_set_ndeg_from_residu(chebfi,lambda_minus,lambda_plus,occ,DivResults,ndeg_filter_max,ndeg_filter)
2981 : :
2982 : : integer,intent(in) :: ndeg_filter_max
2983 : : integer,intent(out) :: ndeg_filter
2984 : : type(chebfi_t), intent(inout) :: chebfi
2985 : : type(xgBlock_t), intent(in) :: occ
2986 : : type(xgBlock_t), intent(in) :: DivResults
2987 : : real(dp), intent(in) :: lambda_minus, lambda_plus
2988 : :
2989 : : logical :: test1,test2,test3
2990 : : integer :: iband_tot,iband
2991 : : integer :: bandpp,ierr,ndeg_filter_tolwfr,ndeg_filter_decrease,nbdbuf,ndeg_filter_all,shift
2992 : 2312 : integer,allocatable :: ndeg_filter_bands(:)
2993 : : type(xgBlock_t) :: occBlock,occ_reshaped
2994 : : type(xg_t) :: residu
2995 : 2312 : real(dp),pointer :: residu_(:,:),occ_(:,:)
2996 : : real(dp) :: eig_iband,res_iband,occ_iband
2997 : 2312 : real(dp),pointer :: eig(:,:)
2998 : : ! character(len=500) :: msg
2999 : :
3000 : 2312 : bandpp = chebfi%bandpp
3001 : :
3002 : : !Compute residu here for oracle, use X_next as a work space
3003 : : ! X_next = S|Psi>
3004 : 2312 : call xgBlock_copy(chebfi%xBXColsRows,chebfi%X_next)
3005 : : ! X_next = - eig * S|Psi>
3006 : 2312 : call xgBlock_ymax(chebfi%X_next,DivResults,0,1)
3007 : : ! X_next = H|Psi> - eig * S|Psi>
3008 : 2312 : call xgBlock_add(chebfi%X_next,chebfi%xAXColsRows)
3009 : : ! resid = |X_next|^2
3010 : 2312 : call xg_init(residu,SPACE_R,bandpp,1)
3011 : 2312 : call xgBlock_colwiseNorm2(chebfi%X_next, residu%self,comm_loc=xmpi_comm_null)
3012 : :
3013 : 2312 : occ_reshaped = occ
3014 : 2312 : shift=xmpi_comm_rank(chebfi%comm_cols)*bandpp
3015 : 2312 : call xgBlock_reshape(occ_reshaped,1,chebfi%neigenpairs)
3016 : 2312 : call xgBlock_setBlock(occ_reshaped,occBlock,1,bandpp,fcol=1+shift)
3017 : 2312 : call xgBlock_reshape(occBlock,bandpp,1)
3018 [ + + ]: 2312 : if (chebfi%nbdbuf==-101) then
3019 : 952 : call xgBlock_apply_diag(residu%self,occBlock,1)
3020 : : end if
3021 : :
3022 [ - + + - : 6936 : ABI_MALLOC(ndeg_filter_bands,(bandpp))
+ - + - -
+ - + ]
3023 : :
3024 : : ! DivResults could be complex (with null imaginary part), so bandpp has to be in cols, not rows
3025 : 2312 : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=bandpp)
3026 : 2312 : call xgBlock_reverseMap(residu%self,residu_,rows=1,cols=bandpp)
3027 : 2312 : call xgBlock_reverseMap(occBlock,occ_,rows=1,cols=bandpp)
3028 : :
3029 [ + + ]: 2312 : if (chebfi%nbdbuf>0) then
3030 : : nbdbuf = chebfi%nbdbuf
3031 [ + - ]: 952 : else if (chebfi%nbdbuf==-101) then
3032 : 952 : nbdbuf = 0
3033 : : end if
3034 : :
3035 [ + + ]: 12104 : do iband=1, bandpp
3036 : 9792 : eig_iband = eig(1,iband)
3037 : 9792 : res_iband = residu_(1,iband)
3038 : 9792 : occ_iband = occ_(1,iband)
3039 : 9792 : iband_tot = iband + shift
3040 : 9792 : test1 = res_iband<chebfi%tolerance ! band already converged
3041 : 9792 : test2 = iband_tot>chebfi%neigenpairs-nbdbuf ! band in the buffer
3042 [ + + + + ]: 9792 : test3 = chebfi%nbdbuf==-101.and.occ_iband<chebfi%oracle_min_occ ! occupancy is too low
3043 [ + + - + ]: 12104 : if (test1.or.test2.or.test3) then
3044 : 5750 : ndeg_filter_bands(iband) = 0
3045 : : else
3046 : : !ndeg_filter necessary to converge to tolerance
3047 : 4042 : ndeg_filter_tolwfr = cheb_oracle1(eig_iband, lambda_minus, lambda_plus, chebfi%tolerance / res_iband, 1000)
3048 [ + + ]: 4042 : if (chebfi%oracle==1) then
3049 : 832 : ndeg_filter_bands(iband) = MIN(ndeg_filter_max, ndeg_filter_tolwfr, chebfi%ndeg_filter)
3050 [ + - ]: 3210 : else if (chebfi%oracle==2) then
3051 : : !ndeg_filter necessary to decrease residual by a constant factor
3052 : 3210 : ndeg_filter_decrease = cheb_oracle1(eig_iband, lambda_minus, lambda_plus, chebfi%oracle_factor, 15)
3053 : 3210 : ndeg_filter_bands(iband) = MIN(ndeg_filter_max, ndeg_filter_tolwfr, ndeg_filter_decrease)
3054 : : else
3055 : 0 : ABI_ERROR('Wrong value for chebfi%oracle')
3056 : : end if
3057 : : end if
3058 : : end do
3059 [ + + ]: 12104 : ndeg_filter = MAXVAL(ndeg_filter_bands)
3060 : 2312 : call xmpi_max(ndeg_filter,ndeg_filter_all,chebfi%comm_cols,ierr)
3061 : 2312 : ndeg_filter=ndeg_filter_all
3062 : :
3063 : 2312 : call xg_free(residu)
3064 [ + - ]: 2312 : ABI_SFREE(ndeg_filter_bands)
3065 : :
3066 : 4624 : end subroutine chebfi_set_ndeg_from_residu
3067 : : !!***
3068 : :
3069 : : !----------------------------------------------------------------------
3070 : :
3071 : : !!****f* m_chebfi2/chebfi_getSubspaceResidual
3072 : : !! NAME
3073 : : !! chebfi_getSubspaceResidual
3074 : : !!
3075 : : !! FUNCTION
3076 : : !! || (I-QQ'B)Aq_j ||^2 for j=1,..,bandpp
3077 : : !! stored to resid
3078 : :
3079 : 0 : subroutine chebfi_getSubspaceResidual(chebfi, resid)
3080 : :
3081 : : implicit none
3082 : :
3083 : : type(chebfi_t), intent(inout) :: chebfi
3084 : : type(xgBlock_t), intent(inout) :: resid
3085 : :
3086 : : type(xg_t) :: M, R
3087 : : integer :: nrows, ncols
3088 : : integer :: space_buf
3089 : :
3090 : : ! *********************************************************************
3091 : :
3092 : 0 : nrows = rows(chebfi%X) ! todo this should be chebfi%spacedim
3093 : 0 : ncols = chebfi%neigenpairs ! todo active part only
3094 : :
3095 : 0 : if (chebfi%space==SPACE_C) then
3096 : 0 : space_buf = SPACE_C
3097 [ # # ]: 0 : else if (chebfi%space==SPACE_CR) then
3098 : 0 : space_buf = SPACE_R
3099 : : else
3100 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
3101 : : end if
3102 : :
3103 : 0 : call xg_init(M, space_buf, ncols, ncols, comm=chebfi%spacecom, gpu_option=chebfi%gpu_option)
3104 : 0 : call xg_init(R, space_buf, nrows, ncols, comm=chebfi%spacecom, gpu_option=chebfi%gpu_option)
3105 : :
3106 : : ! Compute M = X^T*AX
3107 : : ! sum all process contribution
3108 : 0 : call xgBlock_gemm('t','n',1.0d0,chebfi%X,chebfi%AX%self,0.d0,M%self,comm=chebfi%spacecom)
3109 : :
3110 : : ! Compute R = AX - BX*M
3111 : 0 : call xgBlock_copy(chebfi%AX%self, R%self)
3112 : 0 : call xgBlock_gemm('n','n',-1.0d0,chebfi%BX%self,M%self,1.0d0,R%self)
3113 : :
3114 : : ! sum along rows communicator
3115 : 0 : call xgBlock_colwiseNorm2(R%self, resid)
3116 : :
3117 : 0 : call xg_free(M)
3118 : 0 : call xg_free(R)
3119 : :
3120 [ # # ]: 0 : end subroutine chebfi_getSubspaceResidual
3121 : : !!***
3122 : :
3123 : : !----------------------------------------------------------------------
3124 : :
3125 : : !!****f* m_chebfi2/chebfi_swapConvergedVectors
3126 : : !! NAME
3127 : : !! chebfi_swapConvergedVectors
3128 : : !!
3129 : : !! FUNCTION
3130 : : !! In-place contiguous-in-memory locking of column vectors in row distribution.
3131 : : !! Algorithm is based on swapping: at the end first k are locked last are active.
3132 : : !! | locked columns | active columns |
3133 : : !! | 1 ... k | k+1 ... m |
3134 : : !! TODO logic could also be applied to divide Xext to slices.
3135 : : !!
3136 : :
3137 : 0 : subroutine chebfi_swapConvergedVectors(chebfi, resid, tol, n_locked)
3138 : :
3139 : : implicit none
3140 : :
3141 : : type(chebfi_t), intent(inout) :: chebfi
3142 : : type(xgBlock_t), intent(in) :: resid
3143 : : real(dp), intent(in) :: tol
3144 : : integer, intent(out) :: n_locked
3145 : :
3146 : : integer :: i, nrows
3147 : : integer :: j, left
3148 : 0 : real(dp), pointer :: resid_vals(:,:)
3149 : 0 : logical, allocatable :: mask(:)
3150 : 0 : logical, allocatable :: is_locked(:)
3151 : 0 : integer, allocatable :: idx(:)
3152 : : integer :: space_buf
3153 : : type(xg_t) :: T_swap ! temporary column buffer
3154 : :
3155 : : ! *********************************************************************
3156 : :
3157 : 0 : nrows = rows(resid)
3158 : :
3159 : 0 : if (chebfi%space==SPACE_C) then
3160 : 0 : space_buf = SPACE_C
3161 [ # # ]: 0 : else if (chebfi%space==SPACE_CR) then
3162 : 0 : space_buf = SPACE_R
3163 : : else
3164 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
3165 : : end if
3166 : :
3167 [ # # # # : 0 : ABI_MALLOC(mask, (nrows))
# # # # #
# # # ]
3168 [ # # # # : 0 : ABI_MALLOC(is_locked, (nrows))
# # # # ]
3169 : :
3170 : 0 : call xg_init(T_swap, space_buf, nrows, 1, comm=chebfi%spacecom, gpu_option=chebfi%gpu_option)
3171 : :
3172 : 0 : call xgBlock_reverseMap(resid, resid_vals, rows=nrows, cols=1)
3173 : :
3174 [ # # # # : 0 : mask = resid_vals(:,1) < tol
# # # # ]
3175 [ # # # # : 0 : idx = pack([(i, i=1,nrows)], mask)
# # # # ]
3176 : 0 : n_locked = size(idx)
3177 : :
3178 : 0 : write(std_out,*) 'n_locked=', n_locked
3179 : 0 : write(std_out,*) 'locked indices=', idx(:)
3180 : 0 : write(std_out,*) 'residuals=', resid_vals(:,1)
3181 : 0 : flush(std_out)
3182 : :
3183 [ # # ]: 0 : is_locked = .false.
3184 [ # # ]: 0 : is_locked(idx) = .true.
3185 : :
3186 : : ! partition via swapping (in-place)
3187 : 0 : left = 1
3188 [ # # ]: 0 : do j = 1, nrows
3189 [ # # ]: 0 : if (is_locked(j)) then
3190 [ # # ]: 0 : if (j /= left) then
3191 : 0 : call xgBlock_colwiseSwap(chebfi%X, j, left, T_swap%self)
3192 : 0 : call xgBlock_colwiseSwap(chebfi%AX%self, j, left, T_swap%self)
3193 : 0 : call xgBlock_colwiseSwap(chebfi%BX%self, j, left, T_swap%self)
3194 : :
3195 : : ! keep mask consistent after swap
3196 : 0 : is_locked(j) = is_locked(left)
3197 : 0 : is_locked(left) = .true.
3198 : : end if
3199 : 0 : left = left + 1
3200 : : end if
3201 : : end do
3202 : :
3203 [ # # ]: 0 : ABI_FREE(mask)
3204 [ # # ]: 0 : ABI_FREE(is_locked)
3205 : 0 : call xg_free(T_swap)
3206 : :
3207 [ # # # # ]: 0 : end subroutine chebfi_swapConvergedVectors
3208 : : !!***
3209 : :
3210 : : !----------------------------------------------------------------------
3211 : :
3212 : : !!****f* m_chebfi2/chebfi_deflateWrtLocked
3213 : : !! NAME
3214 : : !! chebfi_deflateWrtLocked
3215 : : !!
3216 : :
3217 : 0 : subroutine chebfi_deflateWrtLocked(chebfi, n_locked, X_locked, BX_locked)
3218 : :
3219 : : type(chebfi_t), intent(inout) :: chebfi
3220 : : integer, intent(in) :: n_locked
3221 : : type(xgBlock_t), intent(inout) :: X_locked
3222 : : type(xgBlock_t), intent(inout) :: BX_locked
3223 : :
3224 : : integer :: nrows
3225 : : integer :: space_buf
3226 : : integer :: n_active
3227 : : type(xg_t) :: buffer
3228 : : type(xgBlock_t) :: X_active, AX_active, BX_active
3229 : :
3230 : : ! *********************************************************************
3231 : :
3232 : 0 : if (chebfi%space==SPACE_C) then
3233 : 0 : space_buf = SPACE_C
3234 [ # # ]: 0 : else if (chebfi%space==SPACE_CR) then
3235 : 0 : space_buf = SPACE_R
3236 : : else
3237 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
3238 : : end if
3239 : :
3240 : 0 : nrows = rows(chebfi%X)
3241 : 0 : n_active = chebfi%neigenpairs - n_locked
3242 : 0 : call xg_init(buffer,space_buf,n_locked,n_active,comm=chebfi%spacecom,gpu_option=chebfi%gpu_option)
3243 : :
3244 : 0 : write(std_out,*) 'start gemm with sum over comm'; flush(std_out)
3245 : 0 : write(std_out,*) 'getid', xgBlock_getid(buffer%self); flush(std_out)
3246 : 0 : write(std_out,*) 'getid', xgBlock_getid(buffer%self,comm=xmpi_comm_null); flush(std_out)
3247 : 0 : write(std_out,*) 'getid', xgBlock_getid(buffer%self,comm=chebfi%spacecom); flush(std_out)
3248 : :
3249 : : ! buffer = BX0^T*X
3250 : : ! sum all process contribution
3251 : 0 : call xgBlock_setBlock(chebfi%X, X_active, nrows, n_active, fcol=n_locked+1)
3252 : 0 : call xgBlock_gemm('t','n',1.0d0,BX_locked,X_active,0.d0,buffer%self, comm=chebfi%spacecom)
3253 : :
3254 : 0 : write(std_out,*) 'exit gemm with sum over comm'; flush(std_out)
3255 : :
3256 : 0 : write(std_out,*) '@deflate W', rows(buffer%self), cols(buffer%self)
3257 : 0 : write(std_out,*) '@deflate A', rows(BX_locked), cols(BX_locked)
3258 : 0 : write(std_out,*) '@deflate B', rows(X_active), cols(X_active)
3259 : 0 : flush(std_out)
3260 : :
3261 : : ! X = X - X0*buffer
3262 : 0 : call xgBlock_gemm('n','n',-1.0d0,X_locked,buffer%self,1.0d0,X_active)
3263 : : ! AX = AX - X0*buffer
3264 : 0 : call xgBlock_setBlock(chebfi%AX%self, AX_active, nrows, n_active, fcol=n_locked+1)
3265 : 0 : call xgBlock_gemm('n','n',-1.0d0,X_locked,buffer%self,1.0d0,AX_active)
3266 : : ! BX = BX - X0*buffer
3267 : 0 : call xgBlock_setBlock(chebfi%BX%self, BX_active, nrows, n_active, fcol=n_locked+1)
3268 : 0 : call xgBlock_gemm('n','n',-1.0d0,X_locked,buffer%self,1.0d0,BX_active)
3269 : :
3270 : 0 : call xg_free(buffer)
3271 : :
3272 [ # # ]: 0 : end subroutine chebfi_deflateWrtLocked
3273 : : !!***
3274 : :
3275 : : !----------------------------------------------------------------------
3276 : :
3277 : : !!****f* m_chebfi2/chebfi_deflateActiveWrtLocked
3278 : : !! NAME
3279 : : !! chebfi_deflateActiveWrtLocked
3280 : : !!
3281 : : !! FUNCTION
3282 : : !! Deflate active vectors with respect to locked vectors in linalg representation
3283 : :
3284 : 0 : subroutine chebfi_deflateActiveWrtLocked(chebfi, n_locked, X_part, AX_part, BX_part)
3285 : :
3286 : : type(chebfi_t), intent(inout) :: chebfi
3287 : : integer, intent(in) :: n_locked
3288 : : type(bandPartitionData_t), intent(inout) :: X_part
3289 : : type(bandPartitionData_t), intent(inout) :: AX_part
3290 : : type(bandPartitionData_t), intent(inout) :: BX_part
3291 : :
3292 : : integer :: space_buf
3293 : : integer :: n_active
3294 : : type(xg_t) :: buffer
3295 : :
3296 : : ! *********************************************************************
3297 : :
3298 : 0 : if (chebfi%space==SPACE_C) then
3299 : 0 : space_buf = SPACE_C
3300 [ # # ]: 0 : else if (chebfi%space==SPACE_CR) then
3301 : 0 : space_buf = SPACE_R
3302 : : else
3303 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
3304 : : end if
3305 : :
3306 : 0 : n_active = chebfi%neigenpairs - n_locked
3307 : 0 : call xg_init(buffer,space_buf,n_locked,n_active,comm=chebfi%spacecom,gpu_option=chebfi%gpu_option)
3308 : :
3309 : 0 : write(std_out,*) '@deflate W', rows(buffer%self), cols(buffer%self)
3310 : 0 : write(std_out,*) '@deflate A', rows(BX_part%linalg_locked), cols(BX_part%linalg_locked)
3311 : 0 : write(std_out,*) '@deflate B', rows(X_part%linalg_active), cols(X_part%linalg_active)
3312 : 0 : flush(std_out)
3313 : :
3314 : : ! buffer = BX0^T*X
3315 : : call xgBlock_gemm('t','n',1.0d0,BX_part%linalg_locked,X_part%linalg_active,0.d0,&
3316 : 0 : buffer%self,comm=chebfi%spacecom)
3317 : :
3318 : : ! sum all process contribution
3319 : : ! X = X - X0*buffer
3320 : 0 : call xgBlock_gemm('n','n',-1.0d0,X_part%linalg_locked,buffer%self,1.0d0,X_part%linalg_active)
3321 : : ! AX = AX - X0*buffer
3322 : 0 : call xgBlock_gemm('n','n',-1.0d0,X_part%linalg_locked,buffer%self,1.0d0,AX_part%linalg_active)
3323 : : ! BX = BX - X0*buffer
3324 : 0 : call xgBlock_gemm('n','n',-1.0d0,X_part%linalg_locked,buffer%self,1.0d0,BX_part%linalg_active)
3325 : :
3326 : 0 : call xg_free(buffer)
3327 : :
3328 [ # # ]: 0 : end subroutine chebfi_deflateActiveWrtLocked
3329 : : !!***
3330 : :
3331 : : !----------------------------------------------------------------------
3332 : :
3333 : : !!****f* m_chebfi2/bandPartitionInfo_initActiveSubcomm
3334 : : !! NAME
3335 : : !! bandPartitionInfo_initActiveSubcomm
3336 : : !!
3337 : : !! FUNCTION
3338 : : !! Every process can treat between 1 and bandpp bands
3339 : : !! with p processes we cover p x bandpp bands at maximum.
3340 : : !! Solve problem find minimal p to cover n_active bands.
3341 : : !! Constraints: use at least 1 process and at max nproc.
3342 : : !!
3343 : : !! INPUT
3344 : : !! comm global communicator
3345 : : !!
3346 : : !! OUTPUT
3347 : : !! mpiInfo%n_active_bandpp stores the result per mpi rank
3348 : : !! mpiInfo%rank_active true if process treats active bands
3349 : :
3350 : 0 : subroutine bandPartitionInfo_initActiveSubcomm(mpiInfo, comm, bandpp_max)
3351 : :
3352 : : implicit none
3353 : :
3354 : : type(bandPartitionInfo_t), intent(inout) :: mpiInfo
3355 : : integer, intent(in) :: comm
3356 : : integer, intent(in) :: bandpp_max
3357 : :
3358 : 0 : integer, allocatable :: ranks(:)
3359 : : integer :: nproc
3360 : : integer :: nproc_active, min_p, ii
3361 : :
3362 : : ! *********************************************************************
3363 : :
3364 : 0 : nproc = xmpi_comm_size(comm)
3365 [ # # ]: 0 : min_p = ceiling(1.d0 * mpiInfo%n_active / bandpp_max)
3366 : 0 : nproc_active = min(nproc, max(1, min_p))
3367 : :
3368 [ # # # # : 0 : ABI_MALLOC(ranks,(0:nproc_active-1))
# # # # #
# # # ]
3369 [ # # # # : 0 : ranks(0:nproc_active-1)=(/((ii),ii=0,nproc_active-1)/)
# # ]
3370 : 0 : mpiInfo%comm_active=xmpi_subcomm(comm,nproc_active,ranks)
3371 [ # # ]: 0 : ABI_FREE(ranks)
3372 : :
3373 : : ! todo @debug all processes of slice must create the same active subcomm
3374 : 0 : write(std_out,*) '@initActiveSubcomm comm_active=', mpiInfo%comm_active
3375 : 0 : flush(std_out)
3376 : :
3377 : :
3378 : 0 : end subroutine bandPartitionInfo_initActiveSubcomm
3379 : : !!***
3380 : :
3381 : : !----------------------------------------------------------------------
3382 : :
3383 : : !!****f* m_chebfi2/bandPartitionData_allocActiveColsrows
3384 : : !! NAME
3385 : : !! bandPartitionData_allocActiveColsrows
3386 : :
3387 : 0 : subroutine bandPartitionData_allocActiveColsrows(bpart, chebfi, mpiInfo, nspinor)
3388 : :
3389 : : implicit none
3390 : :
3391 : : type(bandPartitionData_t), intent(inout) :: bpart
3392 : : type(chebfi_t), intent(in) :: chebfi
3393 : : type(bandPartitionInfo_t), intent(inout) :: mpiInfo
3394 : : integer, intent(in) :: nspinor
3395 : :
3396 : : integer :: n_active
3397 : : integer :: comm_active
3398 : : real(dp) :: tsec(2)
3399 : :
3400 : : ! *********************************************************************
3401 : :
3402 : 0 : n_active = mpiInfo%n_active
3403 : 0 : comm_active = mpiInfo%comm_active
3404 : :
3405 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_TRANSPOSE)
3406 : 0 : call timab(tim_transpose,1,tsec)
3407 : :
3408 : : call xgTransposer_constructor(bpart%transposer_active,bpart%linalg_active,bpart%colsrows_active,&
3409 : : nspinor,STATE_LINALG,TRANS_ALL2ALL,chebfi%comm_rows,comm_active,0,0,chebfi%me_g0_fft,&
3410 : 0 : gpu_option=chebfi%gpu_option,gpu_thread_limit=chebfi%gpu_thread_limit)
3411 : :
3412 : 0 : bpart%transposer_active%gpu_kokkos_nthrd = chebfi%gpu_kokkos_nthrd
3413 : :
3414 : : ABI_NVTX_END_RANGE()
3415 : 0 : call timab(tim_transpose,2,tsec)
3416 : :
3417 : 0 : end subroutine bandPartitionData_allocActiveColsrows
3418 : : !!***
3419 : :
3420 : : !----------------------------------------------------------------------
3421 : :
3422 : : !!****f* m_chebfi2/bandPartitionData_setLinalg
3423 : : !! NAME
3424 : : !! bandPartitionData_setLinalg
3425 : : !!
3426 : : !! FUNCTION
3427 : : !! Set active and locked vectors in linalg representation
3428 : :
3429 : 0 : subroutine bandPartitionData_setLinalg(X_in, bpart, n_locked)
3430 : :
3431 : : implicit none
3432 : :
3433 : : type(xgBlock_t), intent(in) :: X_in
3434 : : type(bandPartitionData_t), intent(inout) :: bpart
3435 : : integer, intent(in) :: n_locked
3436 : :
3437 : : integer :: nrows
3438 : : integer :: n_active
3439 : :
3440 : : ! *********************************************************************
3441 : :
3442 : 0 : nrows = rows(X_in)
3443 : 0 : n_active = cols(X_in) - n_locked
3444 : :
3445 : 0 : call xgBlock_setBlock(X_in, bpart%linalg_locked, nrows, n_locked)
3446 : 0 : call xgBlock_setBlock(X_in, bpart%linalg_active, nrows, n_active, fcol=n_locked+1)
3447 : :
3448 : 0 : end subroutine bandPartitionData_setLinalg
3449 : : !!***
3450 : :
3451 [ # # # # : 0 : end module m_chebfi2
# # # # #
# # # #
# ]
3452 : : !!***
|
