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