Branch data Line data Source code
1 : : !!****f* ABINIT/m_slice
2 : : !! NAME
3 : : !! m_slice
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the types and routines used to apply the Spectrum Slicing
7 : : !! method. It mainly defines 'slice' datatypes and associated methods.
8 : : !!
9 : : !! Main features:
10 : : !! - uses 'xgTools' implementation for matrix data structure.
11 : : !! - adopts 'chebfi' functionalities for most matrix calculations.
12 : : !! - uses polynomial filtering, Chebyshev lowpass, next Chebyshev-Jackson.
13 : : !! - uses scheduler and resource allocator for slice tasks.
14 : : !! - uses logic for asynchronous slice memory avoiding race condition in read/write.
15 : : !! - applies Rayleigh-Ritz for individual slices in parallel or sequentially.
16 : : !!
17 : : !! NOTES
18 : : !! Dependence with other modules and hierarchy in this directory:
19 : : !! m_slice uses:
20 : : !! |- m_polynomial_filter | various math compute routines
21 : : !! |- m_trace_estimation | various math compute routines
22 : : !! |- m_chebfi2 ! the Chebyshev recursion
23 : : !! |- m_slice_task | Low-level
24 : : !!
25 : : !! COPYRIGHT
26 : : !! Copyright (C) 2018-2026 ABINIT group (IML)
27 : : !! This file is distributed under the terms of the
28 : : !! gnu general public license, see ~abinit/COPYING
29 : : !! or http://www.gnu.org/copyleft/gpl.txt .
30 : : !! for the initials of contributors, see ~abinit/doc/developers/contributors.txt .
31 : : !!
32 : : !! SOURCE
33 : :
34 : : #if defined HAVE_CONFIG_H
35 : : #include "config.h"
36 : : #endif
37 : :
38 : : #include "abi_common.h"
39 : :
40 : : ! nvtx related macro definition
41 : : #include "nvtx_macros.h"
42 : :
43 : : module m_slice
44 : :
45 : : use defs_basis
46 : : use defs_abitypes
47 : : use m_abicore
48 : : use m_errors
49 : : use m_time, only : timab
50 : : use m_sort, only: sort_dp
51 : :
52 : : use m_cgtools
53 : : use m_xg
54 : : use m_xgTransposer
55 : : use m_xg_ortho_RR
56 : : use m_chebfi2
57 : :
58 : : use m_trace_estimation
59 : : use m_slice_task
60 : : use m_polynomial_filter
61 : :
62 : : use m_xmpi
63 : : use m_xomp
64 : : #ifdef HAVE_OPENMP
65 : : use omp_lib
66 : : #endif
67 : :
68 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
69 : : use m_gpu_toolbox, only : CPU_DEVICE_ID, gpu_device_synchronize
70 : : #endif
71 : :
72 : : #if defined(HAVE_GPU_MARKERS)
73 : : use m_nvtx_data
74 : : #endif
75 : :
76 : : implicit none
77 : :
78 : : private
79 : :
80 : : ! Spectral interval decomposition strategy (spectral_cut input variable)
81 : : ! todo delete this obsolete
82 : : !---------------------------------------------------
83 : : integer, parameter :: DIVIDE_INTERVAL_WIDTH = 0
84 : : integer, parameter :: DIVIDE_NUMBER_OF_VECTORS = 1
85 : : integer, parameter :: DIVIDE_SPECTRAL_GAPS = 2
86 : :
87 : : ! Load balance criteria for fair resource allocation (paral_slice input variable)
88 : : ! todo should be either 0 or 1
89 : : !---------------------------------------------------
90 : : integer, parameter :: FAIR_BANDPP = 1 ! bandpp (unweigted)
91 : : integer, parameter :: FAIR_BANDPP_WDEG = 2 ! bandpp weighted by degree
92 : : integer, parameter :: EVEN_SLICES = 3 ! all slices have the same num of procs
93 : : integer, parameter :: SEQUENTIAL_SLICE = 4 ! sequential slice treatment
94 : :
95 : : integer, parameter :: ENABLE_PARAL = 0
96 : : integer, parameter :: DISABLE_PARAL = 1
97 : :
98 : : ! Timers
99 : : !---------------------------------------------------
100 : : integer, parameter :: tim_swap = 1761
101 : : integer, parameter :: tim_RR_q = 1759
102 : : integer, parameter :: tim_barrier = 1764
103 : : integer, parameter :: tim_copy = 1765
104 : : integer, parameter :: tim_Bortho_X = 1641
105 : : integer, parameter :: tim_getAX_BX = 1754
106 : : integer, parameter :: tim_invovl = 1755
107 : : integer, parameter :: tim_slice_sched = 2161
108 : : integer, parameter :: tim_slice_1 = 2162
109 : : integer, parameter :: tim_slice_2 = 2163
110 : : integer, parameter :: tim_slice_3 = 2164
111 : : integer, parameter :: tim_slice_X = 2165
112 : : integer, parameter :: tim_lanczos = 2167
113 : : integer, parameter :: tim_trace = 2168
114 : :
115 : : ! Public 'slice' datatype
116 : : !-------------------------------------------------
117 : : type, public :: slice_t
118 : :
119 : : ! MPI-related common to all-slice phase
120 : : integer :: nproc ! number of available processes
121 : : integer :: comm_rows ! xmpi_comm_self ...
122 : : integer :: comm_cols ! same as spacecom
123 : : integer :: spacecom ! same as comm_cols
124 : : integer :: bandpp ! nb of bands per process in colsrows representation
125 : : integer :: me_g0 ! 1 if this processors treats G=0, 0 otherwise
126 : : integer :: me_g0_fft ! 1 if this processors treats G=0 in FFT, 0 otherwise
127 : :
128 : : ! Eigenpair parameters
129 : : integer :: neigenpairs ! total number of bands (=number of eigenpairs)
130 : : integer :: total_spacedim ! total number of plane-waves
131 : : integer :: spacedim ! nb of plane-waves per process in linalg representation
132 : : integer :: nslice ! number of spectral slices
133 : : integer :: space ! real or complex eigenvectors
134 : : integer :: nbdbuf ! fixme for the moment not clear if useful
135 : :
136 : : ! GPU-related
137 : : integer :: gpu_kokkos_nthrd
138 : : integer :: gpu_thread_limit
139 : :
140 : : ! Various options
141 : : integer :: gpu_option ! enable GPU
142 : : integer :: paral_kgb ! enable parallel (k-points, G basis, bands)
143 : : integer :: paral_slice ! how to allocate resources for parallel slices
144 : : integer :: spectral_cut ! how to decompose spectrum into slices
145 : :
146 : : ! Various model parameters
147 : : logical :: paw ! use PAW or not
148 : : integer :: ndeg_filter ! lowpass degree of polynomial filter
149 : : real(dp) :: tolerance ! tolerance on the residu to stop the minimization
150 : : real(dp) :: ramp ! bandpass filter tolerance
151 : : real(dp) :: ecut ! Ecut Fermi level
152 : : real(dp) :: mineig ! slightly larger than mineig_global
153 : : real(dp) :: mineig_global ! guaranteed lower spectral bound
154 : : real(dp) :: maxeig_global ! guaranteed upper spectral bound
155 : :
156 : : ! Arrays related to polynomial filtering
157 : : integer, allocatable :: neigen_per_slice(:) ! total number of eigenpairs per slice
158 : : integer, allocatable :: poly_degrees(:) ! polynomial filter degrees per slice
159 : : real(dp), allocatable :: low_bounds(:) ! lower bounds in spectral partition (disjoint)
160 : : real(dp), allocatable :: upp_bounds(:) ! upper bounds in spectral partition (disjoint)
161 : :
162 : : end type slice_t
163 : :
164 : : ! Public methods associated to 'slice' datatype
165 : : !-------------------------------------------------
166 : : public :: slice_init ! initialize slice datatype object
167 : : public :: slice_free ! free slice datatype object
168 : : public :: slice_run ! run Spectrum Slicing for active slice task
169 : :
170 : : CONTAINS
171 : : !=====================================================================
172 : : !!***
173 : :
174 : : !!****f* m_slice/slice_init
175 : : !! NAME
176 : : !! slice_init
177 : : !!
178 : : !! FUNCTION
179 : : !! Initialize a 'slice' datastructure.
180 : : !!
181 : : !! INPUTS
182 : : !! nslice= number of spectral slices
183 : : !! neigenpairs= number of requested eigenvectors/eigenvalues
184 : : !! spacedim= space dimension for one vector
185 : : !! tolerance= tolerance criterion on the residu to stop the minimization
186 : : !! ecut= plane-wave cut-off energy
187 : : !! paral_kgb= flag controlling (k,g,bands) parallelization
188 : : !! bandpp= number of 'bands' handled by a processor
189 : : !! ndeg_filter= polynomial degree of the polynomial filter (.i.e. number of H applications)
190 : : !! is_lowpass= flag. True: use lowpass Chebyshev, false: use bandpass Heaviside expanded on Chebyshev
191 : : !! space= defines in which space we are (columns, rows, etc.)
192 : : !! eigenProblem= type of eigenpb: 1 (A*x = (lambda)*B*x), 2 (A*B*x = (lambda)*x), 3 (B*A*x = (lambda)*x)
193 : : !! spacecom= MPI communicator
194 : : !! me_g0= 1 if this processors treats G=0, 0 otherwise
195 : : !! me_g0_ftt= 1 if this processors treats G=0 in FFT, 0 otherwise
196 : : !! paw= flag. TRUE if current calculation ses the PAW approach
197 : : !! comm_rows= "rows" communicator
198 : : !! comm_cols= "cols" communicator
199 : : !! gpu_option= flag. Enable GPU if true
200 : : !! gpu_kokkos_nthrd= number of OpenMP offloaded threads used
201 : : !! gpu_thread_limit= maximum number of OpenMP offloaded threads
202 : : !!
203 : : !! SOURCE
204 : :
205 : 18 : subroutine slice_init(slice,nslice,neigenpairs,spacedim,tolerance,paral_kgb,&
206 : : paral_slice,ndeg_filter,nbdbuf,ramp,ecut,bandpp,space,spacecom,me_g0,me_g0_fft,&
207 : : paw,comm_rows,comm_cols,spectral_cut,gpu_option,gpu_kokkos_nthrd,gpu_thread_limit)
208 : :
209 : : implicit none
210 : :
211 : : ! Arguments ------------------------------------
212 : : integer , intent(in ) :: bandpp
213 : : integer , intent(in ) :: nslice
214 : : integer , intent(in ) :: me_g0
215 : : integer , intent(in ) :: me_g0_fft
216 : : integer , intent(in ) :: neigenpairs
217 : : integer , intent(in ) :: comm_rows
218 : : integer , intent(in ) :: comm_cols
219 : : integer , intent(in ) :: paral_kgb
220 : : integer , intent(in ) :: paral_slice
221 : : integer , intent(in ) :: space
222 : : integer , intent(in ) :: spacecom
223 : : integer , intent(in ) :: spacedim
224 : : integer , intent(in ) :: ndeg_filter
225 : : integer , intent(in ) :: nbdbuf
226 : : integer , intent(in ) :: spectral_cut
227 : : integer , intent(in ) :: gpu_option
228 : : logical , intent(in ) :: paw
229 : : real(dp) , intent(in ) :: ramp
230 : : real(dp) , intent(in ) :: tolerance
231 : : real(dp) , intent(in ) :: ecut
232 : : type(slice_t), intent(inout) :: slice
233 : : integer , intent(in ), optional :: gpu_kokkos_nthrd
234 : : integer , intent(in ), optional :: gpu_thread_limit
235 : :
236 : : ! Local variables --------------------------------
237 : : integer :: total_spacedim, ierr
238 : : !logical :: on_host, on_device
239 : :
240 : : ! *********************************************************************
241 : :
242 : 18 : slice%bandpp = bandpp
243 : 18 : slice%nslice = nslice
244 : 18 : slice%me_g0 = me_g0
245 : 18 : slice%me_g0_fft = me_g0_fft
246 : 18 : slice%neigenpairs = neigenpairs
247 : 18 : slice%comm_rows = comm_rows
248 : 18 : slice%comm_cols = comm_cols
249 : 18 : slice%paral_kgb = paral_kgb
250 : 18 : slice%paral_slice = paral_slice
251 : 18 : slice%space = space
252 : 18 : slice%spacecom = spacecom
253 : 18 : slice%spacedim = spacedim
254 : 18 : slice%ndeg_filter = ndeg_filter
255 : 18 : slice%nbdbuf = nbdbuf
256 : 18 : slice%spectral_cut = spectral_cut
257 : 18 : slice%gpu_option = gpu_option
258 : 18 : slice%paw = paw
259 : 18 : slice%ramp = ramp
260 : 18 : slice%tolerance = tolerance
261 : 18 : slice%ecut = ecut
262 : :
263 : 18 : slice%gpu_kokkos_nthrd = 1
264 [ + - ]: 18 : if (present(gpu_kokkos_nthrd)) then
265 : 18 : slice%gpu_kokkos_nthrd = gpu_kokkos_nthrd
266 : : end if
267 : 18 : slice%gpu_thread_limit = 0
268 [ + - ]: 18 : if (present(gpu_thread_limit)) then
269 : 18 : slice%gpu_thread_limit = gpu_thread_limit
270 : : end if
271 : :
272 : : ! Total number of rows (used in colsrows representation)
273 [ + + ]: 18 : if (paral_kgb==0) then
274 : 2 : slice%total_spacedim = spacedim
275 : : else
276 : 16 : total_spacedim = spacedim
277 : 16 : call xmpi_sum(total_spacedim, spacecom, ierr)
278 : 16 : slice%total_spacedim = total_spacedim
279 : : end if
280 : :
281 : : ! Total number of processes
282 : 18 : slice%nproc = xmpi_comm_size(spacecom)
283 : :
284 : : ! Arrays
285 : 18 : call slice_allocateAll(slice)
286 : :
287 : 18 : end subroutine slice_init
288 : : !!***
289 : :
290 : : !----------------------------------------------------------------------
291 : :
292 : : !!****f* m_slice/slice_allocateAll
293 : : !! name
294 : : !! slice_allocateAll
295 : :
296 : 18 : subroutine slice_allocateAll(slice)
297 : :
298 : : implicit none
299 : : type(slice_t), intent(inout) :: slice
300 : :
301 : : ! *********************************************************************
302 : :
303 : 18 : call slice_free(slice)
304 : :
305 [ + - - + : 54 : ABI_MALLOC_IFNOT(slice%neigen_per_slice, (slice%nslice))
+ - + - +
- - + - +
- + ]
306 [ + - - + : 54 : ABI_MALLOC_IFNOT(slice%poly_degrees, (slice%nslice))
+ - + - +
- - + - +
- + ]
307 [ + - - + : 54 : ABI_MALLOC_IFNOT(slice%low_bounds, (slice%nslice))
+ - + - +
- - + - +
- + ]
308 [ + - - + : 54 : ABI_MALLOC_IFNOT(slice%upp_bounds, (slice%nslice))
+ - + - +
- - + - +
- + ]
309 : :
310 : 18 : end subroutine slice_allocateAll
311 : : !!***
312 : :
313 : : !----------------------------------------------------------------------
314 : :
315 : : !!****f* m_slice/slice_free
316 : : !! name
317 : : !! slice_free
318 : :
319 : 36 : subroutine slice_free(slice)
320 : :
321 : : implicit none
322 : : type(slice_t), intent(inout) :: slice
323 : :
324 : : ! *********************************************************************
325 : :
326 [ + + ]: 36 : ABI_SFREE(slice%neigen_per_slice)
327 [ + + ]: 36 : ABI_SFREE(slice%poly_degrees)
328 [ + + ]: 36 : ABI_SFREE(slice%low_bounds)
329 [ + + ]: 36 : ABI_SFREE(slice%upp_bounds)
330 : :
331 : 36 : end subroutine slice_free
332 : : !!***
333 : :
334 : : !----------------------------------------------------------------------
335 : :
336 : : !!****f* m_slice/slice_run
337 : : !! NAME
338 : : !! slice_run
339 : : !!
340 : : !! FUNCTION
341 : : !! Run Spectrum slicing on a given set of active vectors.
342 : : !! Computation uses marked resources for current task **only**.
343 : : !!
344 : : !! INPUTS
345 : : !! X0= size (spacedim, neigenpairs)
346 : : !! eigen,residu= size (neigenpairs, 1) (column vectors)
347 : : !! getAX_BX= pointer to the function giving A|X> and B|X>
348 : : !! A is typically the Hamiltonian H, and B the overlap operator S
349 : : !! getBm1X= pointer to the function giving B^-1|X>
350 : : !! B is typically the overlap operator S
351 : : !! nspinor= number of spinorial components of the wavefunctions
352 : : !!
353 : : !! SIDE EFFECTS
354 : : !! slice <type(slice_t)>= memory workspace used for Spectrum slicing
355 : : !!
356 : : !! SOURCE
357 : :
358 : 18 : subroutine slice_run(slice, X0, getAX_BX, getBm1X, eigen, residu, nspinor)
359 : :
360 : : implicit none
361 : :
362 : : ! Arguments
363 : : type(slice_t), target, intent(inout) :: slice
364 : : type(xgBlock_t), intent(inout) :: X0 ! size (spacedim, neigenpairs)
365 : : type(xgBlock_t), intent(inout) :: eigen ! size (neigenpairs,1)
366 : : type(xgBlock_t), intent(inout) :: residu ! size (neigenpairs,1)
367 : : integer, intent(in) :: nspinor
368 : : interface
369 : : subroutine getAX_BX(X,AX,BX)
370 : : use m_xg, only : xgBlock_t
371 : : type(xgBlock_t), intent(inout) :: X
372 : : type(xgBlock_t), intent(inout) :: AX
373 : : type(xgBlock_t), intent(inout) :: BX
374 : : end subroutine getAX_BX
375 : : end interface
376 : : interface
377 : : subroutine getBm1X(X,Bm1X)
378 : : use m_xg, only : xgBlock_t
379 : : type(xgBlock_t), intent(inout) :: X
380 : : type(xgBlock_t), intent(inout) :: Bm1X
381 : : end subroutine getBm1X
382 : : end interface
383 : :
384 : : ! Local scalars
385 : : integer :: tim_slice_me
386 : : integer :: nrows, ncols
387 : : integer :: space_res
388 : : !integer :: neigenpairs, bandpp, ndeg_filter
389 : : integer :: neigenpairs_ext
390 : : integer :: p_offset, ndeg_filter_max
391 : : !integer :: ierr
392 : : !real(dp) :: theta
393 : : real(dp) :: safe, tol ! for slice selection window
394 : : !logical :: is_lowpass, on_host, on_device ! todo use slice%..
395 : : ! ----> variables for logic of slice solver
396 : : type(matrixInfo_t) :: matrixInfo
397 : : type(xg_t) :: DivResults
398 : 18 : type(activeTask_t) :: task
399 : 18 : type(taskScheduler_t) :: scheduler
400 : 18 : type(asyncMemory_t) :: asyncMemory
401 : : ! <-----
402 : : ! Arrays
403 : : !real(dp) :: tsec(2)
404 : 18 : real(dp), allocatable :: moments(:)
405 : 18 : logical, allocatable, target :: mapper(:,:)
406 : :
407 : : ! *********************************************************************
408 : :
409 : : tim_slice_me = tim_slice_X
410 : : select case(task%me_id_slice)
411 : : case(1)
412 : : tim_slice_me = tim_slice_1
413 : : case(2)
414 : : tim_slice_me = tim_slice_2
415 : : case(3)
416 : : tim_slice_me = tim_slice_3
417 : : end select
418 : :
419 : 18 : nrows = slice%spacedim
420 : 18 : ncols = slice%neigenpairs
421 : : if (slice%paral_kgb==1) then ! colsrows distribution
422 : : nrows = slice%total_spacedim
423 : : ncols = slice%bandpp
424 : : end if
425 : :
426 : 18 : if (slice%space==SPACE_C) then
427 : 18 : space_res = SPACE_C
428 [ # # ]: 0 : else if (slice%space==SPACE_CR) then
429 : 0 : space_res = SPACE_R
430 : : else
431 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
432 : : end if
433 : :
434 : 18 : call xg_init(DivResults, space_res, slice%neigenpairs, 1, gpu_option=slice%gpu_option)
435 : :
436 : : call init_matrixInfo(matrixInfo, slice%comm_rows, slice%comm_cols, slice%spacecom, &
437 : : slice%neigenpairs,slice%total_spacedim, slice%spacedim, slice%space, &
438 : : slice%gpu_kokkos_nthrd, slice%gpu_thread_limit,slice%gpu_option, slice%paral_kgb, &
439 : 18 : slice%me_g0, slice%me_g0_fft, nspinor)
440 : :
441 : 18 : write(std_out,*) 'in slice_run'
442 : 18 : write(std_out,*) 'array of eigenvalues', rows(eigen), cols(eigen)
443 : 18 : flush(std_out)
444 : :
445 : : ! ============================= Prepare spectral slices =====================================
446 : :
447 : : !call timab(tim_slice_sched,1,tsec) ! fixme those timers
448 : : !ABI_NVTX_START_RANGE(NVTX_SLICE_SCHEDULE) ! fixme those markers
449 : :
450 : 18 : ndeg_filter_max = 10 ! todo IML autotune based on the presence of oscillations
451 [ - + ]: 18 : ABI_MALLOC(moments, (ndeg_filter_max+1))
452 : :
453 : : !ABI_NVTX_START_RANGE(NVTX_SLICE_RRQ)
454 : : call slice_getSpectralInfo(slice, X0, matrixInfo, moments, DivResults%self, residu, &
455 : 18 : getAX_BX, getBm1X)
456 : : !ABI_NVTX_END_RANGE()
457 : :
458 : : !ABI_NVTX_END_RANGE()
459 : : !call timab(tim_slice_sched,2,tsec)
460 : :
461 : 18 : call slice_splitSpectrum(slice, moments)
462 [ - + ]: 18 : ABI_FREE(moments)
463 : :
464 : : ! Attribute column vectors of X0 to slices using query results
465 [ - + + - : 72 : ABI_MALLOC(mapper, (slice%neigenpairs, slice%nslice))
- + + - +
- + - - +
- + - + ]
466 [ + + + + ]: 6966 : mapper = .false.
467 : 18 : write(std_out,*) 'tolerance in residual window=', slice%tolerance
468 : 18 : flush(std_out)
469 : :
470 : 18 : safe = 2.d0
471 : 18 : tol = 1e-2 ! slice%tolerance
472 : 18 : call slice_applySelectionWindow(slice, mapper, DivResults%self, residu, tol, safe)
473 : :
474 : 18 : write(std_out,*) 'created the following map slices to cols='
475 : 18 : write(std_out,*) mapper(1,:)
476 : 18 : flush(std_out)
477 : :
478 : : ! Include m+p where p offset 10%
479 : 18 : p_offset = 0.1d0
480 : 18 : slice%neigen_per_slice = ceiling(slice%neigen_per_slice*(1.d0+p_offset))
481 [ + + ]: 54 : neigenpairs_ext = sum(slice%neigen_per_slice)
482 : :
483 : 18 : write(std_out,*) 'allocating async memory of size', neigenpairs_ext
484 : 18 : flush(std_out)
485 : :
486 : : ! Allocate extended buffer in Linalg representation. Notice spacecom communicator (global)
487 : 18 : call slice_task_allocAsyncMemory(asyncMemory, matrixInfo, neigenpairs_ext)
488 : :
489 : : ! Init extended buffer with vectors from X or random vectors. This part also uses global comm.
490 : 18 : call slice_task_copyToAsyncMemory(asyncMemory, X0, matrixInfo, mapper, slice%neigen_per_slice)
491 : :
492 : 18 : write(std_out,*) 'getid after init', xgBlock_getid(asyncMemory%XextLinalg)
493 : 18 : flush(std_out)
494 : :
495 : : ! =========================== MPI distribution for slices ==================================
496 : :
497 : : ! MPI phase I: Compute 'process-to-slice(task)' distribution according to paral options
498 : : call slice_task_initSchedule(scheduler, slice%nslice, slice%nproc, slice%neigen_per_slice,&
499 : 18 : slice%paral_kgb, slice%paral_slice)
500 : : ! /IML\ in here there should be some per-rank parameter called slice%neigen_per_slice
501 : : ! that is updated along slice iterations
502 : : ! also nslice should be rank-wise as well
503 : :
504 : 18 : call slice_task_printSchedule(scheduler, std_out)
505 : :
506 : : ! MPI phase II: Compute 'slice columns(load)-to-process' distribution for next task
507 : : ! and apply schedule to distribute the async memory to processes
508 : 18 : call slice_task_initNextTask(scheduler, task, matrixInfo)
509 : 18 : call slice_task_enableAsync(asyncMemory, task, matrixInfo)
510 : :
511 : : !!! ^------- start region containing data on slice only --------^
512 : :
513 : : ! ============================ Active task execution =======================================
514 : :
515 : : ! IML is here
516 : : ! pass explicit slice parameters such as degree, lowb, uppb
517 : : ! call slice_task_runActiveTask(task...)
518 : :
519 : : ! if (paral_slice==yes) then
520 : : !! exit
521 : : !! otherwise, update slice%neigen_per_slice etc
522 : :
523 : : !
524 : : ! call timab(tim_slice_me,1,tsec)
525 : : !
526 : : ! if (slice%paral_kgb==0 .or. slice%paral_slice==DISABLE_PARAL) then
527 : : ! ! execute active tasks sequentially
528 : : !
529 : : ! do islice=1, nslice
530 : : ! if (task%active(islice)) then ! <--- normally it should not be that different..
531 : : ! call schedule_next_task(schedule, neigenpairs)
532 : : ! call init_active_memory(asyncMemory, task, X0, p)
533 : : ! call execute_active_task(task)
534 : : ! call mask_active_task(task, tol) ! mask extendedMem
535 : : ! end do
536 : : !
537 : : ! else
538 : : ! ! execute active tasks in parallel
539 : : !
540 : : ! call init_active_memory(asyncMemory, task, X0, p)
541 : : ! call execute_active_task(task)
542 : : ! call mask_active_task(task, tol) ! mask extendedMem
543 : : !
544 : : ! end if
545 : : !
546 : : ! ! Copy to spectrum memory only when active tasks have finished
547 : : ! call compress_extended_memory(work) ! extendendMem -> spectrumMem
548 : : !
549 : : !
550 : : ! call free_schedule(schedule)
551 : : !
552 : : !
553 : : !
554 : : ! ! Free temporary memory
555 : : ! call chebfi_free(chebfi)
556 : : ! ABI_SFREE(nrowsLinalg)
557 : : !
558 : : ! ! Timer is BEFORE the barrier !!
559 : : ! call timab(tim_slice_me,2,tsec)
560 : : !
561 [ - + ]: 18 : ABI_FREE(mapper)
562 : 18 : call xg_free(DivResults)
563 : 18 : call slice_task_freeAsyncMemory(asyncMemory)
564 : 18 : call slice_task_freeSchedule(scheduler)
565 : 18 : call slice_task_freeActiveTask(task)
566 : :
567 [ + - - + : 36 : end subroutine slice_run
- + - + -
+ - + - +
- + - + -
+ ]
568 : : !!***
569 : :
570 : : !----------------------------------------------------------------------
571 : :
572 : : !!****f* m_slice/slice_getSpectralInfo
573 : : !! NAME
574 : : !! slice_getSpectralInfo
575 : : !!
576 : : !! FUNCTION
577 : : !! Approximate eigenvalue spectrum and query input X for Ritz values and residuals.
578 : : !! This is a Hamiltonian-application-based approximation.
579 : : !! Inner calculations use colsrows representation to apply the Hamiltonian.
580 : : !!
581 : : !! INPUTS
582 : : !! X =eigenvector guess
583 : : !! getAX_BX =Hamiltonian application
584 : : !!
585 : : !! SOURCE
586 : :
587 [ - + ]: 18 : subroutine slice_getSpectralInfo(slice, X, matrixInfo, moments, eigen, resid, getAX_BX, getBm1X)
588 : :
589 : : implicit none
590 : :
591 : : ! Arguments
592 : : type(slice_t), intent(inout) :: slice
593 : : type(xgBlock_t), intent(inout) :: X
594 : : type(matrixInfo_t), intent(in) :: matrixInfo
595 : : real(dp), intent(inout) :: moments(:)
596 : : type(xgBlock_t), intent(inout) :: eigen
597 : : type(xgBlock_t), intent(inout) :: resid
598 : : interface
599 : : subroutine getAX_BX(X,AX,BX)
600 : : use m_xg, only : xgBlock_t
601 : : type(xgBlock_t), intent(inout) :: X
602 : : type(xgBlock_t), intent(inout) :: AX
603 : : type(xgBlock_t), intent(inout) :: BX
604 : : end subroutine getAX_BX
605 : : end interface
606 : : interface
607 : : subroutine getBm1X(X,Bm1X)
608 : : use m_xg, only : xgBlock_t
609 : : type(xgBlock_t), intent(inout) :: X
610 : : type(xgBlock_t), intent(inout) :: Bm1X
611 : : end subroutine getBm1X
612 : : end interface
613 : :
614 : : ! Local variables
615 : : ! Scalars
616 : : integer :: spacedim, tot_spacedim, nspinor
617 : : integer :: ndeg_filter_max, neigenpairs
618 : : integer :: ierr
619 : : !integer :: ideg
620 : : integer :: kmax
621 : : integer :: m_probe
622 : : integer :: my_shift, my_rank
623 : : integer :: space, spacecom, gpu_option
624 : : !integer :: k_sketch
625 : : real(dp) :: lambda_min, res_norm
626 : : !real(dp) :: center, radius
627 : : real(dp) :: lowb!, mineig, maxeig
628 : : real(dp) :: lanczos_lowb, lanczos_lowb_global
629 : : ! Derived types
630 : : #ifdef HAVE_OPENMP_OFFLOAD
631 : : integer :: me_g0
632 : : type(xg_t) :: W_dummy
633 : : integer :: work_size
634 : : #endif
635 : : !type(xg_t) :: BX
636 : : type(xgBlock_t) :: xXColsRows
637 : : type(xgBlock_t) :: eigen_me, resid_me
638 : 18 : type(xgTransposer_t) :: xgTransposerX
639 : : ! Arrays
640 : : #ifdef HAVE_OPENMP_OFFLOAD
641 : : real(dp) :: tsec(2)
642 : : #endif
643 : : ! *********************************************************************
644 : :
645 : 18 : neigenpairs = slice%neigenpairs
646 : 18 : tot_spacedim = slice%total_spacedim
647 : 18 : spacedim = slice%spacedim
648 : 18 : gpu_option = slice%gpu_option
649 : 18 : spacecom = slice%spacecom
650 : 18 : space = slice%space
651 : 18 : nspinor = matrixInfo%nspinor
652 : : #ifdef HAVE_OPENMP_OFFLOAD
653 : : me_g0 = slice%me_g0
654 : : if (slice%paral_kgb==1) then
655 : : me_g0 = slice%me_g0_fft
656 : : end if
657 : : #endif
658 : :
659 : : ! ============== Transpose ==============
660 : 18 : if (slice%paral_kgb==1) then
661 : :
662 : : ! Allocate memory for X in colsrows representation
663 : : call xgTransposer_constructor(xgTransposerX, X, xXColsRows, nspinor, STATE_LINALG,&
664 : : TRANS_ALL2ALL, slice%comm_rows, slice%comm_cols, 0, 0, slice%me_g0_fft,&
665 : 16 : gpu_option=slice%gpu_option, gpu_thread_limit=slice%gpu_thread_limit)
666 : :
667 : 16 : xgTransposerX%gpu_kokkos_nthrd = slice%gpu_kokkos_nthrd
668 : :
669 : : !call xmpi_barrier(slice%spacecom)
670 : : ABI_NVTX_START_RANGE(NVTX_SLICE_TRANSPOSE)
671 : 16 : call xgTransposer_transpose(xgTransposerX, STATE_COLSROWS)
672 : : ABI_NVTX_END_RANGE()
673 : :
674 : : else
675 : :
676 : : ! Use colsrows notion instead of X notion
677 : 2 : call xgBlock_setBlock(X, xXColsRows, spacedim, neigenpairs)
678 : :
679 : : end if
680 : :
681 : : ! Prevent invovl error
682 : : ! dummy invovl calculation to allocate buffers of full size
683 : : ! Workaround: getBm1X error related to invovl allocated buffers.
684 : : ! With this setup we allocate large then use 1 in Lanczos.
685 : : ! The inverse is not possible with current implementation. This order allows to avoid
686 : : ! make_invovl for 1 vector then make_invovl for nband vectors.
687 : : #ifdef HAVE_OPENMP_OFFLOAD
688 : : if (slice%paw) then
689 : : work_size = slice%neigenpairs ! fixme will complain
690 : : call xg_init(W_dummy, slice%space, tot_spacedim, work_size, xmpi_comm_null, &
691 : : me_g0=me_g0, gpu_option=slice%gpu_option)
692 : : call timab(tim_invovl, 1, tsec)
693 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_BM1X)
694 : : call getBm1X(W_dummy%self, W_dummy%self)
695 : : ABI_NVTX_END_RANGE()
696 : : call timab(tim_invovl, 2, tsec)
697 : : call xg_free(W_dummy)
698 : : end if
699 : : #endif
700 : :
701 : 18 : kmax = 30
702 : 18 : call computeBLanczos(matrixInfo, slice%paw, getAX_BX, getBm1X, kmax, lambda_min, res_norm)
703 : :
704 : 18 : lanczos_lowb = lambda_min - res_norm
705 : 18 : call xmpi_min(lanczos_lowb,lanczos_lowb_global,slice%spacecom,ierr)
706 : 18 : lowb = lanczos_lowb_global
707 : :
708 : 18 : write(std_out,*) 'Lanczos lambda_min=', lambda_min
709 : 18 : write(std_out,*) 'Lanczos res_norm =', res_norm
710 : 18 : write(std_out,*) 'Lanczos guarantee =', lanczos_lowb
711 : 18 : write(std_out,*) 'Lanczos guarantee(global) =', lanczos_lowb_global
712 : 18 : flush(std_out)
713 : :
714 : : ! Perform sensitivity study of trace estimation for these parameters
715 : : ! Keep m_probe small allows to reduce noise
716 : 18 : m_probe = 5
717 : 18 : ndeg_filter_max = size(moments) - 1
718 : 18 : write(std_out,*) 'Here I compute Stochastic Trace Estimation'
719 : 18 : write(std_out,*) 'm_probe= ', m_probe
720 : 18 : write(std_out,*) 'ndeg_filter=', ndeg_filter_max
721 : 18 : flush(std_out)
722 : :
723 : : call computeTraceEstimation(matrixInfo, slice%ecut, slice%paw, slice%tolerance, &
724 : 18 : getAX_BX, getBm1X, ndeg_filter_max, m_probe, lanczos_lowb_global, moments)
725 : :
726 : : ! Define working spectrum to be splitted to slices
727 : 18 : slice%mineig = lanczos_lowb
728 : 18 : slice%mineig_global = lanczos_lowb_global
729 : 18 : slice%maxeig_global = slice%ecut
730 : :
731 : 18 : write(std_out,*) 'STE exited'
732 : 18 : flush(std_out)
733 : :
734 : : ! Compute Rayleigh quotients (colsrows distribution) ! fixme timers
735 : : !ABI_NVTX_START_RANGE(NVTX_CHEBFI2_RRQ)
736 : : !call timab(tim_RR_q, 1, tsec)
737 : :
738 [ + + ]: 18 : if (slice%paral_kgb==1) then
739 : 16 : my_rank = xmpi_comm_rank(slice%spacecom)
740 : 16 : my_shift = my_rank * slice%bandpp
741 : 16 : call xgBlock_reshape(eigen, 1, slice%neigenpairs)
742 : 16 : call xgBlock_reshape(resid, 1, slice%neigenpairs)
743 : 16 : call xgBlock_setBlock(eigen, eigen_me, 1, slice%bandpp, fcol=my_shift+1)
744 : 16 : call xgBlock_setBlock(resid, resid_me, 1, slice%bandpp, fcol=my_shift+1)
745 : 16 : call xgBlock_reshape(eigen, slice%neigenpairs, 1)
746 : 16 : call xgBlock_reshape(resid, slice%neigenpairs, 1)
747 : 16 : call xgBlock_reshape(eigen_me, slice%bandpp, 1)
748 : 16 : call xgBlock_reshape(resid_me, slice%bandpp, 1)
749 : : else
750 : 2 : call xgBlock_setBlock(eigen, eigen_me, rows(eigen), 1)
751 : 2 : call xgBlock_setBlock(resid, resid_me, rows(resid), 1)
752 : : end if
753 : :
754 : 18 : call slice_queryCandidates(slice, xXColsRows, eigen_me, resid_me, getAX_BX)
755 : : !call timab(tim_RR_q, 2, tsec)
756 : : !ABI_NVTX_END_RANGE()
757 : :
758 : 18 : call xgBlock_mpi_sum(eigen, comm=slice%spacecom)
759 : 18 : call xgBlock_mpi_sum(resid, comm=slice%spacecom)
760 : :
761 : : ! ============== Transpose ==============
762 [ + + ]: 18 : if (slice%paral_kgb == 1) then
763 : 16 : call xmpi_barrier(slice%spacecom)
764 : : ABI_NVTX_START_RANGE(NVTX_SLICE_TRANSPOSE)
765 : 16 : call xgTransposer_transpose(xgTransposerX, STATE_LINALG)
766 : : ABI_NVTX_END_RANGE()
767 : :
768 : : ! reset buffers to right address
769 [ - + ]: 16 : if (xmpi_comm_size(slice%spacecom) == 1) then
770 : 0 : call xgBlock_setBlock(xXColsRows, X, spacedim, neigenpairs)
771 : : end if
772 : : else
773 : 2 : call xgBlock_setBlock(xXColsRows, X, spacedim, neigenpairs)
774 : : end if
775 : :
776 [ + + ]: 18 : if (slice%paral_kgb == 1) then
777 : 16 : call xgTransposer_free(xgTransposerX)
778 : : end if
779 : :
780 [ + + - + ]: 36 : end subroutine slice_getSpectralInfo
781 : : !!***
782 : :
783 : : !----------------------------------------------------------------------
784 : :
785 : : !!****f* m_slice/slice_queryCandidates
786 : : !! NAME
787 : : !! slice_queryCandidates
788 : : !!
789 : : !! FUNCTION
790 : : !! Compute Rayleigh-Ritz quotients and residuals in colsrows MPI representation.
791 : : !! Return bandpp indicators per MPI process.
792 : : !! These two indicators are useful to rank candidates per slices.
793 : : !!
794 : : !! SOURCE
795 : :
796 : 18 : subroutine slice_queryCandidates(slice, xXColsRows, eigen, resid, getAX_BX)
797 : :
798 : : implicit none
799 : :
800 : : !Arguments ------------------------------------
801 : : type(slice_t), intent(inout) :: slice
802 : : type(xgBlock_t), intent(inout) :: xXColsRows
803 : : type(xgBlock_t), intent(inout) :: eigen
804 : : type(xgBlock_t), intent(inout) :: resid
805 : : interface
806 : : subroutine getAX_BX(X,AX,BX)
807 : : use m_xg, only : xgBlock_t
808 : : type(xgBlock_t), intent(inout) :: X
809 : : type(xgBlock_t), intent(inout) :: AX
810 : : type(xgBlock_t), intent(inout) :: BX
811 : : end subroutine getAX_BX
812 : : end interface
813 : :
814 : : !Local variables-------------------------------
815 : : type(xg_t) :: XAB
816 : : type(xg_t) :: Results1
817 : : type(xg_t) :: Results2
818 : : type(xg_t) :: norml
819 : : type(xgBlock_t) :: xAXColsRows
820 : : type(xgBlock_t) :: xBXColsRows
821 : : type(xgBlock_t) :: xRXColsRows
822 : : integer :: space_res
823 : : integer :: me_g0, nrows, ncols, gpu_option
824 : : real(dp) :: tsec(2)
825 : :
826 : : ! *********************************************************************
827 : :
828 : 18 : if (slice%space==SPACE_C) then
829 : 18 : space_res = SPACE_C
830 [ # # ]: 0 : else if (slice%space==SPACE_CR) then
831 : 0 : space_res = SPACE_R
832 : : else
833 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
834 : : end if
835 : :
836 : 18 : ncols = slice%bandpp
837 : 18 : nrows = rows(xXColsRows)
838 : 18 : me_g0 = slice%me_g0
839 : 18 : gpu_option = slice%gpu_option
840 [ + + ]: 18 : if (slice%paral_kgb==1) then
841 : : ncols = slice%bandpp
842 : 16 : nrows = slice%total_spacedim
843 : 16 : me_g0 = slice%me_g0_fft
844 : : end if
845 : 18 : call xg_init(Results1, space_res, ncols, 1, gpu_option=slice%gpu_option)
846 : 18 : call xg_init(Results2, space_res, ncols, 1, gpu_option=slice%gpu_option)
847 : 18 : call xg_init(norml, SPACE_R, slice%bandpp, 1)
848 : :
849 : 18 : call xg_init(XAB,slice%space,nrows,3*ncols,slice%spacecom,me_g0=me_g0,gpu_option=gpu_option)
850 : 18 : call xg_setBlock(XAB, xAXColsRows, nrows, ncols)
851 : 18 : call xg_setBlock(XAB, xBXColsRows, nrows, ncols, fcol=ncols+1)
852 : 18 : call xg_setBlock(XAB, xRXColsRows, nrows, ncols, fcol=2*ncols+1)
853 : :
854 : : ! Compute A*Psi
855 : 18 : call timab(tim_getAX_BX,1,tsec)
856 : : ABI_NVTX_START_RANGE(NVTX_SLICE_GET_AX_BX)
857 : 18 : call getAX_BX(xXColsRows, xAXColsRows, xBXColsRows)
858 : 18 : call xgBlock_zero_im_g0(xAXColsRows)
859 : 18 : call xgBlock_zero_im_g0(xBXColsRows)
860 : : ABI_NVTX_END_RANGE()
861 : 18 : call timab(tim_getAX_BX,2,tsec)
862 : :
863 : : ! <Psi|H|Psi>
864 : 18 : call xgBlock_colwiseDotProduct(xXColsRows, xAXColsRows, Results1%self, comm_loc=xmpi_comm_null)
865 : :
866 : : ! <Psi|S|Psi>
867 : 18 : call xgBlock_colwiseDotProduct(xXColsRows, xBXColsRows, Results2%self, comm_loc=xmpi_comm_null)
868 : :
869 : : ! eig = <Psi|H|Psi> / <Psi|S|Psi>
870 : 18 : call xgBlock_colwiseDivision(Results1%self, Results2%self, eigen)
871 : :
872 : : ! xRXColsRows = S|Psi>
873 : 18 : call xgBlock_copy(xBXColsRows,xRXColsRows)
874 : : ! xRXColsRows = - eigen * S|Psi>
875 : 18 : call xgBlock_ymax(xRXColsRows,eigen,0,1)
876 : : ! norml = |eigen*S|Psi>|^2
877 : 18 : call xgBlock_colwiseNorm2(xBXColsRows, norml%self, comm_loc=xmpi_comm_null)
878 : : ! xRXColsRows = H|Psi> - eigen * S|Psi>
879 : 18 : call xgBlock_add(xRXColsRows,xAXColsRows)
880 : : ! resid = |xRXColsRows|^2
881 : 18 : call xgBlock_colwiseNorm2(xRXColsRows, resid, comm_loc=xmpi_comm_null)
882 : :
883 : : ! compute relative resid=resid/(|eigen|*|BX|)^2
884 : : ! scales with the eigenvalue
885 : 18 : call xgBlock_colwiseDivision(resid, norml%self, resid)
886 : :
887 : 18 : call xg_free(Results1)
888 : 18 : call xg_free(Results2)
889 : 18 : call xg_free(XAB)
890 : 18 : call xg_free(norml)
891 : :
892 [ + - ]: 18 : end subroutine slice_queryCandidates
893 : : !!***
894 : :
895 : : !----------------------------------------------------------------------
896 : :
897 : : !!****f* m_slice/slice_unitTest
898 : : !! NAME
899 : : !! slice_unitTest
900 : : !!
901 : : !! FUNCTION
902 : : !! Given a set of vectors X computes the columnwise norm squared, on:
903 : : !! - CPU sequential,
904 : : !! - CPU MPI,
905 : : !! - GPU sequential (if GPU available),
906 : : !! - GPU MPI (if GPU available).
907 : : !! Checks that all four are the same.
908 : : !!
909 : : !! SOURCE
910 : :
911 : : function slice_unitTest(X) result(ierr)
912 : :
913 : : implicit none
914 : :
915 : : !Arguments ------------------------------------
916 : : type(xgBlock_t), intent(in) :: X
917 : : integer :: ierr
918 : :
919 : : !Local variables-------------------------------
920 : : integer :: nrows,ncols
921 : : integer :: comm
922 : : integer :: gpu_option
923 : : real(dp) :: id_cpu, id_from_gpu, id_to_gpu
924 : :
925 : : ! *********************************************************************
926 : :
927 : : ! Get nrows, nrcols, gpu_option of X
928 : : call xgBlock_getSize(X,nrows,ncols)
929 : : call xgBlock_get_gpu_option(X,gpu_option)
930 : : call xgBlock_get_communicator(X,comm)
931 : : ! above this is just comm(X)
932 : : write(std_out,*) '<--*--> Start unitary test for matrix'
933 : : write(std_out,*) 'nrows =', nrows
934 : : write(std_out,*) 'ncols =', ncols
935 : : write(std_out,*) 'gpu_option =', gpu_option
936 : : write(std_out,*) 'communicator=', comm
937 : :
938 : : id_cpu = -1.0d0
939 : : id_from_gpu = -1.0d0
940 : : id_to_gpu = -1.0d0
941 : :
942 : : if (gpu_option==ABI_GPU_DISABLED) then
943 : : id_cpu = xgBlock_getid(X, comm)
944 : : write(std_out,*) 'id_cpu =', id_cpu
945 : : end if
946 : :
947 : : #if defined HAVE_GPU && defined HAVE_OPENMP_OFFLOAD
948 : : if (gpu_option/=ABI_GPU_DISABLED) then
949 : : call xgBlock_copy_from_gpu(X)
950 : : id_from_gpu = xgBlock_getid(X, comm)
951 : : call xgBlock_copy_to_gpu(X)
952 : : write(std_out,*) 'id_from_gpu=', id_from_gpu
953 : : else
954 : : call xgBlock_copy_to_gpu(X)
955 : : call xgBlock_copy_from_gpu(X)
956 : : id_to_gpu = xgBlock_getid(X, comm)
957 : : call xgBlock_copy_to_gpu(X)
958 : : call xgBlock_copy_from_gpu(X)
959 : : write(std_out,*) 'id_to_gpu =', id_to_gpu
960 : : end if
961 : : #endif
962 : :
963 : : ! TODO getid on xmpi_comm_null
964 : : ! then sum across comm
965 : : ! should give the same as
966 : : ! getid on comm
967 : : ! this is useful for MPI
968 : :
969 : : ierr = 0
970 : : if (id_cpu < 0) ierr = -1
971 : : if (id_from_gpu < 0) ierr = -1
972 : : if (id_to_gpu < 0) ierr = -1
973 : :
974 : : end function slice_unitTest
975 : : !!***
976 : :
977 : : !----------------------------------------------------------------------
978 : :
979 : : !!****f* m_slice/slice_splitSpectrum
980 : : !! NAME
981 : : !! slice_splitSpectrum
982 : : !!
983 : : !! FUNCTION
984 : : !! Split working spectrum [a,b) mapped to [-1,1) with center and radius.
985 : : !!
986 : : !! INPUT
987 : : !! a= slice%mineig_global
988 : : !! b= slice%maxeig_global
989 : : !!
990 : : !! SIDE EFFECTS
991 : : !! slice%low_bounds
992 : : !! slice%upp_bounds
993 : : !! slice%neigen_per_slice
994 : : !!
995 : : !! SOURCE
996 : :
997 [ - + ]: 18 : subroutine slice_splitSpectrum(slice, moments)
998 : :
999 : : implicit none
1000 : :
1001 : : ! Arguments ------------------------------------
1002 : : type(slice_t), intent(inout) :: slice
1003 : : real(dp), intent(inout) :: moments(:)
1004 : :
1005 : : ! Local variables --------------------------------
1006 : : integer :: uppb_loc, ib
1007 : : integer :: ngrid_fine, ngrid_coarse
1008 : : integer :: half_index, i_left, i_right
1009 : : integer :: i_split, num_moments
1010 : : logical :: gap_left, gap_right
1011 : : real(dp) :: a, b
1012 : : real(dp) :: center, radius
1013 : : real(dp) :: step_fine, step_coarse
1014 : : real(dp) :: partial_mass, half_mass
1015 : : real(dp) :: b_scaled
1016 : : real(dp) :: mass_left, mass_right
1017 : : real(dp) :: lambda_plus_wanted
1018 : 18 : real(dp), allocatable :: cumm_eigen_count(:)
1019 : 18 : real(dp), allocatable :: bgrid_fine(:)
1020 : 18 : real(dp), allocatable :: bgrid_coarse(:)
1021 : 18 : real(dp), allocatable :: work(:)
1022 : : character(len=500) :: msg
1023 : :
1024 : : ! *********************************************************************
1025 : :
1026 [ - + ]: 18 : if (slice%nslice/=2) then
1027 : 0 : ABI_ERROR("spectral bidirectional split not yet implemented for 3 slices")
1028 : : end if
1029 : :
1030 : 18 : ngrid_coarse = 10 ! coarse, just to find uppb, hardcoded
1031 : 18 : ngrid_fine = 30 ! used for cumulative eigenvalue count, hardcoded
1032 : :
1033 : 18 : a = slice%mineig_global
1034 : 18 : b = slice%maxeig_global
1035 : :
1036 : 18 : num_moments = size(moments)
1037 [ - + + - : 54 : ABI_MALLOC(work, (num_moments))
+ - + - -
+ - + -
+ ]
1038 [ - + ]: 18 : ABI_MALLOC(bgrid_coarse, (ngrid_coarse))
1039 [ - + ]: 18 : ABI_MALLOC(bgrid_fine, (ngrid_fine))
1040 [ - + ]: 18 : ABI_MALLOC(cumm_eigen_count, (ngrid_fine))
1041 : :
1042 : 18 : center = (a + b) / 2.d0
1043 : 18 : radius = (b - a) / 2.d0
1044 : :
1045 : : ! #########################################
1046 : : ! ########## Coarse resolution ############
1047 : : ! #########################################
1048 : :
1049 : : ! scan with lowpass
1050 : : ! first pass uppb is actually unknown
1051 : 18 : step_coarse = (b - a) / (ngrid_coarse - 1)
1052 [ + - + + : 594 : bgrid_coarse = (/ ( a + (ib-1)*step_coarse, ib=1,ngrid_coarse ) /)
+ - + + +
- - + - -
+ + ]
1053 : 18 : uppb_loc = -1
1054 [ + - ]: 36 : do ib=1, ngrid_coarse
1055 : 36 : b_scaled = (bgrid_coarse(ib) - center) / radius
1056 : 36 : partial_mass = get_eigenvalue_count(b_scaled, moments, work)
1057 : :
1058 : 36 : write(std_out,*) ib, 'scan: <=', bgrid_coarse(ib), 'mass=', partial_mass
1059 : 36 : flush(std_out)
1060 : :
1061 [ + + ]: 36 : if (partial_mass > slice%neigenpairs) then
1062 : : uppb_loc = ib
1063 : : exit
1064 : : end if
1065 : :
1066 : : end do
1067 : :
1068 : 18 : lambda_plus_wanted = (bgrid_coarse(uppb_loc-1) + bgrid_coarse(uppb_loc)) / 2.d0
1069 : 18 : write(std_out,*) 'found upp bound in', lambda_plus_wanted
1070 : 18 : write(std_out,*) 'estimated mass=', get_eigenvalue_count((lambda_plus_wanted-center)/radius, moments, work)
1071 : 18 : write(std_out,*) 'starting adaptive refinement ..'
1072 : 18 : flush(std_out)
1073 : :
1074 : : ! #########################################
1075 : : ! ########### Fine resolution #############
1076 : : ! #########################################
1077 : :
1078 : : ! Now compute eigenvalue count
1079 : 18 : step_fine = (lambda_plus_wanted - a) / (ngrid_fine - 1)
1080 [ + - + + : 1674 : bgrid_fine = (/ (a + (ib-1)*step_fine, ib=1,ngrid_fine) /)
+ - + + +
- - + - -
+ + ]
1081 [ + + ]: 558 : do ib=1, ngrid_fine
1082 : 540 : b_scaled = (bgrid_fine(ib) - center) / radius
1083 : 540 : partial_mass = get_eigenvalue_count(b_scaled, moments, work)
1084 : 540 : cumm_eigen_count(ib) = partial_mass
1085 [ - + - - ]: 540 : if (partial_mass < 0.d0 .and. abs(partial_mass) > 1.d0) then
1086 : 0 : write(msg,'(a)') "fine resolution for eigenvalue count failed due to oscillations. ",&
1087 : 0 : "Reduce ndeg_filter_max (hard-coded) as a solution"
1088 : 0 : ABI_ERROR(msg)
1089 : : end if
1090 : 558 : write(std_out,*) ib, 'scan: <=', bgrid_fine(ib), 'mass=', partial_mass
1091 : : end do
1092 : : !! todo if not found then add points in the fine grid...
1093 : :
1094 : : ! Prepare: detect gap existence in the interior of slice
1095 : : ! define value of smallest_gap
1096 : : ! todo
1097 : :
1098 : : ! step 1 cut in half balanced mass
1099 : 18 : half_mass = cumm_eigen_count(ngrid_fine)/2.d0
1100 [ + - + - : 576 : half_index = minloc(abs(half_mass - cumm_eigen_count), dim=1)
+ + + + ]
1101 : : ! step 2 adjust so that it is on constant mass (predicts gap)
1102 : : ! bidirectional search left and right
1103 : 18 : write(std_out,*) 'bidirectional search from index=', half_index
1104 : 18 : write(std_out,*) 'half mass=', half_mass
1105 : 18 : flush(std_out)
1106 : 18 : i_left = half_index
1107 : 18 : i_right = half_index
1108 : 18 : gap_left = .false.
1109 : 18 : gap_right = .false.
1110 [ + + - + ]: 36 : do while((.not.gap_right .and. .not.gap_left) .and. (i_left >= 2 .and. i_right <=ngrid_fine-1))
1111 : : ! this is if gap exists. If it does not exist.. must minimize using smallest_gap
1112 : 18 : gap_left = abs(cumm_eigen_count(i_left) - cumm_eigen_count(i_left-1)) < 1e-4
1113 : 18 : gap_right = abs(cumm_eigen_count(i_right) - cumm_eigen_count(i_right+1)) < 1e-4
1114 : 18 : i_left = i_left - 1
1115 : 18 : i_right = i_right + 1
1116 : : end do
1117 [ + - ]: 18 : if (gap_right) then
1118 : 18 : i_split = i_right
1119 : : end if
1120 [ - + ]: 18 : if (gap_left) then
1121 : 0 : i_split = i_left
1122 : : end if
1123 [ - + ]: 18 : if (.not.gap_right .and. .not.gap_left) then
1124 : : ! todo treat this case
1125 : 0 : ABI_ERROR("spectrum has no gap..")
1126 : : end if
1127 : 18 : mass_left = get_eigenvalue_count((bgrid_fine(i_split)-center)/radius,moments,work)
1128 : 18 : mass_right = slice%neigenpairs - mass_left
1129 : 18 : write(std_out,*) 'i_split val=', i_split, bgrid_fine(i_split)
1130 : 18 : write(std_out,*) 'mass left=', mass_left
1131 : 18 : write(std_out,*) 'mass right=', mass_right
1132 : 18 : flush(std_out)
1133 : :
1134 : : ! Store results into the slice
1135 [ + - ]: 18 : slice%neigen_per_slice(1) = ceiling(mass_left)
1136 [ + - ]: 18 : slice%neigen_per_slice(2) = ceiling(mass_right)
1137 : :
1138 : 18 : slice%low_bounds(1) = slice%mineig
1139 : 18 : slice%low_bounds(2) = bgrid_fine(i_split)
1140 : :
1141 : 18 : slice%upp_bounds(1) = bgrid_fine(i_split)
1142 : 18 : slice%upp_bounds(2) = bgrid_fine(ngrid_fine)
1143 : :
1144 : 18 : ABI_FREE(bgrid_fine)
1145 : 18 : ABI_FREE(bgrid_coarse)
1146 : 18 : ABI_FREE(cumm_eigen_count)
1147 [ - + ]: 18 : ABI_FREE(work)
1148 : :
1149 : 18 : end subroutine slice_splitSpectrum
1150 : : !!***
1151 : :
1152 : : !----------------------------------------------------------------------
1153 : :
1154 : : !!****f* m_slice/slice_applySelectionWindow
1155 : : !! NAME
1156 : : !! slice_applySelectionWindow
1157 : : !!
1158 : : !! FUNCTION
1159 : : !! Input: k=m+p ou m la masse et p oversampling parameter
1160 : : !! step 1) calculer les valeurs de Ritz
1161 : : !! step 2) calculer les residus
1162 : : !! step 3) mettre les valeurs de Ritz dans part_low_bounds, part_upp_bounds sous condition que res<tol
1163 : : !! step 4) sinon faire un sketch de taille K du reste
1164 : : !!
1165 : : !! INPUT
1166 : : !! safe is between 2 and 10
1167 : : !!
1168 : : !! OUTPUT
1169 : : !! mapper contains the column indices per slice (allowing repetitions)
1170 : : !!
1171 : : !! SOURCE
1172 : :
1173 [ - + ]: 18 : subroutine slice_applySelectionWindow(slice, mapper, eigen, resid, tol, safe)
1174 : :
1175 : : implicit none
1176 : :
1177 : : type(slice_t), intent(inout) :: slice
1178 : : logical, intent(inout) :: mapper(:,:)
1179 : : type(xgBlock_t), intent(inout) :: eigen
1180 : : type(xgBlock_t), intent(inout) :: resid
1181 : : real(dp), intent(in) :: tol
1182 : : real(dp), intent(in) :: safe
1183 : :
1184 : : integer :: nband, iband, islice, m, k_kept
1185 : : real(dp) :: uppb, lowb, uppb_plus, lowb_minus
1186 : : real(dp) :: theta, relaxf, relres
1187 : : real(dp), pointer :: resid_vals(:,:) => null()
1188 : : complex(dp), pointer :: thetas(:,:) => null()
1189 : :
1190 : : ! *********************************************************************
1191 : :
1192 : 18 : nband = slice%neigenpairs
1193 : :
1194 : 18 : call xgBlock_reverseMap(eigen, thetas, rows=nband, cols=1)
1195 : 18 : call xgBlock_reverseMap(resid, resid_vals, rows=nband, cols=1)
1196 : :
1197 : : ! apply selection window per slice
1198 [ + + ]: 54 : do islice=1,slice%nslice
1199 : 36 : m = slice%neigen_per_slice(islice)
1200 : 36 : lowb = slice%low_bounds(islice)
1201 : 36 : uppb = slice%upp_bounds(islice)
1202 : 36 : write(std_out,*) islice, 'slice bounds', lowb, uppb
1203 : :
1204 : : ! keep candidates in expanded interval
1205 : 36 : k_kept = 0
1206 [ + + ]: 6948 : do iband=1, nband
1207 : 6912 : relres = resid_vals(iband, 1)
1208 [ - + ]: 6948 : if (relres < tol) then
1209 : : ! Soft selection:
1210 : : ! relaxed bounds for converged eigenpairs near boundaries
1211 : 0 : relaxf = safe * relres
1212 : 0 : lowb_minus = lowb - relaxf
1213 : 0 : uppb_plus = uppb + relaxf
1214 : 0 : theta = real(thetas(iband, 1))
1215 : 0 : write(std_out,*) iband, 'theta=', theta, 'resid=', relres
1216 : 0 : flush(std_out)
1217 : :
1218 : : ! accept if lambda in relaxed interval
1219 [ # # # # ]: 0 : if (theta < uppb_plus .and. theta > lowb_minus) then
1220 : 0 : k_kept = k_kept + 1
1221 : 0 : mapper(iband, islice) = .true.
1222 : : end if
1223 : : end if
1224 : : end do
1225 : :
1226 : : ! todo m not set must modify computeTraceEstimation and separate the split spectrum from that
1227 : : ! because split spectrum will store info on slice object that is invisible from trace
1228 : 36 : write(std_out,*) 'selection window kept', k_kept, 'for eigendimension', m
1229 : 54 : flush(std_out)
1230 : : end do
1231 : :
1232 : 18 : end subroutine slice_applySelectionWindow
1233 : : !!***
1234 : :
1235 [ # # # # : 0 : end module m_slice
# # # # #
# # # # #
# # # # #
# ]
1236 : : !!***
|
