Branch data Line data Source code
1 : : !!****f* ABINIT/m_slice_cprj
2 : : !! NAME
3 : : !! m_slice_cprj
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the types and routines used to apply the
7 : : !! Spectrum Slicing method (2021 implementation using xG abstraction layer)
8 : : !! It mainly defines a 'slice' datatypes and associated methods.
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2023-2025 ABINIT group (LB,IML)
12 : : !! This file is distributed under the terms of the
13 : : !! gnu general public license, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !! for the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_slice_cprj
26 : :
27 : : use defs_basis
28 : : use defs_abitypes
29 : : use m_abicore
30 : : use m_errors
31 : : use m_time, only : timab
32 : : use m_sort, only: sort_dp
33 : :
34 : : use m_cgtools
35 : : use m_xg
36 : : use m_xgTransposer
37 : : use m_xg_ortho_RR
38 : : use m_xg_nonlop
39 : :
40 : : use m_polynomial_filter, only: buildChebyshevJacksonCoeffs
41 : : use m_trace_estimation, only: smallestTridiagEigenpair
42 : :
43 : : use m_xmpi
44 : : use m_xomp
45 : : #ifdef HAVE_OPENMP
46 : : use omp_lib
47 : : #endif
48 : :
49 : : implicit none
50 : :
51 : : private
52 : :
53 : : !Several (private) parameters
54 : : !-------------------------------------------------
55 : :
56 : : integer, parameter :: tim_init = 2171
57 : : integer, parameter :: tim_free = 2172
58 : : integer, parameter :: tim_cprj = 2173
59 : : integer, parameter :: tim_invovl = 2174
60 : : integer, parameter :: tim_residu = 2175
61 : : integer, parameter :: tim_RR = 2176
62 : : integer, parameter :: tim_transpose = 2177
63 : : integer, parameter :: tim_RR_q = 2178
64 : : integer, parameter :: tim_postinvovl = 2179
65 : : integer, parameter :: tim_swap = 2180
66 : : integer, parameter :: tim_amp_f = 2181
67 : : integer, parameter :: tim_barrier = 2182
68 : : integer, parameter :: tim_copy = 2183
69 : : integer, parameter :: tim_ax_k = 2184
70 : : integer, parameter :: tim_ax_v = 2185
71 : : integer, parameter :: tim_ax_nl = 2186
72 : : integer, parameter :: tim_enl = 2187
73 : : integer, parameter :: tim_ortho = 2188
74 : :
75 : : integer, parameter :: tim_slice1_fi = 2190
76 : : integer, parameter :: tim_slice1_rr = 2191
77 : : integer, parameter :: tim_slice1_pr = 2192
78 : : integer, parameter :: tim_slice2_fi = 2193
79 : : integer, parameter :: tim_slice2_rr = 2194
80 : : integer, parameter :: tim_slice2_pr = 2195
81 : :
82 : : !Public 'slice' datatype
83 : : !-------------------------------------------------
84 : :
85 : : type, public :: slice_t
86 : :
87 : : integer :: space
88 : : integer :: space_cprj
89 : : integer :: spacedim ! Space dimension for one vector
90 : : integer :: cprjdim ! cprj dimension
91 : : integer :: blockdim_cprj ! ncols of cprj for slice only
92 : : integer :: all_blockdim_cprj ! ncols of cprj for all vectors
93 : : integer :: neigenpairs ! Number of eigen values/vectors we want
94 : : integer :: ndeg_filter ! Degree of the polynomial filter as input
95 : : integer :: nslice ! Number of spectral slices
96 : : integer :: slicedim ! Number of eigen values/vectors we want in one slice
97 : : integer :: paral_slice ! slice parallelization strategy (off)
98 : : integer :: spectral_cut ! how to cut the eigenvalue spectral interval (off)
99 : : integer :: nbdbuf ! Number of bands in the buffer
100 : : real(dp) :: tolerance ! Tolerance on the residu to stop the minimization
101 : : real(dp) :: ecut ! Ecut for Chebfi oracle
102 : : real(dp) :: tolfilter ! Polynomial filter wanted amplification
103 : :
104 : : ! MPI-related
105 : : integer :: bandpp
106 : : integer :: total_spacedim ! Maybe not needed
107 : : integer :: spacecom ! Communicator for MPI
108 : :
109 : : logical :: paw
110 : : integer :: eigenProblem !1 (A*x = (lambda)*B*x), 2 (A*B*x = (lambda)*x), 3 (B*A*x = (lambda)*x)
111 : : integer :: me_g0
112 : : integer :: me_g0_fft
113 : :
114 : : !ARRAYS for entire spectum
115 : : type(xgBlock_t) :: AllX ! Block of initial and final solution
116 : : type(xg_t) :: AllAX ! space to save AX Hamiltonian application
117 : :
118 : : ! cprj for entire spectrum
119 : : type(xgBlock_t) :: AllcprjX
120 : : type(xg_t) :: Allcprj_work
121 : :
122 : : !ARRAYS on slice using slice-safe memory
123 : : type(xg_t) :: X_SLICE ! memory independent of the entire spectrum
124 : : type(xgBlock_t) :: X ! pointers to slice memory
125 : : type(xgBlock_t) :: AX !
126 : :
127 : : ! space to hold X_next, X_prev for Chebyshev recursion, slice only
128 : : type(xg_t) :: X_NP
129 : : type(xgBlock_t) :: X_next
130 : : type(xgBlock_t) :: X_prev
131 : : !SWAP POINTERS on slice, also for Chebyshev recursion
132 : : type(xgBlock_t) :: X_swap
133 : : type(xgBlock_t) :: AX_swap
134 : :
135 : : ! cprj for slice only
136 : : type(xgBlock_t) :: cprjX ! pointer to AllcprjX
137 : : type(xgBlock_t) :: cprj_work ! pointer to Allcprj_work
138 : : type(xg_t) :: cprj_work2
139 : :
140 : : ! Independent of number of vectors, common to slice and spectrum
141 : : type(xg_t) :: proj_work
142 : : type(xg_nonlop_t) :: xg_nonlop
143 : :
144 : : type(xgBlock_t) :: eigenvalues
145 : :
146 : : end type slice_t
147 : :
148 : : !Public methods associated to 'slice' datatype
149 : : !-------------------------------------------------
150 : : public :: slice_init
151 : : public :: slice_free
152 : : public :: slice_memInfo
153 : : public :: slice_run_cprj
154 : :
155 : : CONTAINS !========================================================================================
156 : : !!***
157 : :
158 : : !!****f* m_slice_cprj/slice_init
159 : : !! NAME
160 : : !! slice_init
161 : : !!
162 : : !! FUNCTION
163 : : !! Initialize a 'slice' datastructure.
164 : : !!
165 : : !! INPUTS
166 : : !! bandpp= number of 'bands' handled by a processor
167 : : !! eigenProblem= type of eigenpb: 1 (A*x = (lambda)*B*x), 2 (A*B*x = (lambda)*x), 3 (B*A*x = (lambda)*x)
168 : : !! me_g0= 1 if this processors treats G=0, 0 otherwise
169 : : !! neigenpairs= number of requested eigenvectors/eigenvalues
170 : : !! ndeg_filter= polynomial degree of the Chebyshev filter (.i.e. number of H applications)
171 : : !! space= defines in which space we are (columns, rows, etc.)
172 : : !! spacecom= MPI communicator
173 : : !! spacedim= space dimension for one vector
174 : : !! paw= flag. TRUE if current calculation ses the PAW approach
175 : : !! ecut= plane-wave cut-off energy
176 : : !! tolerance= tolerance criterion on the residu to stop the minimization
177 : : !! nslice= number of slices
178 : : !! tolfilter= tolerance for filter amplification (a priori)
179 : : !! paral_slice= slice parallelization
180 : : !! spectral_cut= how to cut the eigenvalue spectrum (a priori)
181 : : !!
182 : : !! OUTPUT
183 : : !!
184 : : !! SIDE EFFECTS
185 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
186 : : !!
187 : : !! SOURCE
188 : :
189 : 0 : subroutine slice_init(slice,neigenpairs,spacedim,cprjdim,tolerance,ecut,bandpp, &
190 : : ndeg_filter,nbdbuf,space,space_cprj,eigenProblem,spacecom,me_g0,paw,&
191 : : nslice,tolfilter,paral_slice,spectral_cut,xg_nonlop,me_g0_fft)
192 : :
193 : : !Arguments ------------------------------------
194 : : integer , intent(in ) :: bandpp
195 : : integer , intent(in ) :: eigenProblem
196 : : integer , intent(in ) :: me_g0
197 : : integer , intent(in ) :: me_g0_fft
198 : : integer , intent(in ) :: neigenpairs
199 : : integer , intent(in ) :: ndeg_filter
200 : : integer , intent(in ) :: nbdbuf
201 : : integer , intent(in ) :: space
202 : : integer , intent(in ) :: space_cprj
203 : : integer , intent(in ) :: spacecom
204 : : integer , intent(in ) :: spacedim
205 : : integer , intent(in ) :: cprjdim
206 : : integer , intent(in ) :: nslice
207 : : integer , intent(in ) :: paral_slice
208 : : integer , intent(in ) :: spectral_cut
209 : : logical , intent(in ) :: paw
210 : : real(dp) , intent(in ) :: ecut
211 : : real(dp) , intent(in ) :: tolerance
212 : : real(dp) , intent(in ) :: tolfilter
213 : : type(xg_nonlop_t), intent(in ) :: xg_nonlop
214 : : type(slice_t) , intent(inout) :: slice
215 : :
216 : : !Local variables-------------------------------
217 : : real(dp) :: tsec(2)
218 : :
219 : : ! *********************************************************************
220 : :
221 : 0 : call timab(tim_init,1,tsec)
222 : :
223 : 0 : slice%space = space
224 : 0 : slice%space_cprj = space_cprj
225 : 0 : slice%neigenpairs = neigenpairs
226 : 0 : slice%spacedim = spacedim
227 : 0 : slice%bandpp = bandpp
228 : 0 : slice%spacecom = spacecom
229 : 0 : slice%cprjdim = cprjdim
230 : 0 : slice%tolerance = 1.0e-20
231 [ # # ]: 0 : if (tolerance > 0.0) then
232 : 0 : slice%tolerance = tolerance
233 : : end if
234 : 0 : slice%ecut = ecut
235 : 0 : slice%ndeg_filter = ndeg_filter
236 : 0 : slice%nbdbuf = nbdbuf
237 : 0 : slice%eigenProblem = eigenProblem
238 : 0 : slice%me_g0 = me_g0
239 : 0 : slice%me_g0_fft = me_g0_fft
240 : 0 : slice%paw = paw
241 [ # # # # : 0 : slice%xg_nonlop = xg_nonlop
# # # # #
# # # # #
# # # # #
# ]
242 : :
243 : : ! slice specific
244 : 0 : slice%nslice = nslice
245 : 0 : slice%tolfilter = tolfilter
246 : 0 : slice%paral_slice = paral_slice
247 : 0 : slice%spectral_cut = spectral_cut
248 : :
249 : : ! cprj depending on number of bands specific
250 : : !slice%blockdim_cprj = slice%slicedim*xg_nonlop%nspinor
251 : 0 : slice%blockdim_cprj = slice%neigenpairs*xg_nonlop%nspinor
252 : 0 : slice%all_blockdim_cprj = neigenpairs*xg_nonlop%nspinor
253 : :
254 : 0 : call slice_allocateAll(slice)
255 : :
256 : 0 : call timab(tim_init,2,tsec)
257 : :
258 : 0 : end subroutine slice_init
259 : : !!***
260 : :
261 : : !----------------------------------------------------------------------
262 : :
263 : : !!****f* m_slice_cprj/slice_allocateAll
264 : : !! NAME
265 : : !! slice_allocateAll
266 : : !!
267 : : !! FUNCTION
268 : : !! Allocate all memory spaces in a 'slice' datastructure.
269 : : !!
270 : : !! INPUTS
271 : : !!
272 : : !! OUTPUT
273 : : !!
274 : : !! SIDE EFFECTS
275 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
276 : : !!
277 : : !! SOURCE
278 : :
279 : 0 : subroutine slice_allocateAll(slice)
280 : :
281 : : ! Arguments ------------------------------------
282 : : type(slice_t) , intent(inout) :: slice
283 : :
284 : : ! Local variables-------------------------------
285 : : ! scalars
286 : : integer :: neigenpairs
287 : : integer :: space,space_cprj
288 : : integer :: spacedim
289 : : integer :: slicedim
290 : : integer :: total_spacedim, ierr
291 : : integer :: nspinor
292 : :
293 : : ! *********************************************************************
294 : :
295 : 0 : space = slice%space
296 : 0 : space_cprj = slice%space_cprj
297 : 0 : spacedim = slice%spacedim
298 : : !slicedim = slice%slicedim
299 : 0 : slicedim = slice%neigenpairs
300 : 0 : neigenpairs = slice%neigenpairs
301 : 0 : nspinor = slice%xg_nonlop%nspinor
302 : :
303 : 0 : call slice_free(slice)
304 : :
305 : 0 : total_spacedim = spacedim
306 : 0 : call xmpi_sum(total_spacedim,slice%spacecom,ierr)
307 : 0 : slice%total_spacedim = total_spacedim
308 : :
309 : : ! transposed arrays (for slice only)
310 : 0 : call xg_init(slice%X_NP,space,total_spacedim,2*slicedim,xmpi_comm_self,me_g0=slice%me_g0_fft)
311 : 0 : call xg_setBlock(slice%X_NP,slice%X_next,total_spacedim,slicedim)
312 : 0 : call xg_setBlock(slice%X_NP,slice%X_prev,total_spacedim,slicedim,fcol=slicedim+1)
313 : :
314 : : ! This is temporary slice workspace that is the same for all slice initially
315 : 0 : call xg_init(slice%X_SLICE,space,total_spacedim,2*slicedim,xmpi_comm_self,me_g0=slice%me_g0_fft)
316 : 0 : call xg_setBlock(slice%X_SLICE,slice%X,total_spacedim,slicedim)
317 : 0 : call xg_setBlock(slice%X_SLICE,slice%AX,total_spacedim,slicedim,fcol=slicedim+1)
318 : : ! note: the last two are necessary to set space and me_g0 for slice%X and slice%AX
319 : :
320 : : ! cprj workspaces (for entire spectrum)
321 : 0 : call xg_init(slice%AllAX,space,spacedim,neigenpairs,slice%spacecom,me_g0=slice%me_g0)
322 : 0 : call xg_init(slice%Allcprj_work,space_cprj,slice%cprjdim,slice%all_blockdim_cprj,slice%spacecom)
323 : :
324 : : ! cprj workspaces (for slice only)
325 : 0 : call xg_init(slice%cprj_work2,space_cprj,slice%cprjdim,slice%blockdim_cprj,slice%spacecom)
326 : :
327 : : ! cprj workspaces common to slices and entire spectrum
328 : 0 : call xg_init(slice%proj_work,space,slice%xg_nonlop%max_npw_k,slice%xg_nonlop%cprjdim,slice%spacecom,me_g0=slice%me_g0)
329 : :
330 : 0 : end subroutine slice_allocateAll
331 : : !!***
332 : :
333 : : !----------------------------------------------------------------------
334 : :
335 : : !!****f* m_slice_cprj/slice_free
336 : : !! NAME
337 : : !! slice_free
338 : : !!
339 : : !! FUNCTION
340 : : !! Destroy a 'slice' datastructure.
341 : : !!
342 : : !! INPUTS
343 : : !!
344 : : !! OUTPUT
345 : : !! arraymem(2)= memory information
346 : : !!
347 : : !! SIDE EFFECTS
348 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
349 : : !!
350 : : !! SOURCE
351 : :
352 : 0 : subroutine slice_free(slice)
353 : :
354 : : !Arguments ------------------------------------
355 : : type(slice_t) , intent(inout) :: slice
356 : :
357 : : ! *********************************************************************
358 : :
359 : 0 : call xg_free(slice%X_NP)
360 : 0 : call xg_free(slice%X_SLICE)
361 : 0 : call xg_free(slice%Allcprj_work)
362 : 0 : call xg_free(slice%cprj_work2)
363 : 0 : call xg_free(slice%proj_work)
364 : 0 : call xg_free(slice%AllAX)
365 : 0 : call xg_free(slice%Allcprj_work)
366 : :
367 : 0 : end subroutine slice_free
368 : : !!***
369 : :
370 : : !----------------------------------------------------------------------
371 : :
372 : : !!****f* m_slice_cprj/slice_memInfo
373 : : !! NAME
374 : : !! slice_memInfo
375 : : !!
376 : : !! FUNCTION
377 : : !! Provides memory information about a 'slice' datastructure.
378 : : !!
379 : : !! INPUTS
380 : : !! bandpp= number of 'bands' handled by a processor
381 : : !! neigenpairs= number of requested eigenvectors/eigenvalues
382 : : !! paral_kgb= flag controlling (k,g,bands) parallelization
383 : : !! space= defines in which space we are (columns, rows, etc.)
384 : : !! spacedim= dimension of MPI communicator
385 : : !! total_spacedim= size of global KGB communicator (typically 'banspinorfft' comm.)
386 : : !!
387 : : !! OUTPUT
388 : : !! arraymem(2)= memory information
389 : : !!
390 : : !! SIDE EFFECTS
391 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
392 : : !!
393 : : !! SOURCE
394 : :
395 [ # # ]: 0 : function slice_memInfo(neigenpairs,spacedim,space,total_spacedim,bandpp) result(arraymem)
396 : :
397 : : !Arguments ------------------------------------
398 : : integer, intent(in ) :: bandpp
399 : : integer, intent(in ) :: neigenpairs
400 : : integer, intent(in ) :: space
401 : : integer, intent(in ) :: spacedim
402 : : integer, intent(in ) :: total_spacedim
403 : :
404 : : !Local variables-------------------------------
405 : : !scalars
406 : : real(dp) :: memX
407 : : real(dp) :: memX_next
408 : : real(dp) :: memX_prev
409 : : real(dp) :: memAX
410 : : real(dp) :: memBX
411 : : !Transposer variables
412 : : real(dp) :: memX_CR
413 : : real(dp) :: memAX_CR
414 : : real(dp) :: memBX_CR
415 : : !slice_rayleighRitz function variables
416 : : real(dp) :: memA_und_X
417 : : real(dp) :: memB_und_X
418 : : real(dp) :: memEigenvalues
419 : : real(dp) :: cplx
420 : : !arrays
421 : : real(dp) :: arraymem(2)
422 : :
423 : : ! *********************************************************************
424 : 0 : cplx = 1
425 [ # # ]: 0 : if ( space == SPACE_C ) cplx = 2 !for now only complex
426 : :
427 : : !Permanent in slice
428 : 0 : memX = cplx * kind(1.d0) * spacedim * neigenpairs
429 : :
430 : : ! if (paral_kgb == 0) then
431 : : ! memX_next = cplx * kind(1.d0) * spacedim * neigenpairs
432 : : ! memX_prev = cplx * kind(1.d0) * spacedim * neigenpairs
433 : : ! else
434 : 0 : memX_next = cplx * kind(1.d0) * total_spacedim * bandpp
435 : 0 : memX_prev = cplx * kind(1.d0) * total_spacedim * bandpp
436 : : ! end if
437 : :
438 : 0 : memAX = cplx * kind(1.d0) * spacedim * neigenpairs
439 : 0 : memBX = cplx * kind(1.d0) * spacedim * neigenpairs
440 : :
441 : : !Transposer colrow array
442 : : ! if (paral_kgb == 1) then
443 : 0 : memX_CR = cplx * kind(1.d0) * total_spacedim * bandpp
444 : 0 : memAX_CR = cplx * kind(1.d0) * total_spacedim * bandpp
445 : 0 : memBX_CR = cplx * kind(1.d0) * total_spacedim * bandpp
446 : : ! else
447 : : ! memX_CR = 0
448 : : ! memAX_CR = 0
449 : : ! memBX_CR = 0
450 : : ! end if
451 : :
452 : : !slice_rayleighRitz function variables
453 : 0 : memA_und_X = cplx * kind(1.d0) * neigenpairs * neigenpairs
454 : 0 : memB_und_X = cplx * kind(1.d0) * neigenpairs * neigenpairs
455 : 0 : memEigenvalues = kind(1.d0) * neigenpairs
456 : :
457 : : arraymem(1) = memX + memX_next + memX_prev + &
458 : 0 : memAX + memBX + memX_CR + memAX_CR + memBX_CR
459 : 0 : arraymem(2) = memA_und_X + memB_und_X + memEigenvalues
460 : :
461 : 0 : end function slice_memInfo
462 : : !!***
463 : :
464 : : !----------------------------------------------------------------------
465 : :
466 : : !!****f* m_slice_cprj/slice_run_cprj
467 : : !! NAME
468 : : !! slice_run_cprj
469 : : !!
470 : : !! FUNCTION
471 : : !! Apply the Spectrum Slicing algorithm on a set of vectors.
472 : : !!
473 : : !! INPUTS
474 : : !! mpi_enreg = information about MPI parallelization
475 : : !! getAX_BX= pointer to the function giving A|X> and B|X>
476 : : !! A is typically the Hamiltonian H, and B the overlap operator S
477 : : !! getBm1X= pointer to the function giving B^-1|X>
478 : : !! B is typically the overlap operator S
479 : : !!
480 : : !! OUTPUT
481 : : !!
482 : : !! SIDE EFFECTS
483 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
484 : : !! eigen= Full eigenvalues (initial values on entry)
485 : : !! residu= residuals, i.e. norm of (A-lambdaB)|X>
486 : : !! X0= Full set of vectors (initial values on entry)
487 : : !!
488 : : !! SOURCE
489 : :
490 : 0 : subroutine slice_run_cprj(slice,X0,cprjX0,getAX,kin,eigen,occ,residu,enl,nspinor)
491 : :
492 : : !Arguments ------------------------------------
493 : : type(slice_t) , intent(inout) :: slice
494 : : integer, intent(in) :: nspinor
495 : : type(xgBlock_t), intent(inout) :: X0
496 : : type(xgBlock_t), intent(inout) :: cprjX0
497 : : type(xgBlock_t), intent(inout) :: eigen
498 : : type(xgBlock_t), intent(in) :: occ
499 : : type(xgBlock_t), intent(inout) :: residu
500 : : type(xgBlock_t), intent(inout) :: enl
501 : : type(xgBlock_t), intent(in ) :: kin
502 : : interface
503 : : subroutine getAX(X,AX)
504 : : use m_xg, only : xgBlock_t
505 : : type(xgBlock_t), intent(inout) :: X
506 : : type(xgBlock_t), intent(inout) :: AX
507 : : end subroutine getAX
508 : : end interface
509 : :
510 : : !Local variables-------------------------------
511 : : !scalars
512 : : integer :: spacedim
513 : : integer :: space_res
514 : : integer :: neigenpairs
515 : : integer :: nslice
516 : : integer :: spectral_cut
517 : : integer :: paral_slice
518 : : integer :: ndeg_filter,ndeg_filter_max
519 : : integer :: deg_i
520 : : integer :: ideg, ierr, j
521 : : integer :: iband, islice
522 : : integer :: slicedim
523 : : !integer :: shift_x,shift_cprj
524 : : !integer :: niter, max_niter_restart
525 : : integer :: count_mask
526 : : integer :: count_rr
527 : : integer :: count_merge
528 : : integer :: count_slice
529 : : integer :: icount
530 : : integer :: blockdim_cprj
531 : : integer :: fcol_in, lcol_in
532 : : integer :: fcol_dummy, lcol_dummy, ncount_out
533 : : integer :: fcol_global
534 : : integer :: ndeg, ndeg_max
535 : : integer :: tim_slice_fi
536 : : integer :: tim_slice_rr
537 : : integer :: tim_slice_pr
538 : : integer :: my_rank
539 : : integer :: trace_degree, trace_rank
540 : : integer :: nstep_spectrum
541 : : integer :: kmax
542 : : integer :: nfilters, nstep_bisect, ishift
543 : : integer :: nvec_approx
544 : : !integer :: ideg_shift
545 : : integer :: nm
546 : : real(dp) :: balance_prev, balance_this
547 : : real(dp) :: conf_tol
548 : : !real(dp) :: tol_step
549 : : !real(dp) :: tol_probe
550 : : real(dp) :: tolerance
551 : : real(dp) :: tolfilter
552 : : real(dp) :: ramp
553 : : real(dp) :: maxeig, maxeig_global
554 : : real(dp) :: mineig, mineig_global
555 : : real(dp) :: lambda_minus, alpha_minus
556 : : real(dp) :: lambda_plus, alpha_plus
557 : : real(dp) :: ls_in, us_in
558 : : real(dp) :: f_u, f_l, f_uw, f_lw
559 : : real(dp) :: overlap_width
560 : : real(dp) :: one_over_r
561 : : real(dp) :: two_over_r
562 : : real(dp) :: center
563 : : real(dp) :: radius
564 : : real(dp) :: ls, us, cdeg, mu, damp
565 : : real(dp) :: min_low_est
566 : : !real(dp) :: amp_ideg
567 : : !real(dp) :: ein_ideg, eout_ideg
568 : : !real(dp) :: trace_est
569 : : real(dp) :: trace_est_slice1!, trace_est_slice2
570 : : real(dp) :: low_bound, upp_bound, min_low_bound
571 : : real(dp) :: min_upp_bound, max_upp_bound
572 : : real(dp) :: lambda_min, res_norm ! lanczos
573 : : real(dp) :: lower_i, upper_i, mid_i, width
574 : : complex(dp) :: meanz
575 : : real(dp) :: norm2_, variance
576 : : real(dp) :: tol12 = 1.0e-12
577 : : type(xg_t) :: Xsum
578 : : type(xg_t) :: DivResults
579 : : type(xg_t) :: norm2_X
580 : : type(xg_t) :: dot_XfX, norm2_fX
581 : : !type(xg_t) :: res_temp
582 : : type(xg_t) :: X0_out, eigen_out
583 : : type(xgBlock_t) :: X0_out_part, eigen_out_part
584 : : type(xgBlock_t) :: X_in, eigen_in
585 : : !type(xgBlock_t) :: DivResults_part
586 : : !type(xgBlock_t) :: X_prev
587 : : !type(xgBlock_t) :: cprjX_prev
588 : : !type(xgBlock_t) :: cprj_work_prev
589 : : !type(xgBlock_t) :: eig_part
590 : : type(xgBlock_t) :: X_col, AX_col
591 : : type(xg_t) :: X_kept, AX_kept
592 : : type(xg_t) :: cprj_work_slice, cprj_work2_slice
593 : : type(xgBlock_t) :: X_kept_col, AX_kept_col
594 : : type(xgBlock_t) :: eigenvalues_slice, residu_slice
595 : : !type(xg_t) :: X_part
596 : : !type(xg_t) :: AX_part
597 : : !type(xg_t) :: cprjX_part
598 : : integer,parameter :: gpu_option=ABI_GPU_DISABLED
599 : : !arrays
600 : : real(dp) :: tsec(2)
601 : 0 : integer, allocatable :: permute_cols(:)
602 : : !integer, allocatable :: sorted_idx(:) ! same as permute_cols but used elsewhere
603 : 0 : integer, allocatable :: probe_idx(:)
604 : 0 : integer, allocatable :: nb_vec_slices(:)
605 : 0 : real(dp), allocatable :: upper_bound_slices(:)
606 : 0 : real(dp), allocatable :: rayleigh_quotients(:)
607 : 0 : real(dp), allocatable :: energy_filters_left(:,:)
608 : 0 : real(dp), allocatable :: energy_filters_right(:,:)
609 : 0 : real(dp), allocatable :: lower_bounds(:)
610 : 0 : real(dp), allocatable :: upper_bounds(:)
611 : 0 : real(dp), allocatable, target :: cja(:)
612 : 0 : real(dp), allocatable, target :: energy_interval(:)
613 : 0 : real(dp), allocatable :: confi_interval_left(:), confi_interval_right(:)
614 : : real(dp), pointer :: probe(:) => null()
615 : : real(dp), pointer :: probe_XfX(:,:) => null()
616 : : real(dp), pointer :: X0_norm2(:) => null()
617 : : !real(dp), pointer :: lambda_apost(:) => null()
618 : : real(dp), pointer :: lambda_apost_slice(:) => null()
619 : : real(dp), pointer :: theta_(:,:) => null()
620 : : real(dp), pointer :: cheby_moments(:,:) => null()
621 : : real(dp), pointer :: resid(:) => null()
622 : : !Pointers similar to old Chebfi
623 : 0 : type(xg_nonlop_t) :: xg_nonlop
624 : : type(xg_t) :: xgX
625 : : complex(dp), pointer :: Xprobe(:,:) => null()
626 : :
627 : : ! *********************************************************************
628 : :
629 : : ! ITEST
630 : 0 : write(std_out,*) 'inside slice_run'
631 : 0 : write(std_out,*) 'nslice=', slice%nslice
632 : 0 : flush(std_out)
633 : : ! ITEST
634 : :
635 : 0 : write(std_out,*) 'occ array not used', rows(occ), cols(occ)
636 : 0 : flush(std_out)
637 : :
638 : : ! Warning; the entire code assumes this for simplicity and debugging purposes
639 [ # # ]: 0 : ABI_CHECK(slice%bandpp == slice%neigenpairs, "slice_cprj not implemented in MPI")
640 : :
641 : 0 : ABI_WARNING("experimental code is not debugged")
642 : :
643 : : ! Read scalar variables
644 : 0 : nslice = slice%nslice
645 : 0 : slicedim = slice%slicedim
646 : 0 : tolfilter = slice%tolfilter
647 : 0 : spectral_cut = slice%spectral_cut
648 : 0 : paral_slice = slice%paral_slice
649 : 0 : spacedim = slice%spacedim
650 : 0 : neigenpairs = slice%neigenpairs
651 : 0 : ndeg_filter = slice%ndeg_filter
652 : 0 : tolerance = slice%tolerance
653 [ # # # # : 0 : xg_nonlop = slice%xg_nonlop
# # # # #
# # # # #
# # # # #
# ]
654 : :
655 : : ! Set space of results (with symmetry or not)..
656 [ # # ]: 0 : if (slice%space==SPACE_C) then
657 : 0 : space_res = SPACE_C
658 [ # # ]: 0 : else if (slice%space==SPACE_CR) then
659 : 0 : space_res = SPACE_R
660 : : else
661 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
662 : : end if
663 : :
664 : : ! Allocations
665 : : ! DivResults stores Rayleigh quotients
666 : 0 : call xg_init(DivResults, space_res, neigenpairs, 1)
667 [ # # # # : 0 : ABI_MALLOC(permute_cols, (neigenpairs))
# # # # #
# # # ]
668 [ # # # # : 0 : ABI_MALLOC(probe_idx, (neigenpairs))
# # # # ]
669 [ # # # # : 0 : ABI_MALLOC(rayleigh_quotients, (neigenpairs))
# # # # #
# # # ]
670 [ # # # # : 0 : ABI_MALLOC(confi_interval_left, (neigenpairs))
# # # # ]
671 [ # # # # : 0 : ABI_MALLOC(confi_interval_right, (neigenpairs))
# # # # ]
672 [ # # # # : 0 : ABI_MALLOC(nb_vec_slices, (nslice))
# # # # #
# # # ]
673 [ # # # # : 0 : ABI_MALLOC(upper_bound_slices, (nslice))
# # # # #
# # # ]
674 [ # # # # : 0 : ABI_MALLOC(lower_bounds, (nslice))
# # # # ]
675 [ # # # # : 0 : ABI_MALLOC(upper_bounds, (nslice))
# # # # ]
676 : :
677 : : ! Memory used to store results of slice merging.
678 : 0 : call xg_init(X0_out, slice%space, spacedim, neigenpairs, slice%spacecom, me_g0=slice%me_g0)
679 : 0 : call xg_init(eigen_out, SPACE_R, 1, neigenpairs)
680 : :
681 : 0 : call xg_init(norm2_fX, SPACE_R, neigenpairs, 1)
682 : 0 : call xg_init(norm2_X,SPACE_R,neigenpairs,1)
683 : 0 : call xg_init(dot_XfX, space_res, neigenpairs, 1)
684 : :
685 : : ! Set initial vectors from input guess
686 : 0 : slice%eigenvalues = eigen
687 : :
688 : : ! workspace named "All" will store input and output solution
689 : 0 : slice%AllX = X0
690 : 0 : slice%AllcprjX = cprjX0
691 : :
692 : : ! Compute cprjX for all X (in colsrows)
693 : 0 : call timab(tim_cprj,1,tsec)
694 : 0 : call xg_nonlop_getcprj(xg_nonlop,slice%AllX,slice%AllcprjX,slice%proj_work%self)
695 : 0 : call timab(tim_cprj,2,tsec)
696 : :
697 : : !Compute A * Psi
698 : 0 : call timab(tim_AX_v,1,tsec)
699 : 0 : call getAX(slice%AllX,slice%AllAX%self)
700 : 0 : call timab(tim_AX_v,2,tsec)
701 : 0 : call timab(tim_AX_k,1,tsec)
702 : 0 : call xgBlock_add_diag(slice%AllX,kin,nspinor,slice%AllAX%self)
703 : 0 : call timab(tim_AX_k,2,tsec)
704 : 0 : call timab(tim_AX_nl,1,tsec)
705 : 0 : call xg_nonlop_getHX(xg_nonlop,slice%AllAX%self,slice%AllcprjX,slice%Allcprj_work%self,slice%proj_work%self)
706 : 0 : call timab(tim_AX_nl,2,tsec)
707 : :
708 : : !
709 : : ! ------------------------------------------------------------
710 : : ! Compute Rayleigh quotient for every band
711 : : ! ------------------------------------------------------------
712 : : !
713 : :
714 : 0 : call timab(tim_RR_q, 1, tsec)
715 : 0 : call slice_rayleighRitzQuotients(slice, maxeig, mineig, DivResults%self)
716 : 0 : call xmpi_max(maxeig,maxeig_global,slice%spacecom,ierr)
717 : 0 : call xmpi_min(mineig,mineig_global,slice%spacecom,ierr)
718 : 0 : call timab(tim_RR_q, 2, tsec)
719 : :
720 : : ! Recover column indices of sorted Rayleigh quotients in increasing order
721 : 0 : call xgBlock_reverseMap(DivResults%self, theta_, rows=1, cols=neigenpairs)
722 [ # # ]: 0 : rayleigh_quotients(1:neigenpairs) = theta_(1,1:neigenpairs)
723 : :
724 : : ! ITEST
725 : 0 : write(std_out,*) 'rayleigh quotients='
726 : 0 : call xgBlock_print(DivResults%self, std_out)
727 : 0 : flush(std_out)
728 : : ! ITEST
729 : :
730 : : ! Compute |X|^2 colwise L2-norm (before any filter)
731 : 0 : call xgBlock_colwiseNorm2(slice%AllX,norm2_X%self,comm_loc=xmpi_comm_null)
732 : 0 : call xgBlock_reverseMap_1d(norm2_X%self,X0_norm2)
733 : :
734 : : ! ITEST
735 : : !write(std_out,*) 'norm2(squared) ||X||='
736 : : !call xgBlock_print(norm2_X%self,std_out)
737 : : !flush(std_out)
738 : : ! ITEST
739 : : ! Results: |X|=1 so we don't have to normalize everything after..
740 : :
741 : : !
742 : : ! ------------------------------------------------------------
743 : : ! Approximate lowest eigenvalue
744 : : ! ------------------------------------------------------------
745 : : !
746 : :
747 : 0 : kmax = 100
748 : 0 : call computeBLanczos(slice, getAX, kin, spacedim, kmax, lambda_min, res_norm, gpu_option)
749 : 0 : min_low_bound = lambda_min - res_norm
750 : :
751 : 0 : write(std_out,*) 'Lanczos lambda_min=', lambda_min
752 : 0 : write(std_out,*) 'Lanczos res_norm =', res_norm
753 : 0 : write(std_out,*) 'Lanczos guarantee =', min_low_bound
754 : 0 : flush(std_out)
755 : :
756 [ # # ]: 0 : if (res_norm > 0.1d0) then
757 : 0 : ABI_WARNING("Lanczos has residual > 0.1 may need greater kmax to guarantee lower bound")
758 : : end if
759 : :
760 : : !
761 : : ! ------------------------------------------------------------
762 : : ! Split spectrum to slices based on spectral gaps
763 : : ! ------------------------------------------------------------
764 : : !
765 : :
766 : 0 : nstep_spectrum = 30
767 [ # # # # : 0 : min_upp_bound = maxval(rayleigh_quotients)
# # # # ]
768 : 0 : my_rank = xmpi_comm_rank(slice%spacecom)
769 : 0 : max_upp_bound = slice%ecut
770 : 0 : min_low_est = min_low_bound - min(5*res_norm,0.1d0) ! experimental
771 : 0 : trace_rank = neigenpairs ! FIXME for the moment changing this produces a bug
772 : : trace_degree = 50
773 : :
774 : 0 : write(std_out,*) 'splitting spectrum to slices within global', min_low_est, max_upp_bound
775 : 0 : write(std_out,*) ' wanted', min_low_bound, min_upp_bound
776 : 0 : flush(std_out)
777 : : ! call splitSpectrumToSlices(slice, nslice, trace_rank, trace_degree, nstep_spectrum, &
778 : : ! min_low_bound, min_upp_bound, min_low_est, max_upp_bound, my_rank, getAX, kin, &
779 : : ! nb_vec_slices, upper_bound_slices, & ! output
780 : : ! gpu_option=gpu_option)
781 : :
782 : : ! Au31 system
783 : : !upper_bound_slices(1) = 0.10d0 ! TODO hardcoded!
784 : : ! Alu system
785 : 0 : upper_bound_slices(1) = 1.37d0 ! hardcoded
786 : :
787 : : ! Slice 1:
788 : 0 : my_rank = xmpi_comm_rank(slice%spacecom)
789 : 0 : low_bound = min_low_bound ! [a,b)
790 : 0 : upp_bound = upper_bound_slices(1)
791 : 0 : trace_degree = 80
792 : :
793 : : !call computeTraceEstimation(slice, trace_rank, trace_degree, low_bound, upp_bound,&
794 : : ! min_low_est, max_upp_bound, trace_est_slice1, getAX, kin, my_rank, gpu_option=gpu_option)
795 : 0 : write(std_out,*) 'trace estimation for slice1, deg=', trace_est_slice1, trace_degree
796 : 0 : flush(std_out)
797 : :
798 : : ! Slice 2:
799 : 0 : my_rank = xmpi_comm_rank(slice%spacecom)
800 : 0 : low_bound = upper_bound_slices(1) ! [a,b)
801 : 0 : upp_bound = maxval(rayleigh_quotients)
802 : 0 : overlap_width = (upp_bound - low_bound)/10.0
803 : 0 : low_bound = low_bound - overlap_width
804 : : min_low_est = min_low_est
805 : 0 : max_upp_bound = slice%ecut
806 : 0 : trace_degree = 50
807 : 0 : trace_rank = neigenpairs ! FIXME for the moment changing this produces a bug
808 : :
809 : : !call computeTraceEstimation(slice, trace_rank, trace_degree, low_bound, upp_bound,&
810 : : ! min_low_est, max_upp_bound, trace_est_slice2, getAX, kin, my_rank, gpu_option=gpu_option)
811 : : !write(std_out,*) 'trace estimation for slice2, deg=', trace_est_slice2, trace_degree
812 : : !flush(std_out)
813 : :
814 : : ! attention slice%X est modifie n'est plus X0..
815 : : !! normalement il faut remettre slice%X a valeurs de AllX --->
816 : :
817 : :
818 : : !! ------------------------------------------------------------
819 : : !!
820 : : !! - Spectral trees -
821 : : !!
822 : : !! Probe nband with low degree just to prune and reduce to n.
823 : : !! Then actually filter in parallel using high optimal degree.
824 : : !! Spirit of two-step multiresolution.
825 : : !!
826 : : !! Spectral multigrid method. Tree traversal method.
827 : : !! Try to detect number of vectors with spectral trees (hierarchy).
828 : : !! Pay high filter degrees only where spectrum mass lies.
829 : : !!
830 : : !! ------------------------------------------------------------
831 : :
832 : : ! Reset pointers to dimensions of nband (assumes sequential slices)
833 : 0 : call xg_setBlock(slice%X_SLICE,slice%X,slice%total_spacedim,neigenpairs)
834 : 0 : call xg_setBlock(slice%X_SLICE,slice%AX,slice%total_spacedim,neigenpairs,fcol=neigenpairs+1)
835 : :
836 : : ! Initial slice workspaces are independent entire arrays
837 : 0 : call xgBlock_copy(slice%AllX, slice%X)
838 : 0 : call xgBlock_copy(slice%AllAX%self, slice%AX)
839 : :
840 : : ! the two following ones will be recomputed so whatever
841 : 0 : slice%cprjX = slice%AllcprjX
842 : 0 : slice%cprj_work = slice%Allcprj_work%self
843 : :
844 : : ! reinitialize pointers to workspaces ...
845 : 0 : call xg_setBlock(slice%X_NP,slice%X_next,slice%total_spacedim,slice%neigenpairs)
846 : 0 : call xg_setBlock(slice%X_NP,slice%X_prev,slice%total_spacedim,slice%neigenpairs,fcol=slice%neigenpairs+1)
847 : :
848 : 0 : ndeg_filter_max = 80 ! low degree
849 : : ! if energy oscillates instead of being monotonous this means we have to increase
850 : : ! degree
851 : :
852 : 0 : write(std_out,*) 'Spectral trees %%%%%%%%%'
853 : 0 : write(std_out,*) 'ndeg_filter_max=', ndeg_filter_max
854 : 0 : write(std_out,*) 'min_low_est=', min_low_est
855 : 0 : write(std_out,*) 'max_upp_bound=', max_upp_bound
856 : 0 : write(std_out,*) ' ********************** '
857 : 0 : flush(std_out)
858 : :
859 : 0 : call xg_init(xgX, slice%space, spacedim, neigenpairs, slice%spacecom, me_g0=slice%me_g0)
860 : 0 : call xgBlock_reverseMap(xgX%self, Xprobe, spacedim, neigenpairs)
861 [ # # ]: 0 : call generateRademacherMatrix(Xprobe, spacedim, neigenpairs, my_rank)
862 : :
863 : 0 : nm = spacedim * neigenpairs
864 [ # # # # ]: 0 : meanz = sum(Xprobe) / dcmplx(nm,0.0d0)
865 [ # # # # ]: 0 : norm2_ = sum(abs(Xprobe)**2) / nm
866 : 0 : variance = norm2_ - abs(meanz)**2
867 : :
868 : 0 : write(std_out,*) "mean = ", meanz
869 : 0 : write(std_out,*) "E|z|^2 = ", norm2_
870 : 0 : write(std_out,*) "variance = ", variance
871 : 0 : call xgBlock_copy(xgX%self, slice%X)
872 : 0 : write(std_out,*) xgBlock_getid(slice%X)
873 : 0 : flush(std_out)
874 : :
875 : : call computeChebyshevMoments(slice, getAX, kin, min_low_est, max_upp_bound, &
876 : 0 : ndeg_filter_max, cheby_moments, gpu_option)
877 : :
878 : 0 : call xg_free(xgX)
879 : :
880 : : ! Ugly loop to set slice intervals
881 [ # # ]: 0 : do islice=1, nslice
882 : 0 : lambda_minus = upper_bound_slices(1)
883 : 0 : alpha_minus = upper_bound_slices(1)
884 [ # # ]: 0 : if (islice==1) then
885 : 0 : alpha_minus = min_low_bound
886 : 0 : alpha_plus = lambda_minus
887 : : else
888 [ # # # # : 0 : alpha_plus = maxval(rayleigh_quotients)
# # # # ]
889 : : end if
890 : 0 : overlap_width = (alpha_plus - alpha_minus)/10.0
891 : 0 : lambda_plus = alpha_plus+overlap_width
892 : 0 : lambda_minus = alpha_minus-overlap_width
893 [ # # ]: 0 : if (islice==1) then
894 : 0 : lambda_minus = alpha_minus ! otherwise it is outside the center
895 : 0 : lambda_plus = alpha_plus
896 : : end if
897 : 0 : lower_bounds(islice) = lambda_minus
898 : 0 : upper_bounds(islice) = lambda_plus
899 : : end do
900 : :
901 : 0 : nfilters = nslice
902 : 0 : nstep_bisect = 8 ! number of bisection steps
903 : :
904 [ # # # # : 0 : ABI_MALLOC(energy_filters_left, (neigenpairs, nstep_bisect))
# # # # #
# # # #
# ]
905 [ # # # # : 0 : ABI_MALLOC(energy_filters_right, (neigenpairs, nstep_bisect))
# # # # #
# ]
906 [ # # # # : 0 : ABI_MALLOC(energy_interval, (neigenpairs))
# # # # ]
907 [ # # ]: 0 : ABI_MALLOC(cja, (ndeg_filter_max+1))
908 : :
909 : 0 : upper_bounds(2) = 3.0d0 ! hardcoded FIXME auto
910 : 0 : center = (max_upp_bound + min_low_est)*0.5
911 : 0 : radius = (max_upp_bound - min_low_est)*0.5
912 : 0 : lower_i = lower_bounds(1)
913 : 0 : upper_i = upper_bounds(2)
914 : 0 : width = (upper_i - lower_i) / (nstep_bisect + 1)
915 : :
916 : 0 : write(std_out,*)
917 : 0 : write(std_out,*) 'lower_i=', lower_i
918 : 0 : write(std_out,*) 'upper_i=', upper_i
919 : :
920 : : ! total mass
921 : : call buildChebyshevJacksonCoeffs((lower_i-center)/radius, (upper_i-center)/radius, &
922 : 0 : ndeg_filter_max, cja)
923 : : call computeFilterEnergy(neigenpairs, ndeg_filter_max, cja, cheby_moments, &
924 : 0 : energy_interval, nvec_approx)
925 : 0 : write(std_out,*) 'nvec estimate in total', lower_i, upper_i
926 : 0 : write(std_out,*) nvec_approx
927 : 0 : flush(std_out)
928 : :
929 : : ! shifted bisection
930 [ # # ]: 0 : do ishift = 1, nstep_bisect
931 : :
932 : 0 : mid_i = lower_i + ishift * width
933 : 0 : deg_i = ndeg_filter_max
934 : :
935 : 0 : write(std_out,*) '========================================'
936 : 0 : write(std_out,*) 'ishift=', ishift
937 : 0 : write(std_out,*) 'mid_i =', mid_i
938 : :
939 : : ! Slice Left [a,b)
940 : 0 : call buildChebyshevJacksonCoeffs((lower_i-center)/radius, (mid_i-center)/radius, deg_i, cja)
941 : : call computeFilterEnergy(neigenpairs, deg_i, cja, cheby_moments, &
942 : 0 : energy_interval, nvec_approx)
943 : :
944 [ # # ]: 0 : energy_filters_left(:,ishift) = energy_interval(:)
945 : :
946 : 0 : write(std_out,*) 'nvec estimate left from X0 probe=', nvec_approx
947 : 0 : flush(std_out)
948 : :
949 : : ! Slice Right [b,c)
950 : 0 : call buildChebyshevJacksonCoeffs((mid_i-center)/radius, (upper_i-center)/radius, deg_i, cja)
951 : : call computeFilterEnergy(neigenpairs, deg_i, cja, cheby_moments, &
952 : 0 : energy_interval, nvec_approx)
953 : :
954 [ # # ]: 0 : energy_filters_right(:,ishift) = energy_interval(:)
955 : :
956 : 0 : write(std_out,*) 'nvec estimate right from X0 probe=', nvec_approx
957 : 0 : flush(std_out)
958 : :
959 : : !! Uniform mass bisection
960 : : !! si le degre est assez eleve alors la masse totale est
961 : : !! - constante
962 : : !! - egale a nband
963 : : !!
964 : : !! Actuellement la masse totale est 88<192. Ca veut dire qu'on pert de la masse.
965 : : !! Pour ca il faut augmenter l'overlap a droite et a gauche a mon avis.
966 : : !! todo nvec_approx sign flip means that we found important spectral
967 : : !! mass. Do not split there. also follow how the sign flip moves with degree
968 : : !! if it tends to go left or right. Take convergence into account by taking the limit.
969 : :
970 : : end do
971 : :
972 : : ! for every band print its nfilter scores
973 : 0 : write(std_out,*)
974 : 0 : write(std_out,*) 'j= ', nstep_bisect, 'scores'
975 [ # # ]: 0 : do j = 1, neigenpairs
976 : 0 : write(std_out,*) 'eigenvalue', j
977 [ # # ]: 0 : do ishift = 1, nstep_bisect
978 : 0 : balance_this = energy_filters_left(j,ishift) - energy_filters_right(j,ishift)
979 [ # # ]: 0 : if (ishift > 1) then
980 [ # # ]: 0 : if (balance_prev * balance_this < 0) then
981 : 0 : write(std_out,*) 'sign flip! significant mass between=', ishift-1, ishift
982 : : end if
983 : : end if
984 : 0 : write(std_out,*) energy_filters_left(j,ishift), energy_filters_right(j,ishift)
985 : 0 : balance_prev = balance_this
986 : : end do
987 : : end do
988 : 0 : write(std_out,*)
989 : 0 : flush(std_out)
990 : :
991 : : ! sign flip indicates presence of importance spectral mass
992 : : ! we should NOT cut in the (bi,bi+1) that contains the largest spectral mass
993 : : !! adding more nstep_bisect refines the interval (bi,bi+1) containing the large
994 : : !! spectral mass
995 : :
996 [ # # ]: 0 : ABI_FREE(energy_filters_left)
997 [ # # ]: 0 : ABI_FREE(energy_filters_right)
998 [ # # ]: 0 : ABI_FREE(cja)
999 : :
1000 : : !! TODO do the tree traversal
1001 : : !! like subdivide [A,B) and compute energy is each step
1002 : :
1003 : : !! ------------------------------------------------------------
1004 : : !!
1005 : : !! - Main slice loop -
1006 : : !!
1007 : : !! How to optimize degree using energies:
1008 : : !! Increase degree until the band becomes "energetically isolated"
1009 : : !! and further degree increase changes nothing important.
1010 : : !!
1011 : : !! TODO compute Chebyshev recursion for all filters at the same
1012 : : !! time with ONE recursion and MULTIPLE Xsum per slice. Then
1013 : : !! proceed to Rayleigh-Ritz if slice degree reached. Local operation.
1014 : : !!
1015 : : !! ------------------------------------------------------------
1016 : :
1017 : 0 : count_merge = 0
1018 : 0 : fcol_global = 1
1019 : :
1020 [ # # ]: 0 : do islice=1, nslice
1021 : :
1022 : 0 : write(std_out,*)
1023 : 0 : write(std_out,*) '====================Slice=================', islice
1024 : 0 : flush(std_out)
1025 : :
1026 : : !! ------------------------------------------------------------
1027 : : !!
1028 : : !! - Initialize slice workspaces from global workspaces
1029 : : !! Global to Slice operation -
1030 : : !!
1031 : : !! ------------------------------------------------------------
1032 : :
1033 : : ! Reset pointers to dimensions of nband (assumes sequential slices)
1034 : 0 : call xg_setBlock(slice%X_SLICE,slice%X,slice%total_spacedim,neigenpairs)
1035 : 0 : call xg_setBlock(slice%X_SLICE,slice%AX,slice%total_spacedim,neigenpairs,fcol=neigenpairs+1)
1036 : :
1037 : : ! Initial slice workspaces are independent entire arrays
1038 : 0 : call xgBlock_copy(slice%AllX, slice%X)
1039 : 0 : call xgBlock_copy(slice%AllAX%self, slice%AX)
1040 : :
1041 : : ! the two following ones will be recomputed so whatever
1042 : 0 : slice%cprjX = slice%AllcprjX
1043 : 0 : slice%cprj_work = slice%Allcprj_work%self
1044 : :
1045 : : ! reinitialize pointers to workspaces ...
1046 : 0 : call xg_setBlock(slice%X_NP,slice%X_next,slice%total_spacedim,slice%neigenpairs)
1047 : 0 : call xg_setBlock(slice%X_NP,slice%X_prev,slice%total_spacedim,slice%neigenpairs,fcol=slice%neigenpairs+1)
1048 : :
1049 : : !! ------------------------------------------------------------
1050 : : !!
1051 : : !! - Scalar polynomial tuning -
1052 : : !!
1053 : : !! Defining the degree, interval bounds, overlap..
1054 : : !! TODO move outside and before islice loop and store parameters
1055 : : !! into arrays
1056 : : !!
1057 : : !! ------------------------------------------------------------
1058 : :
1059 : : ! ongoing, hardcoded depends on previous code
1060 : 0 : lambda_minus = upper_bound_slices(1)
1061 : 0 : alpha_minus = upper_bound_slices(1)
1062 : :
1063 [ # # ]: 0 : if (islice==1) then
1064 : 0 : alpha_minus = min_low_bound
1065 : 0 : alpha_plus = lambda_minus
1066 : : else
1067 [ # # # # : 0 : alpha_plus = maxval(rayleigh_quotients)
# # # # ]
1068 : : end if
1069 : :
1070 : : ! overlapping between slices
1071 : 0 : overlap_width = (alpha_plus - alpha_minus)/10.0
1072 : 0 : lambda_plus = alpha_plus+overlap_width
1073 : 0 : lambda_minus = alpha_minus-overlap_width
1074 [ # # ]: 0 : if (islice==1) then
1075 : 0 : lambda_minus = alpha_minus ! otherwise it is outside the center
1076 : 0 : lambda_plus = alpha_plus
1077 [ # # ]: 0 : else if (islice==nslice) then
1078 : 0 : lambda_plus = alpha_plus ! there is nothing righwise of max anyway..
1079 : : ! lambda_plus too wide slows down convergence
1080 : : end if
1081 : :
1082 : : ! ITEST
1083 : 0 : write(std_out,*) 'global spectrum=', min_low_est, slice%ecut
1084 : 0 : write(std_out,*) 'wanted slice=', alpha_minus, alpha_plus
1085 : 0 : write(std_out,*) 'with overlap=', lambda_minus, lambda_plus
1086 : 0 : flush(std_out)
1087 : : ! ITEST
1088 : :
1089 : 0 : lower_bounds(islice) = lambda_minus
1090 : 0 : upper_bounds(islice) = lambda_plus
1091 : :
1092 : 0 : center = (slice%ecut + min_low_est)*0.5
1093 : 0 : radius = (slice%ecut - min_low_est)*0.5
1094 : 0 : ls = (lambda_minus - center) / radius
1095 : 0 : us = (lambda_plus - center) / radius
1096 : :
1097 : : ! Optimize polynomial degree with given tolerance
1098 : 0 : ramp = slice%tolfilter
1099 : 0 : ndeg_max = 200
1100 : 0 : ls_in = (alpha_minus - center) / radius ! scaled point inside slice
1101 : 0 : us_in = (alpha_plus - center) / radius ! scaled point inside slice
1102 : 0 : ndeg = 8
1103 : 0 : f_lw = 1.d0; f_uw = 1.d0; f_l = 1.d0; f_u = 1.d0
1104 : :
1105 : : ! Amplification is f(l-w)/f(l) and f(u+w)/f(u) (between 0 and 1)
1106 [ # # ]: 0 : if (islice==1) then
1107 [ # # # # ]: 0 : do while ( f_uw/f_u > ramp .and. ndeg < ndeg_max )
1108 : 0 : ndeg = ndeg + 1
1109 : 0 : f_u = bandpassIndicator_sca(us_in,ls,us,ndeg)
1110 : 0 : f_uw = bandpassIndicator_sca(us ,ls,us,ndeg)
1111 : : end do
1112 : : else
1113 [ # # # # : 0 : do while ( (f_lw/f_l > ramp .or. f_uw/f_u > ramp) .and. ndeg < ndeg_max )
# # ]
1114 : 0 : ndeg = ndeg + 1
1115 : 0 : f_l = bandpassIndicator_sca(ls_in,ls,us,ndeg)
1116 : 0 : f_lw = bandpassIndicator_sca(ls ,ls,us,ndeg)
1117 : 0 : f_u = bandpassIndicator_sca(us_in,ls,us,ndeg)
1118 : 0 : f_uw = bandpassIndicator_sca(us ,ls,us,ndeg)
1119 : : end do
1120 : : end if
1121 : :
1122 : 0 : ndeg = 90
1123 : :
1124 : 0 : ndeg_filter = ndeg
1125 : :
1126 : : ! ITEST
1127 : 0 : write(std_out,*) 'left/right amplif factor f(out)/f(in)=', f_lw/f_l, f_uw/f_u
1128 : 0 : write(std_out,*) 'minimal polynomial degree=', ndeg_filter
1129 : 0 : flush(std_out)
1130 : : ! ITEST
1131 : :
1132 : : !! ------------------------------------------------------------
1133 : : !!
1134 : : !! - Polynomial degree loop -
1135 : : !!
1136 : : !! ------------------------------------------------------------
1137 : :
1138 [ # # ]: 0 : if (islice==1) then
1139 : 0 : tim_slice_fi = tim_slice1_fi
1140 : : else
1141 : 0 : tim_slice_fi = tim_slice2_fi
1142 : : end if
1143 : :
1144 : 0 : call timab(tim_slice_fi,1,tsec)
1145 : :
1146 : : ! Initialize Chebyshev expansion of indicator function of order ndeg_filter
1147 : 0 : call xg_init(Xsum,slice%space,slice%total_spacedim,neigenpairs,slice%spacecom,gpu_option=gpu_option)
1148 : 0 : cdeg = Pi/(ndeg_filter+2)
1149 : 0 : mu = 1.d0/Pi*(ACOS(ls)-ACOS(us))
1150 : 0 : damp = 1.d0 ! Jackson damping
1151 : 0 : one_over_r = 1.d0/radius
1152 : 0 : two_over_r = 2.d0/radius
1153 : 0 : call xgBlock_saxpy(Xsum%self, mu*damp, slice%X)
1154 : :
1155 [ # # ]: 0 : do ideg = 0, ndeg_filter - 1
1156 : :
1157 : : ! ITEST
1158 : 0 : write(std_out,*) 'polynomial degree=', ideg
1159 : 0 : flush(std_out)
1160 : : ! ITEST
1161 : :
1162 : 0 : call timab(tim_cprj,1,tsec)
1163 : 0 : call xg_nonlop_getcprj(xg_nonlop,slice%AX,slice%cprjX,slice%proj_work%self)
1164 : 0 : call timab(tim_cprj,2,tsec)
1165 : :
1166 : 0 : call slice_computeNextOrderChebfiPolynom(slice, ideg, center, one_over_r, two_over_r)
1167 : :
1168 : 0 : call timab(tim_swap,1,tsec)
1169 : 0 : call slice_swapInnerBuffers(slice, slice%total_spacedim, neigenpairs)
1170 : 0 : call timab(tim_swap,2,tsec)
1171 : :
1172 : : ! Accumulate X with weight in Xsum for bandpass filters
1173 : 0 : mu = 2/Pi * (SIN((ideg+1)*ACOS(ls)) - SIN((ideg+1)*ACOS(us)))/(ideg+1)
1174 : : damp = ((1 - (ideg+1)/(ndeg_filter+2))*SIN(cdeg)*COS((ideg+1)*cdeg) + &
1175 : 0 : 1/(ndeg_filter+2)*COS(cdeg)*SIN((ideg+1)*cdeg))/SIN(cdeg)
1176 : 0 : call xgBlock_saxpy(Xsum%self, mu*damp, slice%X)
1177 : :
1178 [ # # ]: 0 : if (ideg==ndeg_filter - 1) then
1179 : :
1180 : : ! store final expansion Xsum to X
1181 : 0 : call xgBlock_copy(Xsum%self, slice%X)
1182 : :
1183 : : end if
1184 : :
1185 : : !A * Psi
1186 : 0 : call timab(tim_AX_v,1,tsec)
1187 : 0 : call getAX(slice%X,slice%AX)
1188 : 0 : call timab(tim_AX_v,2,tsec)
1189 : 0 : call timab(tim_AX_k,1,tsec)
1190 : 0 : call xgBlock_add_diag(slice%X,kin,nspinor,slice%AX)
1191 : 0 : call timab(tim_AX_k,2,tsec)
1192 : 0 : call timab(tim_cprj,1,tsec)
1193 : 0 : call xg_nonlop_getcprj(xg_nonlop,slice%X,slice%cprjX,slice%proj_work%self)
1194 : 0 : call timab(tim_cprj,2,tsec)
1195 : 0 : call timab(tim_AX_nl,1,tsec)
1196 : 0 : call xg_nonlop_getHX(xg_nonlop,slice%AX,slice%cprjX,slice%cprj_work,slice%proj_work%self)
1197 : 0 : call timab(tim_AX_nl,2,tsec)
1198 : :
1199 : : end do ! End polynomial degree loop
1200 : :
1201 : 0 : call timab(tim_slice_fi,2,tsec)
1202 : :
1203 : : !! ------------------------------------------------------------
1204 : : !!
1205 : : !! - Vector pruning for subspace basis -
1206 : : !!
1207 : : !! The purpose of this step is to select the nvec vectors that
1208 : : !! have the greatest energy norm for the slice. No tolerance.
1209 : : !! We start from all vectors then we keep only nvec of them.
1210 : : !! For every column vector x with npw rows,
1211 : : !! * f(x) = ||f(M)x|| (option 1)
1212 : : !! * f(x) = x^T f(M) x (option 2)
1213 : : !!
1214 : : !! ------------------------------------------------------------
1215 : :
1216 [ # # ]: 0 : if (islice==1) then
1217 : 0 : tim_slice_pr = tim_slice1_pr
1218 : : else
1219 : 0 : tim_slice_pr = tim_slice2_pr
1220 : : end if
1221 : :
1222 : 0 : call timab(tim_slice_pr,1,tsec)
1223 : :
1224 : : ! Step 1
1225 : : ! =============
1226 : : ! Compute probe
1227 : : ! Note: option 1 gives betten results than option 2 so far
1228 : : ! =============
1229 : :
1230 [ # # ]: 0 : if (slice%spectral_cut == 1) then
1231 : : ! norm2_fX = <fX,fX> colwise L2-norm
1232 : 0 : call xgBlock_colwiseNorm2(slice%X, norm2_fX%self, comm_loc=xmpi_comm_null)
1233 : 0 : call xgBlock_reverseMap_1d(norm2_fX%self, probe)
1234 [ # # ]: 0 : else if (slice%spectral_cut == 2) then
1235 : : ! dot_XfX = <X,fX> colwise L2-dot product
1236 : 0 : call xgBlock_colwiseDotProduct(X0, slice%X, dot_XfX%self, comm_loc=xmpi_comm_null)
1237 : 0 : call xgBlock_reverseMap(dot_XfX%self, probe_XfX, rows=1, cols=neigenpairs)
1238 : 0 : probe => probe_XfX(1,1:neigenpairs)
1239 : : end if
1240 : :
1241 : : !write(std_out,*)
1242 : : !write(std_out,*) 'probe=', probe(:)
1243 : : !flush(std_out)
1244 : :
1245 : : ! Step 2
1246 : : ! =============
1247 : : ! Probe pruning
1248 : : ! =============
1249 : :
1250 [ # # ]: 0 : if (islice==1) then
1251 : : !count_mask = min(ceiling(trace_est_slice1), neigenpairs)
1252 : : !count_mask = ceiling(trace_est_slice1) + 20
1253 : 0 : count_mask = 150
1254 [ # # ]: 0 : else if (islice==2) then
1255 : : !count_mask = ceiling(neigenpairs - trace_est_slice1*0.8)
1256 : 0 : count_mask = 150
1257 : : end if
1258 : :
1259 [ # # ]: 0 : probe = -probe
1260 [ # # # # : 0 : probe_idx(1:neigenpairs) = (/ (iband, iband=1,neigenpairs) /)
# # ]
1261 [ # # ]: 0 : call sort_dp(neigenpairs, probe, probe_idx, tol12)
1262 [ # # ]: 0 : probe = -probe
1263 : :
1264 : 0 : write(std_out,*) 'kept probes', probe(1:count_mask)
1265 : :
1266 : : ! TODO
1267 : : ! deal with extra vectors: if great probes are found outside the kept ones maybe include them
1268 : : ! appending vectors should be within a loop here
1269 : :
1270 : : ! Step 3
1271 : : ! ==============================================
1272 : : ! Store kept vectors in contiguous memory layout
1273 : : ! ==============================================
1274 : :
1275 : : ! Allocate slice subspace memory, this is contiguous !!
1276 : 0 : call xg_init(X_kept,slice%space,slice%total_spacedim,count_mask,xmpi_comm_self,me_g0=slice%me_g0_fft)
1277 : 0 : call xg_init(AX_kept,slice%space,slice%total_spacedim,count_mask,xmpi_comm_self,me_g0=slice%me_g0_fft)
1278 : :
1279 [ # # ]: 0 : do icount=1,count_mask
1280 : :
1281 : 0 : iband = probe_idx(icount)
1282 : 0 : call xgBlock_setBlock(X_kept%self, X_kept_col, slice%total_spacedim, 1, fcol=icount)
1283 : 0 : call xgBlock_setBlock(AX_kept%self, AX_kept_col, slice%total_spacedim, 1, fcol=icount)
1284 : 0 : call xgBlock_setBlock(slice%X, X_col, slice%total_spacedim, 1, fcol=iband)
1285 : 0 : call xgBlock_setBlock(slice%AX, AX_col, slice%total_spacedim, 1, fcol=iband)
1286 : :
1287 : 0 : call timab(tim_copy, 1, tsec)
1288 : 0 : call xgBlock_copy(X_col, X_kept_col)
1289 : 0 : call xgBlock_copy(AX_col, AX_kept_col)
1290 : 0 : call timab(tim_copy, 2, tsec)
1291 : :
1292 : : end do
1293 : :
1294 : : ! reset pointers to temporary (kept)
1295 : 0 : slice%X = X_kept%self
1296 : 0 : slice%AX = AX_kept%self
1297 : :
1298 : : ! content is not important, but dimensions
1299 : 0 : call xgBlock_setBlock(slice%AllcprjX, slice%cprjX, slice%cprjdim, count_mask*nspinor)
1300 : :
1301 : : ! Recompute cprj to be sure
1302 : 0 : call timab(tim_cprj,1,tsec)
1303 : 0 : call xg_nonlop_getcprj(xg_nonlop,slice%X,slice%cprjX,slice%proj_work%self)
1304 : 0 : call timab(tim_cprj,2,tsec)
1305 : :
1306 : 0 : call timab(tim_slice_pr,2,tsec)
1307 : :
1308 : : !! ------------------------------------------------------------
1309 : : !!
1310 : : !! - Apply Rayleigh-Ritz step -
1311 : : !!
1312 : : !! ------------------------------------------------------------
1313 : :
1314 [ # # ]: 0 : if (islice==1) then
1315 : 0 : tim_slice_rr = tim_slice1_rr
1316 : : else
1317 : 0 : tim_slice_rr = tim_slice2_rr
1318 : : end if
1319 : :
1320 : : ! Number of vectors on which we apply rr
1321 : 0 : count_rr = count_mask
1322 : :
1323 : : ! restrict eigenvalue array for size consistency (essentially keep nonzero entries)
1324 : 0 : call xgblock_reshape(slice%eigenvalues, 1, neigenpairs)
1325 : 0 : call xgblock_setblock(slice%eigenvalues, eigenvalues_slice, rows=1, cols=count_rr)
1326 : 0 : call xgblock_reshape(eigenvalues_slice, count_rr, 1)
1327 : 0 : call xgblock_reshape(slice%eigenvalues, neigenpairs, 1)
1328 : :
1329 : : ! Orthonormalize (ca fait aucune difference)
1330 : : !call xg_Borthonormalize_cprj(xg_nonlop,slice%X,slice%cprjX,ierr,tim_ortho,&
1331 : : ! gpu_option,count_rr*xg_nonlop%nspinor,AX=slice%AX)
1332 : : !=blocksize_cprj
1333 : :
1334 : : ! Apply Rayleigh Ritz on slice (refinement)
1335 : : ! prtvol = 15015015 to print condition number of overlap matrix
1336 : : call xg_RayleighRitz_cprj(xg_nonlop,slice%X,slice%cprjX,slice%AX,eigenvalues_slice,&
1337 : 0 : ierr,15015015,tim_slice_rr,ABI_GPU_DISABLED,solve_ax_bx=.true.)
1338 : :
1339 [ # # ]: 0 : if ( ierr /= 0 ) then
1340 : 0 : ABI_BUG("RayleighRitz did not work")
1341 : : end if
1342 : :
1343 : : ! restart!
1344 : : ! call xg_RayleighRitz_cprj(xg_nonlop,slice%X,slice%cprjX,slice%AX,eigenvalues_slice,&
1345 : : ! ierr,15015015,tim_slice_rr,ABI_GPU_DISABLED,solve_ax_bx=.true.)
1346 : :
1347 : : ! if ( ierr /= 0 ) then
1348 : : ! ABI_BUG("RayleighRitz did not work")
1349 : : ! end if
1350 : :
1351 : : ! ITEST
1352 : 0 : write(std_out,*) 'converged eigenval='
1353 : 0 : call xgBlock_print(eigenvalues_slice, std_out)
1354 : 0 : flush(std_out)
1355 : : ! ITEST
1356 : :
1357 : : ! Restrict dimension of residual array
1358 : 0 : call xgBlock_reshape(residu, 1, neigenpairs)
1359 : 0 : call xgBlock_setBlock(residu, residu_slice, rows=1, cols=count_rr)
1360 : 0 : call xgBlock_reshape(residu_slice, count_rr, 1)
1361 : 0 : call xgBlock_reshape(residu, neigenpairs, 1)
1362 : :
1363 : : !! ------------------------------------------------------------
1364 : : !!
1365 : : !! - Compute residual HX-eSX on slice -
1366 : : !!
1367 : : !! ------------------------------------------------------------
1368 : :
1369 : 0 : blockdim_cprj = count_rr*xg_nonlop%nspinor
1370 : 0 : call xg_init(cprj_work_slice,slice%space_cprj,slice%cprjdim,blockdim_cprj,slice%spacecom)
1371 : 0 : call xg_init(cprj_work2_slice,slice%space_cprj,slice%cprjdim,blockdim_cprj,slice%spacecom)
1372 : :
1373 : : ! Compute H-eSX
1374 [ # # ]: 0 : if (slice%paw) then
1375 : 0 : call timab(tim_AX_nl,1,tsec)
1376 : : call xg_nonlop_getHmeSX(xg_nonlop,slice%X,slice%cprjX,slice%AX,eigenvalues_slice,&
1377 : 0 : cprj_work_slice%self,cprj_work2_slice%self,no_H=.True.)
1378 : 0 : call timab(tim_AX_nl,2,tsec)
1379 : : end if
1380 : :
1381 : 0 : call xg_free(cprj_work_slice)
1382 : 0 : call xg_free(cprj_work2_slice)
1383 : :
1384 : : ! Compute residual norm squared
1385 : 0 : call timab(tim_residu, 1, tsec)
1386 [ # # ]: 0 : if (.not.slice%paw) then
1387 : 0 : call xgBlock_yxmax(slice%AX,eigenvalues_slice,slice%X)
1388 : : end if
1389 : 0 : call xgBlock_colwiseNorm2(slice%AX, residu_slice)
1390 : 0 : call timab(tim_residu, 2, tsec)
1391 : :
1392 : 0 : call xgBlock_reverseMap_1d(residu_slice, resid)
1393 : :
1394 : : ! ITEST
1395 : 0 : write(std_out,*) 'Slice ', islice, ': colwiseNorm2 residu='
1396 : 0 : call xgBlock_print(residu_slice, std_out)
1397 : 0 : flush(std_out)
1398 : : ! ITEST
1399 : :
1400 : : !! ------------------------------------------------------------
1401 : : !!
1402 : : !! - Diagnostic for full slices -
1403 : : !!
1404 : : !! ------------------------------------------------------------
1405 : :
1406 : 0 : call xgBlock_reverseMap_1d(eigenvalues_slice, lambda_apost_slice)
1407 : :
1408 : : ! count the confidence intervals within the slice and outside the slice
1409 : :
1410 [ # # ]: 0 : confi_interval_left(1:count_rr) = lambda_apost_slice + sqrt(resid)
1411 [ # # ]: 0 : confi_interval_right(1:count_rr) = lambda_apost_slice - sqrt(resid)
1412 : :
1413 [ # # ]: 0 : if (islice==1) then
1414 : 0 : write(std_out,*) 'wanted is <', lambda_minus
1415 : : else
1416 : 0 : write(std_out,*) 'wanted is ', alpha_minus, alpha_plus
1417 : : end if
1418 : 0 : flush(std_out)
1419 : :
1420 : : !! TODO etape suivante: une fois qu'on a diagnostiquer une mauvaise convergence
1421 : : !! dans une slice on pourrait faire une procedure de restart pour corriger
1422 : : !! l'erreur soit en ajoutant plus de vecteurs soit jsp a reflechir
1423 : :
1424 : :
1425 : : !! ------------------------------------------------------------
1426 : : !!
1427 : : !! - Slice to Global operation -
1428 : : !!
1429 : : !! Merge to X0_out, eigen_out using Rayleigh value as criterion
1430 : : !!
1431 : : !! ------------------------------------------------------------
1432 : :
1433 : 0 : lcol_in = count_rr
1434 [ # # # # : 0 : fcol_in = maxloc(lambda_apost_slice, dim=1, mask=(lambda_apost_slice < alpha_minus)) + 1
# # # # #
# # # #
# ]
1435 : :
1436 [ # # # # : 0 : conf_tol = maxval(sqrt(resid))
# # # # ]
1437 [ # # # # : 0 : fcol_dummy = maxloc(lambda_apost_slice, dim=1, mask=(lambda_apost_slice+conf_tol < alpha_minus)) + 1
# # # # #
# # # #
# ]
1438 : 0 : write(std_out,*) 'mergeD: fcol slice ', islice, ' ', fcol_in
1439 : 0 : write(std_out,*) 'mergeD: fcol slice ', islice, 'interval', fcol_dummy
1440 [ # # # # ]: 0 : ncount_out = count(lambda_apost_slice+conf_tol < alpha_minus)
1441 : 0 : write(std_out,*) 'mergeD: count out left', islice, ' ', ncount_out, 'out of', count_rr
1442 [ # # # # ]: 0 : ncount_out = count(lambda_apost_slice-conf_tol > alpha_plus)
1443 : 0 : write(std_out,*) 'mergeD: count outright', islice, ' ', ncount_out, 'out of', count_rr
1444 : 0 : write(std_out,*) 'mergeD-----------------------'
1445 : 0 : flush(std_out)
1446 : :
1447 [ # # ]: 0 : if (islice<nslice) then
1448 : :
1449 [ # # # # : 0 : lcol_in = maxloc(lambda_apost_slice, dim=1, mask=(lambda_apost_slice < alpha_plus))
# # # # #
# # # #
# ]
1450 : 0 : write(std_out,*) 'merge: interval on slice ', islice, ':', lcol_in; flush(std_out)
1451 : :
1452 : : else
1453 : :
1454 : 0 : write(std_out,*) 'merge: maximized on slice ', islice, ':', lcol_in
1455 [ # # # # : 0 : lcol_dummy = maxloc(lambda_apost_slice, dim=1, mask=(lambda_apost_slice < alpha_plus))
# # # # #
# # # #
# ]
1456 : 0 : write(std_out,*) 'merge: interval on slice***', islice, ':', lcol_in
1457 [ # # ]: 0 : write(std_out,*) 'merge: resid***************', islice, ':', sqrt(sum(resid(fcol_in:lcol_dummy)))
1458 : 0 : flush(std_out)
1459 : :
1460 [ # # ]: 0 : if (count_merge + lcol_in-fcol_in+1 > neigenpairs) then
1461 : 0 : write(std_out,*) 'merge: pass on slice ', islice, ':', lcol_in; flush(std_out)
1462 : 0 : lcol_in = fcol_in + neigenpairs - count_merge - 1
1463 : 0 : write(std_out,*) 'merge: pass after on slice ', islice, ':', lcol_in, fcol_in, neigenpairs, count_merge
1464 [ # # ]: 0 : write(std_out,*) 'merge: residBBBBBBBBBBBBBBB', islice, ':', sqrt(sum(resid(fcol_in:lcol_in)))
1465 : 0 : flush(std_out)
1466 : : end if
1467 : :
1468 : : end if
1469 : :
1470 : 0 : count_slice = lcol_in - fcol_in + 1
1471 : :
1472 : : ! ITEST
1473 [ # # ]: 0 : write(std_out,*) 'Frobenius norm (inside slice', islice, 'only)=', sqrt(sum(resid(fcol_in:lcol_in)))
1474 : 0 : flush(std_out)
1475 : : ! ITEST
1476 : :
1477 : : ! ITEST
1478 : 0 : write(std_out,*)
1479 : 0 : write(std_out,*) '====================================== Slice', islice
1480 : 0 : write(std_out,*) 'interval bounds', lambda_minus, alpha_minus, alpha_plus
1481 : : !write(std_out,*) 'count_merged=', count_merge
1482 : 0 : write(std_out,*) 'count_slice =', count_slice
1483 : 0 : write(std_out,*) 'count_rr =', count_rr
1484 : 0 : write(std_out,*) 'fcol_in,val =', fcol_in, lambda_apost_slice(fcol_in)
1485 : 0 : write(std_out,*) 'lcol_in,val =', lcol_in, lambda_apost_slice(lcol_in)
1486 : 0 : write(std_out,*) '======================================'
1487 : 0 : write(std_out,*)
1488 : 0 : flush(std_out)
1489 : : ! ITEST
1490 : :
1491 : 0 : count_merge = count_merge + count_slice
1492 [ # # ]: 0 : if (count_merge > neigenpairs) then
1493 : :
1494 : 0 : ABI_WARNING("Attempting to merge more bands than possible")
1495 : :
1496 : : end if
1497 : :
1498 [ # # # # ]: 0 : if (count_merge < neigenpairs .and. islice==nslice) then
1499 : :
1500 : 0 : write(std_out,*) 'missing eigenvalues!'
1501 : 0 : flush(std_out)
1502 : :
1503 : : end if
1504 : :
1505 : 0 : call xgBlock_setBlock(slice%X, X_in, spacedim, count_slice, fcol=fcol_in)
1506 : 0 : call xgBlock_reshape(eigenvalues_slice, 1, count_rr)
1507 : 0 : call xgBlock_setBlock(eigenvalues_slice, eigen_in, 1, count_slice, fcol=fcol_in)
1508 : 0 : call xgBlock_reshape(eigenvalues_slice, count_rr, 1)
1509 : :
1510 : 0 : call xgBlock_setBlock(X0_out%self, X0_out_part, spacedim, count_slice, fcol=fcol_global)
1511 : 0 : call xgBlock_setBlock(eigen_out%self, eigen_out_part, 1, count_slice, fcol=fcol_global)
1512 : :
1513 : 0 : call timab(tim_copy, 1, tsec)
1514 : 0 : call xgBlock_copy(X_in, X0_out_part)
1515 : 0 : call xgBlock_copy(eigen_in, eigen_out_part)
1516 : 0 : call timab(tim_copy, 2, tsec)
1517 : :
1518 : 0 : fcol_global = count_merge + 1
1519 : :
1520 : : !! Free slice memory whose size depends on count_mask, different for every slice
1521 : 0 : call xg_free(Xsum)
1522 : 0 : call xg_free(X_kept)
1523 : 0 : call xg_free(AX_kept)
1524 : :
1525 : : end do ! End loop on slices
1526 : :
1527 : : !! ------------------------------------------------------------
1528 : : !!
1529 : : !! - Final computation of cprjX, residu, enl -
1530 : : !!
1531 : : !! ------------------------------------------------------------
1532 : :
1533 : 0 : call xgBlock_reshape(eigen_out%self, neigenpairs, 1)
1534 : :
1535 : 0 : call timab(tim_copy, 1, tsec)
1536 : 0 : call xgBlock_copy(X0_out%self, X0)
1537 : 0 : call xgBlock_copy(eigen_out%self, eigen)
1538 : 0 : call timab(tim_copy, 2, tsec)
1539 : :
1540 : 0 : slice%eigenvalues = eigen
1541 : 0 : slice%AllX = X0
1542 : 0 : slice%AllcprjX = cprjX0
1543 : :
1544 : 0 : call timab(tim_cprj,1,tsec)
1545 : 0 : call xg_nonlop_getcprj(xg_nonlop,slice%AllX,slice%AllcprjX,slice%proj_work%self)
1546 : 0 : call timab(tim_cprj,2,tsec)
1547 : :
1548 : : ! A * Psi
1549 : 0 : call timab(tim_AX_v,1,tsec)
1550 : 0 : call getAX(slice%AllX,slice%AllAX%self)
1551 : 0 : call timab(tim_AX_v,2,tsec)
1552 : 0 : call timab(tim_AX_k,1,tsec)
1553 : 0 : call xgBlock_add_diag(slice%AllX,kin,nspinor,slice%AllAX%self)
1554 : 0 : call timab(tim_AX_k,2,tsec)
1555 : 0 : call timab(tim_AX_nl,1,tsec)
1556 : 0 : call xg_nonlop_getHX(xg_nonlop,slice%AllAX%self,slice%AllcprjX,slice%Allcprj_work%self,slice%proj_work%self)
1557 : 0 : call timab(tim_AX_nl,2,tsec)
1558 : :
1559 [ # # ]: 0 : if (slice%paw) then
1560 : 0 : call timab(tim_AX_nl,1,tsec)
1561 : : call xg_nonlop_getHmeSX(xg_nonlop,slice%AllX,slice%AllcprjX,slice%AllAX%self,slice%eigenvalues,&
1562 : 0 : slice%Allcprj_work%self,slice%cprj_work2%self,no_H=.True.)
1563 : 0 : call timab(tim_AX_nl,2,tsec)
1564 : : end if
1565 : :
1566 : : ! Compute residual norm squared
1567 : 0 : call timab(tim_residu, 1, tsec)
1568 [ # # ]: 0 : if (.not.slice%paw) then
1569 : 0 : call xgBlock_yxmax(slice%AllAX%self,slice%eigenvalues,slice%AllX)
1570 : : end if
1571 : 0 : call xgBlock_colwiseNorm2(slice%AllAX%self, residu)
1572 : 0 : call timab(tim_residu, 2, tsec)
1573 : :
1574 : : ! ITEST
1575 : 0 : call xgBlock_reverseMap_1d(residu, resid)
1576 [ # # ]: 0 : write(std_out,*) 'Frobenius norm (merged slices)=', sqrt(sum(resid))
1577 : 0 : write(std_out,*) 'resid (merged slices)='
1578 : 0 : call xgBlock_print(residu, std_out)
1579 : 0 : write(std_out,*) 'eigen (merged slices)='
1580 : 0 : call xgBlock_print(eigen, std_out)
1581 : 0 : flush(std_out)
1582 : : ! ITEST
1583 : :
1584 [ # # ]: 0 : if (.not.slice%paw) then
1585 : 0 : call timab(tim_enl,1,tsec)
1586 : 0 : call xg_nonlop_colwiseXHX(xg_nonlop,slice%AllcprjX,slice%Allcprj_work%self,enl)
1587 : 0 : call timab(tim_enl,2,tsec)
1588 : : end if
1589 : :
1590 : : ! Free memory used for slice merging
1591 : 0 : call xg_free(X0_out)
1592 : 0 : call xg_free(eigen_out)
1593 : :
1594 : 0 : call xg_free(norm2_X)
1595 : 0 : call xg_free(norm2_fX)
1596 : 0 : call xg_free(dot_XfX)
1597 : 0 : call xg_free(norm2_fX)
1598 [ # # ]: 0 : ABI_FREE(permute_cols)
1599 [ # # ]: 0 : ABI_FREE(rayleigh_quotients)
1600 [ # # ]: 0 : ABI_FREE(probe_idx)
1601 [ # # ]: 0 : ABI_FREE(confi_interval_left)
1602 [ # # ]: 0 : ABI_FREE(confi_interval_right)
1603 [ # # ]: 0 : ABI_FREE(nb_vec_slices)
1604 : 0 : ABI_FREE(upper_bound_slices)
1605 : 0 : ABI_FREE(lower_bounds)
1606 : 0 : ABI_FREE(upper_bounds)
1607 : :
1608 [ # # # # : 0 : end subroutine slice_run_cprj
# # # # #
# ]
1609 : : !!***
1610 : :
1611 : : !----------------------------------------------------------------------
1612 : :
1613 : : !!****f* m_slice_cprj/slice_orthoXwrtBlocks
1614 : : !! NAME
1615 : : !! slice_orthoXwrtBlocks
1616 : : !!
1617 : : !! FUNCTION
1618 : : !! same as lobpcg_orthoXwrtBlocks but X0 and cprjX0 is given in input
1619 : :
1620 : : subroutine slice_orthoXwrtBlocks(slice,X0,cprjX0,var,cprjvar,islice,cprj_work)
1621 : :
1622 : : type(slice_t) , intent(inout) :: slice
1623 : : type(xgBlock_t), intent(in) :: X0
1624 : : type(xgBlock_t), intent(inout) :: cprjX0
1625 : : type(xgBlock_t), intent(inout) :: var
1626 : : type(xgBlock_t), intent(inout) :: cprjvar
1627 : : type(xgBlock_t), intent(inout) :: cprj_work
1628 : : integer , intent(in ) :: islice
1629 : : integer :: previousBlock
1630 : : integer :: slicedim
1631 : : integer :: spacedim
1632 : : integer :: space_buf
1633 : : integer :: nspinor,cprjdim
1634 : : type(xg_t) :: buffer
1635 : : double precision :: tsec(2)
1636 : : type(xgBlock_t) :: cprjX0_spinor,cprj_work_spinor
1637 : :
1638 : : call timab(tim_ortho,1,tsec)
1639 : :
1640 : : if (islice<2) then
1641 : : ABI_ERROR("islice<2")
1642 : : end if
1643 : :
1644 : : if (cols(cprjvar)/=cols(cprj_work)) then
1645 : : ABI_ERROR("cprjvar and cprj_work should have same number of columns")
1646 : : end if
1647 : : slicedim = slice%slicedim
1648 : : spacedim = slice%spacedim
1649 : : previousBlock = (islice-1)*slice%slicedim
1650 : :
1651 : : cprjdim = slice%xg_nonlop%cprjdim
1652 : : nspinor = slice%xg_nonlop%nspinor
1653 : :
1654 : : space_buf = space(var)
1655 : : if (space(var)==SPACE_CR) then
1656 : : space_buf = SPACE_R
1657 : : end if
1658 : : call xg_init(buffer,space_buf,previousBlock,slicedim,slice%spacecom)
1659 : :
1660 : : ! buffer = X0^T*X
1661 : : call xgBlock_gemm('t','n',1.0d0,X0,var,0.d0,buffer%self,comm=slice%spacecom)
1662 : :
1663 : : ! Add the nonlocal part if paw
1664 : : if (slice%xg_nonlop%paw) then
1665 : : call xg_nonlop_getSX(slice%xg_nonlop,cprjX0,cprjvar,cprj_work,buffer%self)
1666 : : end if
1667 : :
1668 : : ! sum all process contribution
1669 : : ! X = - X0*(BX0^T*X) + X
1670 : : call xgBlock_gemm('n','n',-1.0d0,X0,buffer%self,1.0d0,var)
1671 : :
1672 : : call xgBlock_zero(cprj_work)
1673 : : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,nspinor,COLS2ROWS)
1674 : : call xgBlock_reshape_spinor(cprjX0,cprjX0_spinor,nspinor,COLS2ROWS)
1675 : : call xgBlock_gemm_mpi_cyclic_permutation(cprjX0_spinor,buffer%self,cprj_work_spinor,&
1676 : : & slice%xg_nonlop%me_band,slice%blockdim_cprj/nspinor,comm=slice%xg_nonlop%comm_band)
1677 : : call xgBlock_saxpy(cprjvar,-1.0d0,cprj_work)
1678 : :
1679 : : call xg_free(buffer)
1680 : :
1681 : : call timab(tim_ortho,2,tsec)
1682 : :
1683 : : end subroutine slice_orthoXwrtBlocks
1684 : : !!***
1685 : :
1686 : : !----------------------------------------------------------------------
1687 : :
1688 : : !!****f* m_slice_cprj/slice_rayleighRitzQuotients
1689 : : !! NAME
1690 : : !! slice_rayleighRitzQuotients
1691 : : !!
1692 : : !! FUNCTION
1693 : : !! Compute the Rayleigh-Ritz quotients.
1694 : : !!
1695 : : !! INPUTS
1696 : : !!
1697 : : !! OUTPUT
1698 : : !!
1699 : : !! SIDE EFFECTS
1700 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
1701 : : !! maxeig= highest eigenvalue
1702 : : !! mineig= lowest eigenvalue
1703 : : !! DivResults= Rayleigh-Ritz quotients
1704 : : !!
1705 : : !! SOURCE
1706 : :
1707 : 0 : subroutine slice_rayleighRitzQuotients(slice,maxeig,mineig,DivResults)
1708 : :
1709 : : !Arguments ------------------------------------
1710 : : real(dp), intent(inout) :: maxeig
1711 : : real(dp), intent(inout) :: mineig
1712 : : type(slice_t), intent(inout) :: slice
1713 : : type(xgBlock_t), intent(inout) :: DivResults
1714 : :
1715 : : !Local variables-------------------------------
1716 : : !scalars
1717 : : type(xg_t)::Results1
1718 : : type(xg_t)::Results2
1719 : : type(xg_t)::Results_work
1720 : : !arrays
1721 : : integer :: maxeig_pos(2)
1722 : : integer :: mineig_pos(2)
1723 : : integer :: space_res
1724 : :
1725 : : ! *********************************************************************
1726 : :
1727 : 0 : if (space(slice%AllX)==SPACE_C) then
1728 : 0 : space_res = SPACE_C
1729 [ # # ]: 0 : else if (space(slice%AllX)==SPACE_CR) then
1730 : 0 : space_res = SPACE_R
1731 : : else
1732 : 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
1733 : : end if
1734 : 0 : call xg_init(Results1, space_res, slice%bandpp, 1)
1735 : 0 : call xg_init(Results2, space_res, slice%bandpp, 1)
1736 : :
1737 : 0 : call xgBlock_colwiseDotProduct(slice%AllX,slice%AllAX%self,Results1%self,comm_loc=xmpi_comm_null)
1738 : :
1739 : 0 : call xgBlock_colwiseDotProduct(slice%AllX,slice%AllX,Results2%self,comm_loc=xmpi_comm_null)
1740 [ # # ]: 0 : if (slice%xg_nonlop%paw) then
1741 : 0 : call xg_init(Results_work, space_res, slice%bandpp, 1)
1742 : : call xg_nonlop_colwiseXAX(slice%xg_nonlop,slice%xg_nonlop%Sij%self,slice%AllcprjX,&
1743 : 0 : slice%Allcprj_work%self,Results_work%self)
1744 : 0 : call xgBlock_add(Results2%self,Results_work%self)
1745 : 0 : call xg_free(Results_work)
1746 : : end if
1747 : :
1748 : 0 : call xgBlock_colwiseDivision(Results1%self, Results2%self, DivResults, maxeig, maxeig_pos, mineig, mineig_pos)
1749 : :
1750 : 0 : call xg_free(Results1)
1751 : 0 : call xg_free(Results2)
1752 : :
1753 [ # # ]: 0 : end subroutine slice_rayleighRitzQuotients
1754 : : !!***
1755 : :
1756 : : !----------------------------------------------------------------------
1757 : :
1758 : : !!****f* m_slice_cprj/slice_rayleighRitzQuotientsOnSlice
1759 : : !! NAME
1760 : : !! slice_rayleighRitzQuotientsOnSlice
1761 : : !!
1762 : : !! FUNCTION
1763 : : !! Compute the Rayleigh-Ritz quotients.
1764 : : !!
1765 : : !! INPUTS
1766 : : !!
1767 : : !! OUTPUT
1768 : : !!
1769 : : !! SIDE EFFECTS
1770 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
1771 : : !! maxeig= highest eigenvalue
1772 : : !! mineig= lowest eigenvalue
1773 : : !! DivResults= Rayleigh-Ritz quotients
1774 : : !!
1775 : : !! SOURCE
1776 : :
1777 : : subroutine slice_rayleighRitzQuotientsOnSlice(slice,maxeig,mineig,DivResults)
1778 : :
1779 : : !Arguments ------------------------------------
1780 : : real(dp), intent(inout) :: maxeig
1781 : : real(dp), intent(inout) :: mineig
1782 : : type(slice_t), intent(inout) :: slice
1783 : : type(xgBlock_t), intent(inout) :: DivResults
1784 : :
1785 : : !Local variables-------------------------------
1786 : : !scalars
1787 : : type(xg_t)::Results1
1788 : : type(xg_t)::Results2
1789 : : type(xg_t)::Results_work
1790 : : !arrays
1791 : : integer :: maxeig_pos(2)
1792 : : integer :: mineig_pos(2)
1793 : : integer :: space_res
1794 : :
1795 : : ! *********************************************************************
1796 : :
1797 : : if (space(slice%X)==SPACE_C) then
1798 : : space_res = SPACE_C
1799 : : else if (space(slice%X)==SPACE_CR) then
1800 : : space_res = SPACE_R
1801 : : else
1802 : : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
1803 : : end if
1804 : : call xg_init(Results1, space_res, slice%neigenpairs, 1)
1805 : : call xg_init(Results2, space_res, slice%neigenpairs, 1)
1806 : :
1807 : : call xgBlock_colwiseDotProduct(slice%X,slice%AX,Results1%self,comm_loc=xmpi_comm_null)
1808 : :
1809 : : call xgBlock_colwiseDotProduct(slice%X,slice%X,Results2%self,comm_loc=xmpi_comm_null)
1810 : : if (slice%xg_nonlop%paw) then
1811 : : call xg_init(Results_work, space_res, slice%neigenpairs, 1)
1812 : : call xg_nonlop_colwiseXAX(slice%xg_nonlop,slice%xg_nonlop%Sij%self,slice%cprjX,&
1813 : : slice%cprj_work,Results_work%self)
1814 : : call xgBlock_add(Results2%self,Results_work%self)
1815 : : call xg_free(Results_work)
1816 : : end if
1817 : :
1818 : : call xgBlock_colwiseDivision(Results1%self, Results2%self, DivResults, maxeig, maxeig_pos, mineig, mineig_pos)
1819 : :
1820 : : call xg_free(Results1)
1821 : : call xg_free(Results2)
1822 : :
1823 : : end subroutine slice_rayleighRitzQuotientsOnSlice
1824 : : !!***
1825 : :
1826 : : !----------------------------------------------------------------------
1827 : :
1828 : : !!****f* m_slice_cprj/slice_computeNextOrderChebfiPolynom
1829 : : !! NAME
1830 : : !! slice_computeNextOrderChebfiPolynom
1831 : : !!
1832 : : !! FUNCTION
1833 : : !! From P_n(B-^1.A)|X> (where P_n is the Chebyshev polynom of order n),
1834 : : !! computes P_n+1(B-^1.A)|X>
1835 : : !!
1836 : : !! INPUTS
1837 : : !! ideg=current degree of polynom
1838 : : !! center=filter center
1839 : : !! one_over_r,two_over_r=1/R, 2/R, R being the radius of the filter
1840 : : !! getBm1X= pointer to the function giving B^-1|X>
1841 : : !! B is typically the overlap operator S
1842 : : !!
1843 : : !! OUTPUT
1844 : : !!
1845 : : !! SIDE EFFECTS
1846 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
1847 : : !!
1848 : : !! SOURCE
1849 : :
1850 : 0 : subroutine slice_computeNextOrderChebfiPolynom(slice,ideg,center,one_over_r,two_over_r)
1851 : :
1852 : : !Arguments ------------------------------------
1853 : : real(dp) , intent(in) :: center
1854 : : integer , intent(in) :: ideg
1855 : : real(dp) , intent(in) :: one_over_r
1856 : : real(dp) , intent(in) :: two_over_r
1857 : : type(slice_t) , intent(inout) :: slice
1858 : :
1859 : : !Local variables-------------------------------
1860 : : real(dp) :: tsec(2)
1861 : :
1862 : : ! *********************************************************************
1863 : :
1864 : 0 : call timab(tim_copy, 1, tsec)
1865 : 0 : call xgBlock_copy(slice%AX,slice%X_next)
1866 : 0 : call timab(tim_copy, 2, tsec)
1867 : :
1868 [ # # ]: 0 : if (slice%paw) then
1869 : 0 : call timab(tim_invovl, 1, tsec)
1870 : : call xg_nonlop_getSm1X(slice%xg_nonlop,slice%X_next,slice%cprjX,&
1871 : 0 : & slice%cprj_work,slice%cprj_work2%self,slice%proj_work%self)
1872 : 0 : call timab(tim_invovl, 2, tsec)
1873 : : else
1874 : 0 : call timab(tim_copy, 1, tsec)
1875 : 0 : call xgBlock_copy(slice%AX,slice%X_next)
1876 : 0 : call timab(tim_copy, 2, tsec)
1877 : : end if
1878 : :
1879 : 0 : call timab(tim_postinvovl, 1, tsec)
1880 : 0 : call xgBlock_scale(slice%X, center, 1) !scale by center
1881 : :
1882 : : !(B-1 * A * Psi^i-1 - c * Psi^i-1)
1883 : 0 : call xgBlock_saxpy(slice%X_next, dble(-1.0), slice%X)
1884 : :
1885 : : !Psi^i-1 = 1/c * Psi^i-1
1886 : 0 : call xgBlock_scale(slice%X, 1/center, 1) !counter scale by 1/center
1887 : :
1888 [ # # ]: 0 : if (ideg == 0) then
1889 : 0 : call xgBlock_scale(slice%X_next, one_over_r, 1)
1890 : : else
1891 : 0 : call xgBlock_scale(slice%X_next, two_over_r, 1)
1892 : :
1893 : 0 : call xgBlock_saxpy(slice%X_next, dble(-1.0), slice%X_prev)
1894 : : end if
1895 : :
1896 : 0 : call timab(tim_postinvovl, 2, tsec)
1897 : :
1898 : 0 : end subroutine slice_computeNextOrderChebfiPolynom
1899 : : !!***
1900 : :
1901 : : !----------------------------------------------------------------------
1902 : :
1903 : : !!****f* m_slice_cprj/slice_swapInnerBuffers
1904 : : !! NAME
1905 : : !! slice_swapInnerBuffers
1906 : : !!
1907 : : !! FUNCTION
1908 : : !! Swap buffers inside a 'slice' datastructure.
1909 : : !!
1910 : : !! INPUTS
1911 : : !! ncols= number of requested eigenvectors/eigenvalues
1912 : : !! spacedim= space dimension for one vector
1913 : : !!
1914 : : !! OUTPUT
1915 : : !!
1916 : : !! SIDE EFFECTS
1917 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
1918 : : !!
1919 : : !! SOURCE
1920 : :
1921 : 0 : subroutine slice_swapInnerBuffers(slice,spacedim,ncols)
1922 : :
1923 : : ! Arguments ------------------------------------
1924 : : integer , intent(in ) :: spacedim
1925 : : integer , intent(in ) :: ncols
1926 : : type(slice_t) , intent(inout) :: slice
1927 : :
1928 : : ! *********************************************************************
1929 : :
1930 : 0 : call xgBlock_setBlock(slice%X_prev, slice%X_swap, spacedim, ncols) !X_swap = X_prev
1931 : 0 : call xgBlock_setBlock(slice%X, slice%X_prev, spacedim, ncols) !X_prev = X
1932 : 0 : call xgBlock_setBlock(slice%X_next, slice%X, spacedim, ncols) !X = X_next
1933 : 0 : call xgBlock_setBlock(slice%X_swap, slice%X_next, spacedim, ncols) !X_next = X_swap
1934 : :
1935 : 0 : end subroutine slice_swapInnerBuffers
1936 : : !!***
1937 : :
1938 : : !----------------------------------------------------------------------
1939 : :
1940 : : !!****f* m_slice_cprj/slice_ampfactor
1941 : : !! NAME
1942 : : !! slice_ampfactor
1943 : : !!
1944 : : !! FUNCTION
1945 : : !! Compute amplification factor
1946 : : !!
1947 : : !! INPUTS
1948 : : !! eig (:,:)= eigenvalues
1949 : : !! lambda_minus,lambda_plus=
1950 : : !! ndeg_filter_bands(:)= degree of Spectrum Slicing filter for each band
1951 : : !!
1952 : : !! OUTPUT
1953 : : !!
1954 : : !! SIDE EFFECTS
1955 : : !! residu<type(xgBlock_t)>= vector of residuals
1956 : : !! slice <type(slice_t)>=all data used to apply Spectrum Slicing algorithm
1957 : : !!
1958 : : !! SOURCE
1959 : :
1960 : : subroutine slice_ampfactor(slice,DivResults,lambda_minus,lambda_plus,ndeg_filter_bands)
1961 : :
1962 : : ! Arguments ------------------------------------
1963 : : integer, intent(in ) :: ndeg_filter_bands(:)
1964 : : type(xgBlock_t), intent(in ) :: DivResults
1965 : : real(dp), intent(in ) :: lambda_minus
1966 : : real(dp), intent(in ) :: lambda_plus
1967 : : type(slice_t), intent(inout) :: slice
1968 : :
1969 : : ! Local variables-------------------------------
1970 : : ! scalars
1971 : : integer :: iband
1972 : : !integer :: npw, nband
1973 : : real(dp) :: ampfactor
1974 : : real(dp) :: eig_per_band
1975 : : type(xgBlock_t) :: X_part
1976 : : type(xgBlock_t) :: AX_part
1977 : : real(dp),pointer :: eig(:,:)
1978 : :
1979 : : ! *********************************************************************
1980 : :
1981 : : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=slice%bandpp)
1982 : :
1983 : : do iband = 1, slice%bandpp
1984 : :
1985 : : eig_per_band = eig(1,iband)
1986 : :
1987 : : !cheb_poly1(x, n, a, b)
1988 : : ampfactor = cheb_poly1(eig_per_band, ndeg_filter_bands(iband), lambda_minus, lambda_plus)
1989 : :
1990 : : if(abs(ampfactor) < 1e-3) ampfactor = 1e-3 !just in case, avoid amplifying too much
1991 : :
1992 : : call xgBlock_setBlock(slice%X, X_part, slice%total_spacedim, 1, fcol=iband)
1993 : : call xgBlock_setBlock(slice%AX, AX_part, slice%total_spacedim, 1, fcol=iband)
1994 : :
1995 : : call xgBlock_scale(X_part, 1/ampfactor, 1)
1996 : : call xgBlock_scale(AX_part, 1/ampfactor, 1)
1997 : :
1998 : : end do
1999 : :
2000 : : end subroutine slice_ampfactor
2001 : : !!***
2002 : :
2003 : : subroutine slice_ampfactorMax(slice,DivResults,lambda_minus,lambda_plus,ndeg_filter_bands)
2004 : :
2005 : : ! Arguments ------------------------------------
2006 : : integer, intent(in ) :: ndeg_filter_bands(:)
2007 : : type(xgBlock_t), intent(in ) :: DivResults
2008 : : real(dp), intent(in ) :: lambda_minus
2009 : : real(dp), intent(in ) :: lambda_plus
2010 : : type(slice_t), intent(inout) :: slice
2011 : :
2012 : : ! Local variables-------------------------------
2013 : : ! scalars
2014 : : integer :: iband
2015 : : real(dp) :: ampfactor
2016 : : !type(xgBlock_t) :: X_part
2017 : : !type(xgBlock_t) :: AX_part
2018 : : real(dp),pointer :: eig(:,:)
2019 : :
2020 : : ! *********************************************************************
2021 : :
2022 : : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=cols(DivResults))
2023 : :
2024 : : !cheb_poly1(x, n, a, b)
2025 : : ampfactor = maxval( (/ (cheb_poly1(eig(1,iband), ndeg_filter_bands(iband), lambda_minus, lambda_plus),&
2026 : : iband=1,cols(DivResults)) /) )
2027 : :
2028 : : call xgBlock_scale(slice%X, 1/ampfactor, 1)
2029 : : call xgBlock_scale(slice%AX, 1/ampfactor, 1)
2030 : :
2031 : : end subroutine slice_ampfactorMax
2032 : : !!***
2033 : :
2034 : : !----------------------------------------------------------------------
2035 : :
2036 : : !!****f* m_slice_cprj/slice_probeProximity
2037 : : !! NAME
2038 : : !! slice_probeProximity
2039 : : !!
2040 : : !! FUNCTION
2041 : : !! Compute the principal angles to measure the subspace convergence
2042 : : !! regarding the slice subspace associated to spectral interval [a,b).
2043 : : !! Essentially computes the distance between two subspaces.
2044 : : !! Note that A=Span(a1,..ak) and B=Span(b1,...bk) where families
2045 : : !! are assumed to be linearly independent set of vectors.
2046 : : !! They do not need to be orthogonal!!
2047 : : !!
2048 : : !! IML debug version 08/07:
2049 : : !! probe an interior slice, works for test case only
2050 : :
2051 : : !subroutine slice_probeProximity(slice,islice)
2052 : : !
2053 : : ! implicit none
2054 : : ! type(slice_t), intent(inout) :: slice
2055 : : ! integer, intent(in) :: islice
2056 : : !
2057 : : !
2058 : : !end subroutine slice_probeProximity
2059 : : !!***
2060 : :
2061 : : !----------------------------------------------------------------------
2062 : :
2063 : : !!****f* m_slice_cprj/slice_oracle1
2064 : : !! NAME
2065 : : !! slice_oracle1
2066 : : !!
2067 : : !! FUNCTION
2068 : : !! Compute order of Chebyshev polynom necessary to converge to a given tol
2069 : : !!
2070 : : !! INPUTS
2071 : : !! xx= input variable
2072 : : !! aa= left bound of the interval
2073 : : !! bb= right bound of the interval
2074 : : !! tol= needed precision
2075 : : !! nmax= max number of iterations
2076 : : !!
2077 : : !! OUTPUT
2078 : : !!
2079 : : !! SIDE EFFECTS
2080 : : !!
2081 : : !! SOURCE
2082 : :
2083 : : function cheb_oracle1(xx,aa,bb,tol,nmax) result(nn)
2084 : :
2085 : : ! Arguments ------------------------------------
2086 : : integer :: nn
2087 : : integer, intent(in) :: nmax
2088 : : real(dp), intent(in) :: xx,aa,bb
2089 : : real(dp), intent(in) :: tol
2090 : :
2091 : : ! Local variables-------------------------------
2092 : : integer :: ii
2093 : : real(dp) :: yy,yim1,xred,temp
2094 : :
2095 : : ! *************************************************************************
2096 : :
2097 : : xred = (xx-(aa+bb)/2)/(bb-aa)*2
2098 : : yy = xred
2099 : : yim1 = 1 !ONE
2100 : :
2101 : : nn = nmax
2102 : : if(1/(yy**2) < tol) then
2103 : : nn = 1
2104 : : else
2105 : : do ii=2, nmax-1
2106 : : temp = yy
2107 : : yy = 2*xred*yy - yim1
2108 : : yim1 = temp
2109 : : if(1/(yy**2) < tol) then
2110 : : nn = ii
2111 : : exit
2112 : : end if
2113 : : end do
2114 : : end if
2115 : :
2116 : : end function cheb_oracle1
2117 : : !!***
2118 : :
2119 : : !----------------------------------------------------------------------
2120 : :
2121 : : !!****f* m_slice_cprj/cheb_poly1
2122 : : !! NAME
2123 : : !! cheb_poly1
2124 : : !!
2125 : : !! FUNCTION
2126 : : !! Compute Chebyshev polynomial???
2127 : : !!
2128 : : !! INPUTS
2129 : : !! xx= input variable
2130 : : !! aa= left bound of the interval
2131 : : !! bb= right bound of the interval
2132 : : !! nn=
2133 : : !!
2134 : : !! OUTPUT
2135 : : !!
2136 : : !! SIDE EFFECTS
2137 : : !!
2138 : : !! SOURCE
2139 : :
2140 : : function cheb_poly1(xx,nn,aa,bb) result(yy)
2141 : :
2142 : : ! Arguments ------------------------------------
2143 : : integer, intent(in) :: nn
2144 : : real(dp), intent(in) :: xx, aa, bb
2145 : : real(dp) :: yy
2146 : :
2147 : : ! Local variables-------------------------------
2148 : : integer :: ii
2149 : : real(dp) :: xred,yim1,temp
2150 : :
2151 : : ! *************************************************************************
2152 : :
2153 : : xred = (xx-(aa+bb)/2)/(bb-aa)*2
2154 : : yy = xred
2155 : : yim1 = 1
2156 : : do ii= 2, nn
2157 : : temp = yy
2158 : : yy = 2*xred*yy - yim1
2159 : : yim1 = temp
2160 : : end do
2161 : :
2162 : : end function cheb_poly1
2163 : : !!***
2164 : :
2165 : 0 : function bandpassIndicator_sca(t,a,b,deg) result(f_t)
2166 : :
2167 : : implicit none
2168 : :
2169 : : !Arguments ------------------------------------
2170 : : real(dp), intent(in ) :: t,a,b
2171 : : integer , intent(in ) :: deg
2172 : :
2173 : : real(dp) :: f_t
2174 : :
2175 : : !Local variables-------------------------------
2176 : : real(dp) :: yt0,yt,yt_swap,ck,mu,damp
2177 : : integer :: i
2178 : :
2179 : : ! *********************************************************************
2180 : :
2181 : : ! init cheby of deg=0,1 eval at t
2182 : 0 : yt0 = 1.d0
2183 : 0 : yt = t
2184 : :
2185 : : ! init filter for deg=0
2186 : 0 : ck = Pi/(deg+2)
2187 : 0 : mu = 1/Pi*(ACOS(a)-ACOS(b))
2188 : 0 : damp = 1.d0
2189 : 0 : f_t = mu * damp * yt0
2190 : :
2191 [ # # ]: 0 : do i=1,deg
2192 : :
2193 : : ! Update damping and expansion coefficient
2194 : 0 : mu = 2/Pi * (SIN(i*ACOS(a)) - SIN(i*ACOS(b)))/i
2195 : 0 : damp = ((1 - i/(deg+2))*SIN(ck)*COS(i*ck) + 1/(deg+2)*COS(ck)*SIN(i*ck))/SIN(ck)
2196 : :
2197 : : ! Sum terms
2198 : 0 : f_t = f_t + mu * damp * yt
2199 : :
2200 : : ! Update Chebyshev polynomial
2201 : 0 : yt_swap = yt
2202 : 0 : yt = 2 * t * yt - yt0
2203 : 0 : yt0 = yt_swap
2204 : :
2205 : : end do
2206 : :
2207 : 0 : end function bandpassIndicator_sca
2208 : : !!***
2209 : :
2210 : : !----------------------------------------------------------------------
2211 : :
2212 : : !!****f* m_slice_cprj/generateRademacherMatrix
2213 : : !! NAME
2214 : : !! generateRademacherMatrix
2215 : : !!
2216 : : !! SOURCE
2217 : 0 : subroutine generateRademacherMatrix(V, n, m, rank)
2218 : :
2219 : : implicit none
2220 : :
2221 : : ! input/output
2222 : : integer, intent(in) :: n, m, rank
2223 : : complex(dp), intent(out) :: V(n,m)
2224 : : ! local arguments
2225 : : integer :: nseed, i, j, k
2226 : : integer :: base_seed
2227 : 0 : real(dp) :: R(n,m)
2228 : 0 : integer, allocatable :: seed(:)
2229 : :
2230 : : ! *********************************************************************
2231 : :
2232 : : ! MPI-safe seed: deterministic way to generate a unique seed per MPI rank
2233 : 0 : call random_seed(size=nseed) ! runtime value of nseed
2234 : :
2235 [ # # # # : 0 : ABI_MALLOC(seed, (nseed))
# # # # #
# # # #
# ]
2236 : 0 : base_seed = 123456789
2237 [ # # # # : 0 : seed = mod( base_seed + rank*73856093 + [(i*19349663, i=1,nseed)], 2147483647 )
# # # # #
# # # # #
# # # # ]
2238 : 0 : call random_seed(put=seed)
2239 : :
2240 : 0 : call random_number(R)
2241 [ # # ]: 0 : do j = 1, m
2242 [ # # ]: 0 : do i = 1, n
2243 : 0 : k = int(4.0d0 * R(i,j)) ! 0,1,2,3
2244 : 0 : V(i,j) = dcmplx( cos(0.5d0*pi*k), sin(0.5d0*pi*k))
2245 : : end do
2246 : : end do
2247 : :
2248 [ # # ]: 0 : ABI_FREE(seed)
2249 : :
2250 : 0 : end subroutine generateRademacherMatrix
2251 : : !!***
2252 : :
2253 : : !----------------------------------------------------------------------
2254 : :
2255 : : !!****f* m_slice_cprj/generateGaussianMatrix
2256 : : !! NAME
2257 : : !! generateGaussianMatrix
2258 : : !!
2259 : : !! SOURCE
2260 : : subroutine generateGaussianMatrix(V, n, m, rank)
2261 : :
2262 : : implicit none
2263 : :
2264 : : ! input/output
2265 : : integer, intent(in) :: n, m, rank
2266 : : real(dp), intent(out) :: V(2, n*m)
2267 : :
2268 : : ! local arguments
2269 : : integer :: nseed, i, nm
2270 : : integer :: base_seed
2271 : : integer, allocatable :: seed(:)
2272 : : real(dp) :: u1, u2
2273 : :
2274 : : ! *********************************************************************
2275 : :
2276 : : nm = n * m
2277 : :
2278 : : ! MPI-safe seed: deterministic way to generate a unique seed per MPI rank
2279 : : call random_seed(size = nseed)
2280 : :
2281 : : ABI_MALLOC(seed, (nseed))
2282 : : base_seed = 123456789
2283 : : seed = mod( base_seed + rank*73856093 + [(i*19349663, i=1,nseed)], 2147483647 )
2284 : : call random_seed(put = seed)
2285 : :
2286 : : ! Generate i.i.d. N(0,1) entries (real-valued)
2287 : : i = 1
2288 : : do while (i <= nm)
2289 : : call random_number(u1)
2290 : : call random_number(u2)
2291 : :
2292 : : ! Box-Muller transform
2293 : : V(1, i) = sqrt(-2.0_dp * log(u1)) * cos(2.0_dp * Pi * u2)
2294 : :
2295 : : if (i + 1 <= nm) then
2296 : : V(1, i+1) = sqrt(-2.0_dp * log(u1)) * sin(2.0_dp * Pi * u2)
2297 : : end if
2298 : :
2299 : : i = i + 2
2300 : : end do
2301 : :
2302 : : ! Imaginary part = 0 (consistent with your Rademacher routine)
2303 : : V(2, 1:nm) = 0.0_dp
2304 : :
2305 : : ABI_FREE(seed)
2306 : :
2307 : : end subroutine generateGaussianMatrix
2308 : : !!***
2309 : :
2310 : : !----------------------------------------------------------------------
2311 : :
2312 : : !!****f* m_slice_cprj/applyLowpassFilter
2313 : : !! NAME
2314 : : !! applyLowpassFilter
2315 : : !!
2316 : : !! SOURCE
2317 : : subroutine applyLowpassFilter(slice, getAX, ampl_factors, kin, low_bound, upp_bound, &
2318 : : ndeg_filter, gpu_option)
2319 : :
2320 : : implicit none
2321 : :
2322 : : type(slice_t), intent(inout) :: slice
2323 : : type(xgBlock_t), intent(in) :: kin
2324 : : type(xgBlock_t), intent(in) :: ampl_factors
2325 : : integer, intent(in) :: ndeg_filter
2326 : : real(dp), intent(in) :: low_bound, upp_bound
2327 : : integer, optional, intent(in) :: gpu_option
2328 : : interface
2329 : : subroutine getAX(X,AX)
2330 : : use m_xg, only : xgBlock_t
2331 : : type(xgBlock_t), intent(inout) :: X
2332 : : type(xgBlock_t), intent(inout) :: AX
2333 : : end subroutine getAX
2334 : : end interface
2335 : :
2336 : : integer :: neigenpairs
2337 : : integer :: ideg
2338 : : integer :: nspinor
2339 : : integer :: l_gpu_option
2340 : : real(dp) :: center, radius
2341 : : real(dp) :: one_over_r
2342 : : real(dp) :: two_over_r
2343 : : real(dp) :: tsec(2)
2344 : : integer, allocatable :: ndeg_filter_bands(:) !Oracle variable
2345 : : type(xg_nonlop_t) :: xg_nonlop
2346 : :
2347 : : ! *********************************************************************
2348 : :
2349 : : ! Initialize values
2350 : : neigenpairs = cols(slice%X)
2351 : : xg_nonlop = slice%xg_nonlop
2352 : : nspinor = slice%xg_nonlop%nspinor
2353 : : l_gpu_option = ABI_GPU_DISABLED
2354 : :
2355 : : ! Process input arguments
2356 : : if (present(gpu_option)) then
2357 : : l_gpu_option = gpu_option
2358 : : end if
2359 : :
2360 : : ! Spectral interval to be amplified scaled to [-1,1)
2361 : : center = (upp_bound + low_bound)*0.5
2362 : : radius = (upp_bound - low_bound)*0.5
2363 : : one_over_r = 1.0/radius
2364 : : two_over_r = 2.0/radius
2365 : :
2366 : : ! Loop on degree
2367 : : do ideg = 0, ndeg_filter - 1
2368 : :
2369 : : ! ITEST
2370 : : write(std_out,*) 'polynomial degree=', ideg
2371 : : flush(std_out)
2372 : : ! ITEST
2373 : :
2374 : : call timab(tim_cprj,1,tsec)
2375 : : call xg_nonlop_getcprj(xg_nonlop,slice%AX,slice%cprjX,slice%proj_work%self)
2376 : : call timab(tim_cprj,2,tsec)
2377 : :
2378 : : call slice_computeNextOrderChebfiPolynom(slice, ideg, center, one_over_r, two_over_r)
2379 : :
2380 : : call timab(tim_swap,1,tsec)
2381 : : call slice_swapInnerBuffers(slice, slice%total_spacedim, neigenpairs)
2382 : : call timab(tim_swap,2,tsec)
2383 : :
2384 : : !A * Psi
2385 : : call timab(tim_AX_v,1,tsec)
2386 : : call getAX(slice%X,slice%AX)
2387 : : call timab(tim_AX_v,2,tsec)
2388 : : call timab(tim_AX_k,1,tsec)
2389 : : call xgBlock_add_diag(slice%X,kin,nspinor,slice%AX)
2390 : : call timab(tim_AX_k,2,tsec)
2391 : : call timab(tim_cprj,1,tsec)
2392 : : call xg_nonlop_getcprj(xg_nonlop,slice%X,slice%cprjX,slice%proj_work%self)
2393 : : call timab(tim_cprj,2,tsec)
2394 : : call timab(tim_AX_nl,1,tsec)
2395 : : call xg_nonlop_getHX(xg_nonlop,slice%AX,slice%cprjX,slice%cprj_work,slice%proj_work%self)
2396 : : call timab(tim_AX_nl,2,tsec)
2397 : :
2398 : : end do ! End polynomial degree loop
2399 : :
2400 : : ! Amplify to finalize filter application
2401 : : call timab(tim_amp_f,1,tsec)
2402 : : ABI_MALLOC(ndeg_filter_bands,(neigenpairs))
2403 : : ndeg_filter_bands(:) = ndeg_filter
2404 : : call slice_ampfactor(slice, ampl_factors, low_bound, upp_bound, ndeg_filter_bands)
2405 : : ABI_FREE(ndeg_filter_bands)
2406 : : call timab(tim_amp_f,2,tsec)
2407 : :
2408 : : end subroutine applyLowpassFilter
2409 : : !!***
2410 : :
2411 : : !----------------------------------------------------------------------
2412 : :
2413 : : !!****f* m_slice_cprj/applyBandpassFilter
2414 : : !! NAME
2415 : : !! applyBandpassFilter
2416 : : !!
2417 : : !! SOURCE
2418 : : subroutine applyBandpassFilter(slice, getAX, kin, low_bound, upp_bound, &
2419 : : min_low_bound, max_upp_bound, ndeg_filter, gpu_option)
2420 : :
2421 : : implicit none
2422 : :
2423 : : type(slice_t), intent(inout) :: slice
2424 : : type(xgBlock_t), intent(in) :: kin
2425 : : integer, intent(in) :: ndeg_filter
2426 : : real(dp), intent(in) :: low_bound, upp_bound
2427 : : real(dp), intent(in) :: min_low_bound, max_upp_bound
2428 : : integer, optional, intent(in) :: gpu_option
2429 : : interface
2430 : : subroutine getAX(X,AX)
2431 : : use m_xg, only : xgBlock_t
2432 : : type(xgBlock_t), intent(inout) :: X
2433 : : type(xgBlock_t), intent(inout) :: AX
2434 : : end subroutine getAX
2435 : : end interface
2436 : :
2437 : : integer :: neigenpairs
2438 : : integer :: ideg
2439 : : integer :: nspinor
2440 : : integer :: l_gpu_option
2441 : : real(dp) :: cdeg
2442 : : real(dp) :: center, radius
2443 : : real(dp) :: ls, us
2444 : : real(dp) :: mu, damp
2445 : : real(dp) :: one_over_r
2446 : : real(dp) :: two_over_r
2447 : : real(dp) :: tsec(2)
2448 : : type(xg_t) :: Xsum
2449 : : type(xg_nonlop_t) :: xg_nonlop
2450 : :
2451 : : ! *********************************************************************
2452 : :
2453 : : ! Initialize values
2454 : : neigenpairs = cols(slice%X)
2455 : : xg_nonlop = slice%xg_nonlop
2456 : : nspinor = slice%xg_nonlop%nspinor
2457 : : l_gpu_option = ABI_GPU_DISABLED
2458 : :
2459 : : ! Process input arguments
2460 : : if (present(gpu_option)) then
2461 : : l_gpu_option = gpu_option
2462 : : end if
2463 : :
2464 : : ! Allocate space
2465 : : call xg_init(Xsum,slice%space,slice%total_spacedim,neigenpairs,&
2466 : : slice%spacecom,gpu_option=l_gpu_option)
2467 : :
2468 : : ! Spectral interval to be amplified scaled to [-1,1)
2469 : : center = (max_upp_bound + min_low_bound)*0.5
2470 : : radius = (max_upp_bound - min_low_bound)*0.5
2471 : : ls = (low_bound - center) / radius
2472 : : us = (upp_bound - center) / radius
2473 : : one_over_r = 1.0/radius
2474 : : two_over_r = 2.0/radius
2475 : :
2476 : : ! Initialize Chebyshev expansion of indicator function of order ndeg_filter
2477 : : cdeg = Pi/(ndeg_filter+2)
2478 : : mu = 1.d0/Pi*(ACOS(ls)-ACOS(us))
2479 : : damp = 1.d0 ! Jackson damping
2480 : : call xgBlock_saxpy(Xsum%self, mu*damp, slice%X)
2481 : :
2482 : : ! Loop on degree
2483 : : do ideg = 0, ndeg_filter - 1
2484 : :
2485 : : call timab(tim_cprj,1,tsec)
2486 : : call xg_nonlop_getcprj(xg_nonlop,slice%AX,slice%cprjX,slice%proj_work%self)
2487 : : call timab(tim_cprj,2,tsec)
2488 : :
2489 : : call slice_computeNextOrderChebfiPolynom(slice, ideg, center, one_over_r, two_over_r)
2490 : :
2491 : : call timab(tim_swap,1,tsec)
2492 : : call slice_swapInnerBuffers(slice, slice%total_spacedim, neigenpairs)
2493 : : call timab(tim_swap,2,tsec)
2494 : :
2495 : : ! Accumulate X with weight in Xsum for bandpass filters
2496 : : mu = 2/Pi * (SIN((ideg+1)*ACOS(ls)) - SIN((ideg+1)*ACOS(us)))/(ideg+1)
2497 : : damp = ((1 - (ideg+1)/(ndeg_filter+2))*SIN(cdeg)*COS((ideg+1)*cdeg) + &
2498 : : 1/(ndeg_filter+2)*COS(cdeg)*SIN((ideg+1)*cdeg))/SIN(cdeg)
2499 : : call xgBlock_saxpy(Xsum%self, mu*damp, slice%X)
2500 : :
2501 : : if (ideg==ndeg_filter - 1) then
2502 : :
2503 : : ! store final expansion Xsum to X (inplace rewrites X)
2504 : : call xgBlock_copy(Xsum%self, slice%X)
2505 : :
2506 : : end if
2507 : :
2508 : : !A * Psi
2509 : : call timab(tim_AX_v,1,tsec)
2510 : : call getAX(slice%X,slice%AX)
2511 : : call timab(tim_AX_v,2,tsec)
2512 : : call timab(tim_AX_k,1,tsec)
2513 : : call xgBlock_add_diag(slice%X,kin,nspinor,slice%AX)
2514 : : call timab(tim_AX_k,2,tsec)
2515 : : call timab(tim_cprj,1,tsec)
2516 : : call xg_nonlop_getcprj(xg_nonlop,slice%X,slice%cprjX,slice%proj_work%self)
2517 : : call timab(tim_cprj,2,tsec)
2518 : : call timab(tim_AX_nl,1,tsec)
2519 : : call xg_nonlop_getHX(xg_nonlop,slice%AX,slice%cprjX,slice%cprj_work,slice%proj_work%self)
2520 : : call timab(tim_AX_nl,2,tsec)
2521 : :
2522 : : end do ! End polynomial degree loop
2523 : :
2524 : : ! Free memory
2525 : : call xg_free(Xsum)
2526 : :
2527 : : end subroutine applyBandpassFilter
2528 : : !!***
2529 : :
2530 : : !----------------------------------------------------------------------
2531 : :
2532 : : !!****f* m_slice_cprj/computeChebyshevMoments
2533 : : !! NAME
2534 : : !! computeChebyshevMoments
2535 : : !!
2536 : : !! FUNCTION
2537 : : !! Compute Chebyshev moments up to maximal degree all centered in [A,B)
2538 : : !!
2539 : : !! OUTPUT
2540 : : !! M_n = <X, f_n(B^{-1}AX) X> for n=1,..,ndeg_filter_max
2541 : : !!
2542 : : !! SOURCE
2543 : 0 : subroutine computeChebyshevMoments(slice, getAX, kin, min_low_bound, max_upp_bound, &
2544 : : ndeg_filter_max, cheby_moments, gpu_option)
2545 : :
2546 : : implicit none
2547 : :
2548 : : type(slice_t), intent(inout) :: slice
2549 : : type(xgBlock_t), intent(in) :: kin
2550 : : integer, intent(in) :: ndeg_filter_max
2551 : : real(dp), intent(in) :: min_low_bound, max_upp_bound
2552 : : real(dp), pointer, intent(inout) :: cheby_moments(:,:)
2553 : : integer, optional, intent(in) :: gpu_option
2554 : : interface
2555 : : subroutine getAX(X,AX)
2556 : : use m_xg, only : xgBlock_t
2557 : : type(xgBlock_t), intent(inout) :: X
2558 : : type(xgBlock_t), intent(inout) :: AX
2559 : : end subroutine getAX
2560 : : end interface
2561 : :
2562 : : integer :: neigenpairs
2563 : : integer :: ideg
2564 : : integer :: nspinor
2565 : : integer :: nrows
2566 : : integer :: l_gpu_option
2567 : : real(dp) :: center, radius
2568 : : real(dp) :: one_over_r
2569 : : real(dp) :: two_over_r
2570 : : real(dp) :: tsec(2)
2571 : : type(xg_t) :: Moments, X0
2572 : : type(xgBlock_t) :: moment_ideg
2573 : 0 : type(xg_nonlop_t) :: xg_nonlop
2574 : :
2575 : : ! *********************************************************************
2576 : :
2577 : : ! Initialize values
2578 : 0 : neigenpairs = cols(slice%X)
2579 : 0 : nrows = slice%total_spacedim
2580 [ # # # # : 0 : xg_nonlop = slice%xg_nonlop
# # # # ]
2581 : 0 : nspinor = slice%xg_nonlop%nspinor
2582 : 0 : l_gpu_option = ABI_GPU_DISABLED
2583 : :
2584 : : if (present(gpu_option)) then
2585 : : l_gpu_option = gpu_option
2586 : : end if
2587 : :
2588 : : ! Allocate space
2589 : 0 : call xg_init(Moments, slice%space, neigenpairs, ndeg_filter_max+1) ! M_n=<X0,f_n(A)X0>
2590 : 0 : call xg_init(X0, slice%space, nrows, neigenpairs, slice%spacecom, me_g0=slice%me_g0) ! X0
2591 : :
2592 : 0 : call xgBlock_copy(slice%X, X0%self)
2593 : :
2594 : : ! Spectral interval to be amplified scaled to [-1,1)
2595 : 0 : center = (max_upp_bound + min_low_bound)*0.5
2596 : 0 : radius = (max_upp_bound - min_low_bound)*0.5
2597 : 0 : one_over_r = 1.0/radius
2598 : 0 : two_over_r = 2.0/radius
2599 : :
2600 : : ! Initialize
2601 : 0 : call xgBlock_setBlock(Moments%self, moment_ideg, neigenpairs, 1)
2602 : 0 : call xgBlock_colwiseDotProduct(slice%X, slice%X, moment_ideg, comm_loc=xmpi_comm_null)
2603 : :
2604 : : ! Loop on degree
2605 [ # # ]: 0 : do ideg = 0, ndeg_filter_max - 1
2606 : :
2607 : 0 : call timab(tim_cprj,1,tsec)
2608 : 0 : call xg_nonlop_getcprj(xg_nonlop,slice%AX,slice%cprjX,slice%proj_work%self)
2609 : 0 : call timab(tim_cprj,2,tsec)
2610 : :
2611 : 0 : call slice_computeNextOrderChebfiPolynom(slice, ideg, center, one_over_r, two_over_r)
2612 : :
2613 : : ! slice%X = f_ideg X0
2614 : 0 : call timab(tim_swap,1,tsec)
2615 : 0 : call slice_swapInnerBuffers(slice, slice%total_spacedim, neigenpairs)
2616 : 0 : call timab(tim_swap,2,tsec)
2617 : :
2618 : : ! M_ideg = < X0, f_ideg X0 > in C^nband for every ideg
2619 : 0 : call xgBlock_setBlock(Moments%self, moment_ideg, neigenpairs, 1, fcol=ideg+2)
2620 : 0 : call xgBlock_colwiseDotProduct(X0%self, slice%X, moment_ideg, comm_loc=xmpi_comm_null)
2621 : :
2622 : : !A * Psi
2623 : 0 : call timab(tim_AX_v,1,tsec)
2624 : 0 : call getAX(slice%X,slice%AX)
2625 : 0 : call timab(tim_AX_v,2,tsec)
2626 : 0 : call timab(tim_AX_k,1,tsec)
2627 : 0 : call xgBlock_add_diag(slice%X,kin,nspinor,slice%AX)
2628 : 0 : call timab(tim_AX_k,2,tsec)
2629 : 0 : call timab(tim_cprj,1,tsec)
2630 : 0 : call xg_nonlop_getcprj(xg_nonlop,slice%X,slice%cprjX,slice%proj_work%self)
2631 : 0 : call timab(tim_cprj,2,tsec)
2632 : 0 : call timab(tim_AX_nl,1,tsec)
2633 : 0 : call xg_nonlop_getHX(xg_nonlop,slice%AX,slice%cprjX,slice%cprj_work,slice%proj_work%self)
2634 : 0 : call timab(tim_AX_nl,2,tsec)
2635 : :
2636 : : end do ! End polynomial degree loop
2637 : :
2638 : 0 : call xgBlock_reverseMap(Moments%self, cheby_moments, neigenpairs, ndeg_filter_max+1)
2639 : :
2640 : : ! Free memory
2641 : 0 : call xg_free(X0)
2642 : 0 : call xg_free(Moments)
2643 : :
2644 [ # # # # : 0 : end subroutine computeChebyshevMoments
# # # # #
# ]
2645 : : !!***
2646 : :
2647 : : !----------------------------------------------------------------------
2648 : :
2649 : : !!****f* m_slice_cprj/computeTraceEstimation
2650 : : !! NAME
2651 : : !! computeTraceEstimation
2652 : : !!
2653 : : !! FUNCTION
2654 : : !! Compute Girard-Hutchinson trace estimator
2655 : : !!
2656 : : !! SOURCE
2657 : : subroutine computeTraceEstimation(slice, m_vecs, trace_degree, low_bound, upp_bound, &
2658 : : min_low_bound, max_upp_bound, trace_est, getAX, kin, my_rank, gpu_option)
2659 : :
2660 : : implicit none
2661 : :
2662 : : type(slice_t), intent(inout) :: slice
2663 : : type(xgBlock_t), intent(in) :: kin
2664 : : integer, intent(in) :: my_rank
2665 : : integer, intent(in) :: trace_degree
2666 : : integer, intent(in) :: m_vecs
2667 : : real(dp), intent(out) :: trace_est
2668 : : real(dp), intent(in) :: low_bound, upp_bound
2669 : : real(dp), intent(in) :: min_low_bound, max_upp_bound
2670 : : integer, optional, intent(in) :: gpu_option
2671 : :
2672 : : interface
2673 : : subroutine getAX(X,AX)
2674 : : use m_xg, only : xgBlock_t
2675 : : type(xgBlock_t), intent(inout) :: X
2676 : : type(xgBlock_t), intent(inout) :: AX
2677 : : end subroutine getAX
2678 : : end interface
2679 : :
2680 : : type(xg_t) :: dot_XfX
2681 : : type(xg_t) :: xgX
2682 : : type(xg_nonlop_t) :: xg_nonlop
2683 : : integer :: cprjdim
2684 : : integer :: spacedim
2685 : : integer :: blockdim_cprj
2686 : : integer :: nspinor
2687 : : integer :: l_gpu_option
2688 : : !integer :: idx_i, idx_j
2689 : : real(dp) :: tolerance
2690 : : !real(dp) :: trace_tmp
2691 : : !real(dp) :: normX
2692 : : real(dp), pointer :: accum(:,:) => null()
2693 : : complex(dp), pointer :: X(:,:) => null()
2694 : : real(dp) :: tsec(2)
2695 : :
2696 : : ! *********************************************************************
2697 : :
2698 : : spacedim = slice%total_spacedim
2699 : : tolerance = slice%tolerance
2700 : : cprjdim = slice%cprjdim
2701 : : nspinor = slice%xg_nonlop%nspinor
2702 : : blockdim_cprj = m_vecs*nspinor
2703 : : xg_nonlop = slice%xg_nonlop
2704 : :
2705 : : l_gpu_option = ABI_GPU_DISABLED
2706 : : if (present(gpu_option)) then
2707 : : l_gpu_option = gpu_option
2708 : : end if
2709 : :
2710 : : ! Allocate contiguous memory blocks
2711 : : call xg_init(dot_XfX, slice%space, m_vecs, 1) ! <X,f(A)X>
2712 : : call xg_init(xgX, slice%space, spacedim, m_vecs, slice%spacecom, me_g0=slice%me_g0) ! X
2713 : :
2714 : : call xgBlock_reverseMap(xgX%self, X, spacedim, m_vecs)
2715 : :
2716 : : ! Fill entries of random matrices
2717 : : call generateRademacherMatrix(X, spacedim, m_vecs, my_rank)
2718 : : !call generateGaussianMatrix(X, spacedim, m_vecs, my_rank)
2719 : :
2720 : : ! Initialize X (dimensions then values)
2721 : : call xg_setBlock(slice%X_SLICE, slice%X, spacedim, m_vecs)
2722 : : call xgBlock_copy(xgX%self, slice%X)
2723 : :
2724 : : ! Initialize Xprev, Xnext (dimensions)
2725 : : call xg_setBlock(slice%X_NP,slice%X_next,spacedim,m_vecs)
2726 : : call xg_setBlock(slice%X_NP,slice%X_prev,spacedim,m_vecs,fcol=m_vecs+1)
2727 : :
2728 : : ! Initialize cprj (dimensions and values)
2729 : : call xgBlock_setBlock(slice%AllcprjX, slice%cprjX, slice%cprjdim, blockdim_cprj)
2730 : : call xgBlock_setBlock(slice%Allcprj_work%self, slice%cprj_work, slice%cprjdim, blockdim_cprj)
2731 : :
2732 : : call timab(tim_cprj,1,tsec)
2733 : : call xg_nonlop_getcprj(xg_nonlop,slice%X,slice%cprjX,slice%proj_work%self)
2734 : : call timab(tim_cprj,2,tsec)
2735 : :
2736 : : ! Initialize A * Psi (dimensions and values)
2737 : : call xg_setBlock(slice%X_SLICE, slice%AX, spacedim, m_vecs, fcol=m_vecs+1)
2738 : :
2739 : : call timab(tim_AX_v,1,tsec)
2740 : : call getAX(slice%X,slice%AX)
2741 : : call timab(tim_AX_v,2,tsec)
2742 : : call timab(tim_AX_k,1,tsec)
2743 : : call xgBlock_add_diag(slice%X,kin,nspinor,slice%AX)
2744 : : call timab(tim_AX_k,2,tsec)
2745 : : call timab(tim_AX_nl,1,tsec)
2746 : : call xg_nonlop_getHX(xg_nonlop,slice%AX,slice%cprjX,slice%cprj_work,slice%proj_work%self)
2747 : : call timab(tim_AX_nl,2,tsec)
2748 : :
2749 : : ! Compute f(A) * X
2750 : : call applyBandpassFilter(slice, getAX, kin, low_bound, upp_bound, &
2751 : : min_low_bound, max_upp_bound, trace_degree, l_gpu_option)
2752 : :
2753 : : call xgBlock_colwiseDotProduct(xgX%self, slice%X, dot_XfX%self, comm_loc=xmpi_comm_null)
2754 : : call xgBlock_reverseMap(dot_XfX%self, accum, m_vecs, 1)
2755 : :
2756 : : trace_est = sum(accum) / m_vecs
2757 : :
2758 : : ! Free memory
2759 : : call xg_free(dot_XfX)
2760 : : call xg_free(xgX)
2761 : :
2762 : : end subroutine computeTraceEstimation
2763 : : !!***
2764 : :
2765 : : !----------------------------------------------------------------------
2766 : :
2767 : : !!****f* m_slice_cprj/splitSpectrumToSlices
2768 : : !! NAME
2769 : : !! splitSpectrumToSlices
2770 : : !!
2771 : : !! SOURCE
2772 : : !subroutine splitSpectrumToSlices( &
2773 : : ! slice, nslice, trace_rank, trace_degree, nstep_spectrum, &
2774 : : ! low_bound_wanted, upp_bound_wanted, min_bound, max_bound, my_rank, &
2775 : : ! getAX, kin, nb_vec_slices, upp_bound_slices, gpu_option)
2776 : : subroutine splitSpectrumToSlices( &
2777 : : slice, nslice, nstep_spectrum, &
2778 : : low_bound_wanted, upp_bound_wanted, min_bound, max_bound, &
2779 : : nb_vec_slices, upp_bound_slices, gpu_option)
2780 : :
2781 : : implicit none
2782 : :
2783 : : ! Arguments
2784 : : type(slice_t), intent(inout) :: slice
2785 : : !type(xgBlock_t), intent(in) :: kin
2786 : : integer, intent(in) :: nslice
2787 : : !integer, intent(in) :: trace_rank
2788 : : !integer, intent(in) :: trace_degree
2789 : : integer, intent(in) :: nstep_spectrum
2790 : : !integer, intent(in) :: my_rank
2791 : : real(dp), intent(in) :: low_bound_wanted
2792 : : real(dp), intent(in) :: upp_bound_wanted
2793 : : real(dp), intent(in) :: min_bound
2794 : : real(dp), intent(in) :: max_bound
2795 : : integer, intent(inout) :: nb_vec_slices(nslice)
2796 : : real(dp), intent(inout) :: upp_bound_slices(nslice)
2797 : : integer, optional, intent(in) :: gpu_option
2798 : : !interface
2799 : : ! subroutine getAX(X,AX)
2800 : : ! use m_xg, only : xgBlock_t
2801 : : ! type(xgBlock_t), intent(inout) :: X
2802 : : ! type(xgBlock_t), intent(inout) :: AX
2803 : : ! end subroutine getAX
2804 : : !end interface
2805 : :
2806 : : ! Local variables
2807 : : integer, parameter :: trace_crit = 2
2808 : : integer :: l_gpu_option
2809 : : integer :: trace_sum
2810 : : integer :: neigenpairs
2811 : : integer :: ipart
2812 : : integer :: mid, idx, k, n
2813 : : logical :: found_gap(nstep_spectrum)
2814 : : real(dp) :: width
2815 : : real(dp) :: trace_est
2816 : : real(dp) :: lower_i
2817 : : real(dp) :: upper_i
2818 : : real(dp) :: trace_estim_spectrum(nstep_spectrum)
2819 : : real(dp) :: upper_bound_spectrum(nstep_spectrum)
2820 : :
2821 : : ! *********************************************************************
2822 : :
2823 : : neigenpairs = slice%neigenpairs
2824 : : l_gpu_option = ABI_GPU_DISABLED
2825 : : if (present(gpu_option)) then
2826 : : l_gpu_option = gpu_option
2827 : : end if
2828 : :
2829 : : if (low_bound_wanted < min_bound) then
2830 : : ABI_ERROR('wanted spectrum falls outside given lower bound')
2831 : : end if
2832 : : if (upp_bound_wanted > max_bound) then
2833 : : ABI_ERROR('wanted spectrum falls outside given upper bound')
2834 : : end if
2835 : :
2836 : : ! Split working spectrum into N intervals (N=nstep_spectrum)
2837 : : width = (upp_bound_wanted - low_bound_wanted) / nstep_spectrum
2838 : : trace_sum = 0
2839 : :
2840 : : do ipart=1, nstep_spectrum
2841 : :
2842 : : lower_i = min(max_bound, low_bound_wanted + (ipart - 1) * width)
2843 : : upper_i = min(max_bound, low_bound_wanted + ipart * width)
2844 : :
2845 : : !call computeTraceEstimation(slice, trace_rank, trace_degree, lower_i, upper_i,&
2846 : : ! min_bound, max_bound, trace_est, getAX, kin, my_rank, gpu_option=l_gpu_option)
2847 : :
2848 : : ! TODO keep this version and add second version that reuses Chebyshev
2849 : : ! recursion to avoid repeated calculations...
2850 : :
2851 : : found_gap(ipart) = ( ipart>1 .and. trace_est < trace_crit )
2852 : :
2853 : : write(std_out,*) 'spectrum part : nbvecs', ipart, lower_i, upper_i, ceiling(trace_est)
2854 : : flush(std_out)
2855 : :
2856 : : upper_bound_spectrum(ipart) = upper_i
2857 : : trace_estim_spectrum(ipart) = ceiling(trace_est)
2858 : :
2859 : : trace_sum = trace_sum + ceiling(trace_est)
2860 : :
2861 : : end do
2862 : :
2863 : : if (trace_sum < slice%neigenpairs) then
2864 : : ABI_WARNING("trace estimation missed eigenvalues")
2865 : : end if
2866 : :
2867 : : ! Chose the spectral gap with the index closer to the middle
2868 : : ! same for interval bound closer to the midpoint
2869 : : write(std_out,*) "found_gap=", found_gap
2870 : : write(std_out,*) "upper_bound=", upper_bound_spectrum
2871 : : flush(std_out)
2872 : :
2873 : : found_gap = ( trace_estim_spectrum == minval(trace_estim_spectrum) )
2874 : :
2875 : : write(std_out,*) "found_gap(mod)=", found_gap
2876 : :
2877 : : n = nstep_spectrum
2878 : : mid = (n + 1) / 2
2879 : : idx = 0 ! 0 = not found
2880 : :
2881 : : if (found_gap(mid)) then
2882 : : idx = mid
2883 : : else
2884 : : do k = 1, max(mid-1, n-mid)
2885 : : if (mid-k >= 1) then
2886 : : if (found_gap(mid-k)) then
2887 : : idx = mid-k
2888 : : exit
2889 : : end if
2890 : : end if
2891 : : if (mid+k <= n) then
2892 : : if (found_gap(mid+k)) then
2893 : : idx = mid+k
2894 : : exit
2895 : : end if
2896 : : end if
2897 : : end do
2898 : : end if
2899 : :
2900 : : nb_vec_slices(1) = min(ceiling(sum(trace_estim_spectrum(1:idx+1))), neigenpairs)
2901 : : nb_vec_slices(2) = min(ceiling(sum(trace_estim_spectrum(idx+2:nstep_spectrum))), neigenpairs)
2902 : : write(std_out,*) 'gap idx =', idx
2903 : : write(std_out,*) 'gap value =', upper_bound_spectrum(idx)
2904 : : write(std_out,*) 'trace left =', nb_vec_slices(1)
2905 : : write(std_out,*) 'trace right=', nb_vec_slices(2)
2906 : : flush(std_out)
2907 : :
2908 : : ! output the gap
2909 : : upp_bound_slices(1) = upper_bound_spectrum(idx)
2910 : :
2911 : :
2912 : : end subroutine splitSpectrumToSlices
2913 : : !!***
2914 : :
2915 : : !----------------------------------------------------------------------
2916 : :
2917 : : !!****f* m_slice_cprj/computeBLanczos
2918 : : !! NAME
2919 : : !! computeBLanczos
2920 : : !!
2921 : : !! FUNCTION
2922 : : !! B-Lanczos three-term recurrence (using B-inner product)
2923 : : !!
2924 : : !! SOURCE
2925 : :
2926 : 0 : subroutine computeBLanczos(slice, getAX, kin, n, k, lambda_min, res_norm, gpu_option)
2927 : :
2928 : : implicit none
2929 : :
2930 : : type(slice_t), intent(inout) :: slice
2931 : : type(xgBlock_t), intent(in) :: kin
2932 : : integer, intent(in) :: n, k
2933 : : real(dp), intent(out) :: lambda_min, res_norm
2934 : : integer, optional, intent(in) :: gpu_option
2935 : : interface
2936 : : subroutine getAX(X,AX)
2937 : : use m_xg, only : xgBlock_t
2938 : : type(xgBlock_t), intent(inout) :: X
2939 : : type(xgBlock_t), intent(inout) :: AX
2940 : : end subroutine getAX
2941 : : end interface
2942 : :
2943 : 0 : type(xg_nonlop_t) :: xg_nonlop
2944 : : !type(xg_t) :: xg_Bv
2945 : 0 : real(dp) :: Bv(n), Bv2(2,n)
2946 : 0 : real(dp) :: q(n), q2(2,n), v(n), v2(2,n)
2947 : 0 : real(dp) :: Bm1v2(2,n)
2948 : 0 : real(dp) :: alpha(k), beta(k-1)
2949 : 0 : real(dp) :: q_prev(n)
2950 : 0 : real(dp) :: v_min(n)
2951 : : real(dp) :: beta_prev
2952 : : real(dp) :: normB
2953 : :
2954 : : integer :: i, j
2955 : : integer :: space
2956 : : integer :: spacedim
2957 : : integer :: l_gpu_option
2958 : :
2959 : : ! *********************************************************************
2960 : :
2961 [ # # # # : 0 : xg_nonlop = slice%xg_nonlop
# # # # #
# ]
2962 : 0 : space = slice%space
2963 : 0 : spacedim = slice%spacedim
2964 : :
2965 : 0 : l_gpu_option = ABI_GPU_DISABLED
2966 [ # # ]: 0 : if (present(gpu_option)) then
2967 : 0 : l_gpu_option = gpu_option
2968 : : end if
2969 : :
2970 : : ! ONGOING
2971 : : ! unit tests per elementary operation
2972 : :
2973 : 0 : call random_number(q)
2974 [ # # ]: 0 : do i = 1, n
2975 : 0 : q2(1, i) = q(i)
2976 : 0 : q2(2, i) = 0.0d0
2977 : : end do
2978 [ # # # # : 0 : v2(:,:) = 0.0d0; v(:) = 0.0d0
# # ]
2979 [ # # # # : 0 : Bv2(:,:) = 0.0d0; Bv(:) = 0.0d0
# # ]
2980 [ # # # # ]: 0 : Bm1v2(:,:) = 0.0d0
2981 : :
2982 : : ! Bv = B * q
2983 : 0 : call matmul_op_B(slice, n, q2, Bv2, l_gpu_option)
2984 [ # # ]: 0 : do i = 1, n
2985 : 0 : Bv(i) = Bv2(1, i) ! avoids temporary
2986 : : end do
2987 : :
2988 [ # # ]: 0 : normB = sqrt(dot_product(q, Bv))
2989 [ # # ]: 0 : q = q / normB
2990 : :
2991 [ # # ]: 0 : q_prev = 0.0_dp
2992 : : beta_prev = 0.0_dp
2993 : :
2994 [ # # ]: 0 : do j = 1, k
2995 : : ! v = A * q
2996 [ # # ]: 0 : do i = 1, n
2997 : 0 : q2(1, i) = q(i)
2998 : : end do
2999 : 0 : call matmul_op_A(slice, getAX, kin, n, q2, v2, l_gpu_option)
3000 [ # # ]: 0 : do i = 1, n
3001 : 0 : v(i) = v2(1, i)
3002 : : end do
3003 : :
3004 : : ! alpha_j = q^T * Aq
3005 [ # # ]: 0 : alpha(j) = dot_product(q, v)
3006 : :
3007 : : ! v = B^{-1} * A * q
3008 : 0 : call matmul_op_Binv(slice, n, v2, Bm1v2, l_gpu_option)
3009 [ # # ]: 0 : do i = 1, n
3010 : 0 : v(i) = Bm1v2(1, i)
3011 : : end do
3012 : :
3013 : : ! v = B^{-1} A q - alpha q - beta_prev q_prev
3014 [ # # ]: 0 : v = v - alpha(j)*q
3015 [ # # # # ]: 0 : if (j > 1) v = v - beta_prev*q_prev
3016 : :
3017 : : ! Compute beta_j if j<k
3018 [ # # ]: 0 : if (j < k) then
3019 : : ! Bv = B * v
3020 [ # # ]: 0 : do i = 1, n
3021 : 0 : v2(1, i) = v(i)
3022 : : end do
3023 : 0 : call matmul_op_B(slice, n, v2, Bv2, l_gpu_option)
3024 [ # # ]: 0 : do i = 1, n
3025 : 0 : Bv(i) = Bv2(1, i)
3026 : : end do
3027 : :
3028 [ # # ]: 0 : beta(j) = sqrt(dot_product(v, Bv))
3029 : :
3030 : : ! Update q_prev, q, beta_prev
3031 [ # # ]: 0 : q_prev = q
3032 [ # # ]: 0 : q = v / beta(j)
3033 : : beta_prev = beta(j)
3034 : : end if
3035 : : end do
3036 : :
3037 : : ! Diagonalize T
3038 : 0 : call smallestTridiagEigenpair(k, alpha, beta, lambda_min, v_min)
3039 : :
3040 : : ! residual norm using Lanczos shortcut
3041 : 0 : res_norm = abs(beta(k-1)*v_min(k))
3042 : :
3043 [ # # # # : 0 : end subroutine computeBLanczos
# # # # ]
3044 : : !!***
3045 : :
3046 : : !----------------------------------------------------------------------
3047 : :
3048 : : !!****f* m_slice_cprj/matmul_op_A
3049 : : !! NAME
3050 : : !! matmul_op_A
3051 : : !!
3052 : : !! FUNCTION
3053 : : !! Apply A to Fortran 1D array q2 and store result to Fortran 1D array Aq2
3054 : : !! All in complex.
3055 : : !!
3056 : : !! SOURCE
3057 : :
3058 : 0 : subroutine matmul_op_A(slice, getAX, kin, n, q2, Aq2, gpu_option)
3059 : :
3060 : : implicit none
3061 : :
3062 : : type(slice_t), intent(inout) :: slice
3063 : : type(xgBlock_t), intent(in) :: kin
3064 : : integer, intent(in) :: n
3065 : : real(dp), intent(in) :: q2(2,n) ! 1d complex array (in)
3066 : : real(dp), intent(out) :: Aq2(2,n) ! 1d complex array (out)
3067 : : integer, optional, intent(in) :: gpu_option
3068 : :
3069 : : interface
3070 : : subroutine getAX(X,AX)
3071 : : use m_xg, only : xgBlock_t
3072 : : type(xgBlock_t), intent(inout) :: X
3073 : : type(xgBlock_t), intent(inout) :: AX
3074 : : end subroutine getAX
3075 : : end interface
3076 : :
3077 : : type(xg_t) :: cprjW, cprj_work
3078 : : type(xgBlock_t) :: W, AW, proj_work
3079 : 0 : type(xg_nonlop_t) :: xg_nonlop
3080 : : integer :: l_gpu_option, space, spacecom, space_cprj, cprjdim, nspinor
3081 : : integer :: blockdim_cprj
3082 : : real(dp) :: tsec(2)
3083 : :
3084 : : ! *********************************************************************
3085 : :
3086 : 0 : space = slice%space
3087 : 0 : spacecom = slice%spacecom
3088 : 0 : space_cprj = slice%space_cprj
3089 : 0 : cprjdim = slice%cprjdim
3090 : 0 : nspinor = slice%xg_nonlop%nspinor
3091 [ # # # # : 0 : xg_nonlop = slice%xg_nonlop
# # # # ]
3092 : 0 : proj_work = slice%proj_work%self
3093 : :
3094 : 0 : l_gpu_option = ABI_GPU_DISABLED
3095 : : if (present(gpu_option)) then
3096 : : l_gpu_option = gpu_option
3097 : : end if
3098 : :
3099 : : ! TODO make this optional
3100 : 0 : blockdim_cprj = nspinor ! = number_columns * nspinor
3101 : 0 : call xg_init(cprjW, space_cprj, cprjdim, blockdim_cprj, spacecom)
3102 : 0 : call xg_init(cprj_work, space_cprj, cprjdim, blockdim_cprj, spacecom)
3103 : :
3104 : : ! TODO make this optional otherwise W=slice%X and AW=slice%AX
3105 : 0 : call xgBlock_map(W, q2, space, n, 1, spacecom, gpu_option=gpu_option)
3106 : 0 : call xgBlock_map(AW, Aq2, space, n, 1, spacecom, gpu_option=gpu_option)
3107 : :
3108 : 0 : call timab(tim_cprj,1,tsec)
3109 : 0 : call xg_nonlop_getcprj(xg_nonlop, W, cprjW%self, proj_work)
3110 : 0 : call timab(tim_cprj,2,tsec)
3111 : :
3112 : 0 : call timab(tim_ax_v,1,tsec)
3113 : 0 : call getAX(W,AW)
3114 : 0 : call timab(tim_ax_v,2,tsec)
3115 : :
3116 : 0 : call timab(tim_ax_k,1,tsec)
3117 : 0 : call xgBlock_add_diag(W,kin,nspinor,AW)
3118 : 0 : call timab(tim_ax_k,2,tsec)
3119 : :
3120 : 0 : call timab(tim_AX_nl,1,tsec)
3121 : 0 : call xg_nonlop_getHX(xg_nonlop, AW, cprjW%self, cprj_work%self, proj_work)
3122 : 0 : call timab(tim_AX_nl,2,tsec)
3123 : :
3124 : 0 : call xg_free(cprjW)
3125 : 0 : call xg_free(cprj_work)
3126 : :
3127 [ # # # # : 0 : end subroutine matmul_op_A
# # # # #
# ]
3128 : : !!***
3129 : :
3130 : : !----------------------------------------------------------------------
3131 : :
3132 : : !!****f* m_slice_cprj/matmul_op_B
3133 : : !! NAME
3134 : : !! matmul_op_B
3135 : : !!
3136 : : !! FUNCTION
3137 : : !! Apply B to Fortran 1D array q2 and store result to Fortran 1D array Bq2
3138 : : !! All in complex.
3139 : : !!
3140 : : !! SOURCE
3141 : :
3142 : 0 : subroutine matmul_op_B(slice, n, q2, Bq2, gpu_option)
3143 : :
3144 : : implicit none
3145 : :
3146 : : type(slice_t), intent(inout) :: slice
3147 : : integer, intent(in) :: n
3148 : : real(dp), intent(in) :: q2(2,n) ! 1d complex array (in)
3149 : : real(dp), intent(out) :: Bq2(2,n) ! 1d complex array (out)
3150 : : integer, optional, intent(in) :: gpu_option
3151 : :
3152 : : type(xg_t) :: cprjW, cprj_work
3153 : : type(xgBlock_t) :: W, BW, proj_work
3154 : 0 : type(xg_nonlop_t) :: xg_nonlop
3155 : : integer :: l_gpu_option, space, spacecom, space_cprj, cprjdim, nspinor
3156 : : integer :: blockdim_cprj
3157 : : real(dp) :: tsec(2)
3158 : :
3159 : : ! *********************************************************************
3160 : :
3161 : 0 : space = slice%space
3162 : 0 : spacecom = slice%spacecom
3163 : 0 : space_cprj = slice%space_cprj
3164 : 0 : cprjdim = slice%cprjdim
3165 : 0 : nspinor = slice%xg_nonlop%nspinor
3166 [ # # # # : 0 : xg_nonlop = slice%xg_nonlop
# # # # ]
3167 : 0 : proj_work = slice%proj_work%self
3168 : :
3169 : 0 : l_gpu_option = ABI_GPU_DISABLED
3170 : : if (present(gpu_option)) then
3171 : : l_gpu_option = gpu_option
3172 : : end if
3173 : :
3174 : : ! TODO make this optional
3175 : 0 : blockdim_cprj = nspinor ! = number_columns * nspinor
3176 : 0 : call xg_init(cprjW, space_cprj, cprjdim, blockdim_cprj, spacecom)
3177 : 0 : call xg_init(cprj_work, space_cprj, cprjdim, blockdim_cprj, spacecom)
3178 : :
3179 : : ! TODO make this optional otherwise W=slice%X and AW=slice%AX
3180 : 0 : call xgBlock_map(W, q2, space, n, 1, spacecom, gpu_option=gpu_option)
3181 : 0 : call xgBlock_map(BW, Bq2, space, n, 1, spacecom, gpu_option=gpu_option)
3182 : :
3183 : 0 : call timab(tim_cprj,1,tsec)
3184 : 0 : call xg_nonlop_getcprj(xg_nonlop, W, cprjW%self, proj_work)
3185 : 0 : call timab(tim_cprj,2,tsec)
3186 : :
3187 : : ! BW = S|Psi>
3188 : 0 : call timab(tim_copy, 1, tsec)
3189 : 0 : call xgBlock_copy(W,BW)
3190 : 0 : call timab(tim_copy, 2, tsec)
3191 : :
3192 [ # # ]: 0 : if (slice%paw) then
3193 : 0 : call xg_nonlop_getSX(xg_nonlop, BW, cprjW%self, cprj_work%self, proj_work)
3194 : : end if
3195 : :
3196 : 0 : call xg_free(cprjW)
3197 : 0 : call xg_free(cprj_work)
3198 : :
3199 [ # # # # : 0 : end subroutine matmul_op_B
# # # # #
# ]
3200 : : !!***
3201 : :
3202 : : !----------------------------------------------------------------------
3203 : :
3204 : : !!****f* m_slice_cprj/matmul_op_Binv
3205 : : !! NAME
3206 : : !! matmul_op_Binv
3207 : : !!
3208 : : !! FUNCTION
3209 : : !! Apply Binv to Fortran 1D array q2 and store result to Fortran 1D array Binvq2
3210 : : !! All in complex.
3211 : : !!
3212 : : !! SOURCE
3213 : :
3214 : 0 : subroutine matmul_op_Binv(slice, n, q2, Binvq2, gpu_option)
3215 : :
3216 : : implicit none
3217 : :
3218 : : type(slice_t), intent(inout) :: slice
3219 : : integer, intent(in) :: n
3220 : : real(dp), intent(in) :: q2(2,n) ! 1d complex array (in)
3221 : : real(dp), intent(out) :: Binvq2(2,n) ! 1d complex array (out)
3222 : : integer, optional, intent(in) :: gpu_option
3223 : :
3224 : : type(xg_t) :: cprjW, cprj_work, cprj_work2
3225 : : type(xgBlock_t) :: W, BinvW, proj_work
3226 : 0 : type(xg_nonlop_t) :: xg_nonlop
3227 : : integer :: l_gpu_option, space, spacecom, space_cprj, cprjdim, nspinor
3228 : : integer :: blockdim_cprj
3229 : : real(dp) :: tsec(2)
3230 : :
3231 : : ! *********************************************************************
3232 : :
3233 : 0 : space = slice%space
3234 : 0 : spacecom = slice%spacecom
3235 : 0 : space_cprj = slice%space_cprj
3236 : 0 : cprjdim = slice%cprjdim
3237 : 0 : nspinor = slice%xg_nonlop%nspinor
3238 [ # # # # : 0 : xg_nonlop = slice%xg_nonlop
# # # # ]
3239 : 0 : proj_work = slice%proj_work%self
3240 : :
3241 : 0 : l_gpu_option = ABI_GPU_DISABLED
3242 : : if (present(gpu_option)) then
3243 : : l_gpu_option = gpu_option
3244 : : end if
3245 : :
3246 : 0 : blockdim_cprj = nspinor ! = number_columns * nspinor
3247 : 0 : call xg_init(cprjW, space_cprj, cprjdim, blockdim_cprj, spacecom)
3248 : 0 : call xg_init(cprj_work, space_cprj, cprjdim, blockdim_cprj, spacecom)
3249 : 0 : call xg_init(cprj_work2, space_cprj, cprjdim, blockdim_cprj, spacecom)
3250 : :
3251 : 0 : call xgBlock_map(W, q2, space, n, 1, spacecom, gpu_option=gpu_option)
3252 : 0 : call xgBlock_map(BinvW, Binvq2, space, n, 1, spacecom, gpu_option=gpu_option)
3253 : :
3254 : 0 : call timab(tim_cprj,1,tsec)
3255 : 0 : call xg_nonlop_getcprj(xg_nonlop, W, cprjW%self, proj_work)
3256 : 0 : call timab(tim_cprj,2,tsec)
3257 : :
3258 : : ! BinvW = S^{-1}|Psi>
3259 : 0 : call timab(tim_copy, 1, tsec)
3260 : 0 : call xgBlock_copy(W,BinvW)
3261 : 0 : call timab(tim_copy, 2, tsec)
3262 : :
3263 [ # # ]: 0 : if (slice%paw) then
3264 : 0 : call timab(tim_invovl, 1, tsec)
3265 : : call xg_nonlop_getSm1X(xg_nonlop, BinvW, cprjW%self, cprj_work%self, &
3266 : 0 : cprj_work2%self,proj_work)
3267 : 0 : call timab(tim_invovl, 2, tsec)
3268 : : end if
3269 : :
3270 : 0 : call xg_free(cprjW)
3271 : 0 : call xg_free(cprj_work)
3272 : 0 : call xg_free(cprj_work2)
3273 : :
3274 [ # # # # : 0 : end subroutine matmul_op_Binv
# # # # #
# ]
3275 : : !!***
3276 : :
3277 : : !----------------------------------------------------------------------
3278 : :
3279 : : !!****f* m_slice_cprj/computeFilterEnergy
3280 : : !! NAME
3281 : : !! computeFilterEnergy
3282 : : !!
3283 : : !! SOURCE
3284 : :
3285 : 0 : subroutine computeFilterEnergy(neigenpairs, ndeg_filter, cja, cheby_moments, &
3286 : : energy_interval, nvec_approx)
3287 : :
3288 : : implicit none
3289 : :
3290 : : integer, intent(in) :: ndeg_filter
3291 : : integer, intent(in) :: neigenpairs
3292 : : integer, intent(out) :: nvec_approx
3293 : : real(dp), pointer, intent(in) :: cja(:)
3294 : : real(dp), pointer, intent(in) :: cheby_moments(:,:)
3295 : : real(dp), pointer, intent(in) :: energy_interval(:)
3296 : :
3297 : 0 : real(dp) :: E_re(neigenpairs)
3298 : 0 : real(dp) :: E_im(neigenpairs)
3299 : : real(dp) :: E_tot
3300 : : real(dp) :: energy_magn_j
3301 : : integer :: ideg
3302 : : integer :: j
3303 : :
3304 : : ! *********************************************************************
3305 : :
3306 [ # # ]: 0 : E_re = 0.0d0
3307 [ # # ]: 0 : E_im = 0.0d0
3308 : 0 : E_tot = 0.d0
3309 [ # # ]: 0 : do ideg = 1, ndeg_filter+1
3310 [ # # ]: 0 : do j = 1, neigenpairs
3311 : 0 : E_re(j) = E_re(j) + cja(ideg) * cheby_moments(2*j-1, ideg)
3312 : 0 : E_im(j) = E_im(j) + cja(ideg) * cheby_moments(2*j, ideg)
3313 : : end do
3314 : : end do
3315 [ # # ]: 0 : do j = 1, neigenpairs
3316 : : !energy_magn_j = hypot(E_re(j), E_im(j))
3317 : 0 : energy_magn_j = real(E_re(j))
3318 : 0 : energy_interval(j) = energy_magn_j
3319 : 0 : E_tot = E_tot + energy_magn_j
3320 : : end do
3321 [ # # ]: 0 : nvec_approx = ceiling(E_tot)
3322 : :
3323 : 0 : end subroutine computeFilterEnergy
3324 : : !!***
3325 : :
3326 [ # # # # : 0 : end module m_slice_cprj
# # # # #
# # # # #
# # ]
3327 : : !!***
|
