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 : !!***
|
