Line data Source code
1 : !!****f* ABINIT/m_chebfi2_cprj
2 : !! NAME
3 : !! m_chebfi2_cprj
4 : !!
5 : !! FUNCTION
6 : !! This module contains the types and routines used to apply the
7 : !! Chebyshev filtering method (2021 implementation using xG abstraction layer)
8 : !! It mainly defines a 'chebfi' datatypes and associated methods.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2023-2026 ABINIT group (LB)
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_chebfi2_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 :
33 : use m_cgtools
34 : use m_xg
35 : use m_xgTransposer
36 : use m_xg_ortho_RR
37 : use m_xg_nonlop
38 :
39 : use m_xmpi
40 : use m_xomp
41 : #ifdef HAVE_OPENMP
42 : use omp_lib
43 : #endif
44 :
45 : implicit none
46 :
47 : private
48 :
49 : !Several (private) parameters
50 : !-------------------------------------------------
51 :
52 : integer, parameter :: tim_init = 2061
53 : integer, parameter :: tim_free = 2062
54 : integer, parameter :: tim_cprj = 2063
55 : integer, parameter :: tim_invovl = 2065
56 : integer, parameter :: tim_residu = 2066
57 : integer, parameter :: tim_RR = 2067
58 : integer, parameter :: tim_transpose = 2068
59 : integer, parameter :: tim_RR_q = 2069
60 : integer, parameter :: tim_postinvovl = 2070
61 : integer, parameter :: tim_swap = 2071
62 : integer, parameter :: tim_amp_f = 2072
63 : integer, parameter :: tim_oracle = 2073
64 : integer, parameter :: tim_barrier = 2074
65 : integer, parameter :: tim_copy = 2075
66 : integer, parameter :: tim_ax_k = 2076
67 : integer, parameter :: tim_ax_v = 2077
68 : integer, parameter :: tim_ax_nl = 2078
69 : integer, parameter :: tim_enl = 2079
70 :
71 : !Public 'chebfi' datatype
72 : !-------------------------------------------------
73 :
74 : type, public :: chebfi_t
75 : integer :: space
76 : integer :: space_cprj
77 : integer :: spacedim ! Space dimension for one vector
78 : integer :: cprjdim ! cprj dimension
79 : integer :: total_spacedim ! Maybe not needed
80 : integer :: neigenpairs ! Number of eigen values/vectors we want
81 : integer :: ndeg_filter ! Degree of the polynomial filter
82 : integer :: nbdbuf ! Number of bands in the buffer
83 : integer :: spacecom ! Communicator for MPI
84 : integer :: oracle ! Option to compute ndeg_filter from residuals
85 : real(dp) :: tolerance ! Tolerance on the residu to stop the minimization
86 : real(dp) :: ecut ! Ecut for Chebfi oracle
87 : real(dp) :: oracle_factor ! factor used to decrease residuals
88 : real(dp) :: oracle_min_occ ! threshold on occupancies used for nbdbuf=-101
89 :
90 : integer :: bandpp
91 :
92 : logical :: paw
93 : integer :: eigenProblem !1 (A*x = (lambda)*B*x), 2 (A*B*x = (lambda)*x), 3 (B*A*x = (lambda)*x)
94 : integer :: me_g0
95 : integer :: me_g0_fft
96 :
97 : type(xg_nonlop_t) :: xg_nonlop
98 :
99 : !ARRAYS
100 : type(xgBlock_t) :: X
101 :
102 : type(xg_t) :: X_NP
103 : type(xgBlock_t) :: X_next
104 : type(xgBlock_t) :: X_prev
105 :
106 : type(xg_t) :: AX
107 : type(xgBlock_t) :: cprjX
108 : type(xg_t) :: cprj_work
109 : type(xg_t) :: cprj_work2
110 : type(xg_t) :: proj_work
111 :
112 : type(xgBlock_t) :: xXColsRows
113 : type(xgBlock_t) :: xAXColsRows
114 :
115 : type(xgTransposer_t) :: xgTransposerX
116 : type(xgTransposer_t) :: xgTransposerAX
117 :
118 : type(xgBlock_t) :: eigenvalues
119 :
120 : !SWAP POINTERS
121 : type(xgBlock_t) :: X_swap
122 : type(xgBlock_t) :: AX_swap
123 :
124 : end type chebfi_t
125 :
126 : !Public methods associated to 'chebfi' datatype
127 : !-------------------------------------------------
128 : public :: chebfi_init
129 : public :: chebfi_free
130 : public :: chebfi_memInfo
131 : public :: chebfi_run_cprj
132 :
133 : CONTAINS !========================================================================================
134 : !!***
135 :
136 : !!****f* m_chebfi2_cprj/chebfi_init
137 : !! NAME
138 : !! chebfi_init
139 : !!
140 : !! FUNCTION
141 : !! Initialize a 'chebfi' datastructure.
142 : !!
143 : !! INPUTS
144 : !! bandpp= number of 'bands' handled by a processor
145 : !! eigenProblem= type of eigenpb: 1 (A*x = (lambda)*B*x), 2 (A*B*x = (lambda)*x), 3 (B*A*x = (lambda)*x)
146 : !! me_g0= 1 if this processors treats G=0, 0 otherwise
147 : !! neigenpairs= number of requested eigenvectors/eigenvalues
148 : !! ndeg_filter= polynomial degree of the Chebyshev filter (.i.e. number of H applications)
149 : !! space= defines in which space we are (columns, rows, etc.)
150 : !! spacecom= MPI communicator
151 : !! spacedim= space dimension for one vector
152 : !! paw= flag. TRUE if current calculation ses the PAW approach
153 : !! ecut= plane-wave cut-off energy
154 : !! tolerance= tolerance criterion on the residu to stop the minimization
155 : !!
156 : !! OUTPUT
157 : !!
158 : !! SIDE EFFECTS
159 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
160 : !!
161 : !! SOURCE
162 :
163 15368 : subroutine chebfi_init(chebfi,neigenpairs,spacedim,cprjdim,tolerance,ecut,bandpp, &
164 : ndeg_filter,nbdbuf,space,space_cprj,eigenProblem,spacecom,me_g0,paw,&
165 : oracle,oracle_factor,oracle_min_occ,xg_nonlop,me_g0_fft)
166 :
167 : !Arguments ------------------------------------
168 : integer , intent(in ) :: bandpp
169 : integer , intent(in ) :: eigenProblem
170 : integer , intent(in ) :: me_g0
171 : integer , intent(in ) :: me_g0_fft
172 : integer , intent(in ) :: neigenpairs
173 : integer , intent(in ) :: ndeg_filter
174 : integer , intent(in ) :: nbdbuf
175 : integer , intent(in ) :: space
176 : integer , intent(in ) :: space_cprj
177 : integer , intent(in ) :: spacecom
178 : integer , intent(in ) :: spacedim
179 : integer , intent(in ) :: cprjdim
180 : integer , intent(in ) :: oracle
181 : logical , intent(in ) :: paw
182 : real(dp) , intent(in ) :: ecut
183 : real(dp) , intent(in ) :: tolerance
184 : real(dp) , intent(in ) :: oracle_factor
185 : real(dp) , intent(in ) :: oracle_min_occ
186 : type(xg_nonlop_t), intent(in ) :: xg_nonlop
187 : type(chebfi_t) , intent(inout) :: chebfi
188 :
189 : !Local variables-------------------------------
190 : real(dp) :: tsec(2)
191 :
192 : ! *********************************************************************
193 :
194 15368 : call timab(tim_init,1,tsec)
195 :
196 15368 : chebfi%space = space
197 15368 : chebfi%space_cprj = space_cprj
198 15368 : chebfi%neigenpairs = neigenpairs
199 15368 : chebfi%spacedim = spacedim
200 15368 : chebfi%cprjdim = cprjdim
201 15368 : if (tolerance > 0.0) then
202 2312 : chebfi%tolerance = tolerance
203 : else
204 13056 : chebfi%tolerance = 1.0e-20
205 : end if
206 15368 : chebfi%ecut = ecut
207 15368 : chebfi%bandpp = bandpp
208 15368 : chebfi%ndeg_filter = ndeg_filter
209 15368 : chebfi%nbdbuf = nbdbuf
210 15368 : chebfi%spacecom = spacecom
211 15368 : chebfi%eigenProblem = eigenProblem
212 15368 : chebfi%me_g0 = me_g0
213 15368 : chebfi%me_g0_fft = me_g0_fft
214 15368 : chebfi%paw = paw
215 15368 : chebfi%xg_nonlop = xg_nonlop
216 15368 : chebfi%oracle = oracle
217 15368 : chebfi%oracle_factor = oracle_factor
218 15368 : chebfi%oracle_min_occ = oracle_min_occ
219 :
220 15368 : call chebfi_allocateAll(chebfi)
221 :
222 15368 : call timab(tim_init,2,tsec)
223 :
224 15368 : end subroutine chebfi_init
225 : !!***
226 :
227 : !----------------------------------------------------------------------
228 :
229 : !!****f* m_chebfi2_cprj/chebfi_allocateAll
230 : !! NAME
231 : !! chebfi_allocateAll
232 : !!
233 : !! FUNCTION
234 : !! Allocate all memory spaces in a 'chebfi' datastructure.
235 : !!
236 : !! INPUTS
237 : !!
238 : !! OUTPUT
239 : !!
240 : !! SIDE EFFECTS
241 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
242 : !!
243 : !! SOURCE
244 :
245 30736 : subroutine chebfi_allocateAll(chebfi)
246 :
247 : ! Arguments ------------------------------------
248 : type(chebfi_t) , intent(inout) :: chebfi
249 :
250 : ! Local variables-------------------------------
251 : ! scalars
252 : integer :: neigenpairs
253 : integer :: space,space_cprj
254 : integer :: spacedim
255 : integer :: total_spacedim, ierr
256 : integer :: nspinor
257 :
258 : ! *********************************************************************
259 :
260 15368 : space = chebfi%space
261 15368 : space_cprj = chebfi%space_cprj
262 15368 : spacedim = chebfi%spacedim
263 15368 : neigenpairs = chebfi%neigenpairs
264 15368 : nspinor = chebfi%xg_nonlop%nspinor
265 :
266 15368 : call chebfi_free(chebfi)
267 :
268 15368 : total_spacedim = spacedim
269 15368 : call xmpi_sum(total_spacedim,chebfi%spacecom,ierr)
270 15368 : chebfi%total_spacedim = total_spacedim
271 15368 : call xg_init(chebfi%X_NP,space,total_spacedim,2*chebfi%bandpp,xmpi_comm_self,me_g0=chebfi%me_g0_fft) !transposed arrays
272 15368 : call xg_setBlock(chebfi%X_NP,chebfi%X_next,total_spacedim,chebfi%bandpp)
273 15368 : call xg_setBlock(chebfi%X_NP,chebfi%X_prev,total_spacedim,chebfi%bandpp,fcol=chebfi%bandpp+1)
274 :
275 15368 : call xg_init(chebfi%AX,space,spacedim,neigenpairs,chebfi%spacecom,me_g0=chebfi%me_g0)
276 15368 : call xg_init(chebfi%cprj_work ,space_cprj,chebfi%cprjdim,chebfi%bandpp*nspinor,chebfi%spacecom)
277 15368 : call xg_init(chebfi%cprj_work2,space_cprj,chebfi%cprjdim,chebfi%bandpp*nspinor,chebfi%spacecom)
278 :
279 15368 : call xg_init(chebfi%proj_work,space,chebfi%xg_nonlop%max_npw_k,chebfi%xg_nonlop%cprjdim,chebfi%spacecom,me_g0=chebfi%me_g0)
280 :
281 15368 : end subroutine chebfi_allocateAll
282 : !!***
283 :
284 : !----------------------------------------------------------------------
285 :
286 : !!****f* m_chebfi2_cprj/chebfi_free
287 : !! NAME
288 : !! chebfi_free
289 : !!
290 : !! FUNCTION
291 : !! Destroy a 'chebfi' datastructure.
292 : !!
293 : !! INPUTS
294 : !!
295 : !! OUTPUT
296 : !! arraymem(2)= memory information
297 : !!
298 : !! SIDE EFFECTS
299 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
300 : !!
301 : !! SOURCE
302 :
303 30736 : subroutine chebfi_free(chebfi)
304 :
305 : !Arguments ------------------------------------
306 : type(chebfi_t) , intent(inout) :: chebfi
307 :
308 : ! *********************************************************************
309 :
310 30736 : call xg_free(chebfi%X_NP)
311 :
312 30736 : call xg_free(chebfi%AX)
313 30736 : call xg_free(chebfi%cprj_work)
314 30736 : call xg_free(chebfi%cprj_work2)
315 30736 : call xg_free(chebfi%proj_work)
316 :
317 30736 : end subroutine chebfi_free
318 : !!***
319 :
320 : !----------------------------------------------------------------------
321 :
322 : !!****f* m_chebfi2_cprj/chebfi_memInfo
323 : !! NAME
324 : !! chebfi_memInfo
325 : !!
326 : !! FUNCTION
327 : !! Provides memory information about a 'chebfi' datastructure.
328 : !!
329 : !! INPUTS
330 : !! bandpp= number of 'bands' handled by a processor
331 : !! neigenpairs= number of requested eigenvectors/eigenvalues
332 : !! paral_kgb= flag controlling (k,g,bands) parallelization
333 : !! space= defines in which space we are (columns, rows, etc.)
334 : !! spacedim= dimension of MPI communicator
335 : !! total_spacedim= size of global KGB communicator (typically 'banspinorfft' comm.)
336 : !!
337 : !! OUTPUT
338 : !! arraymem(2)= memory information
339 : !!
340 : !! SIDE EFFECTS
341 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
342 : !!
343 : !! SOURCE
344 :
345 0 : function chebfi_memInfo(neigenpairs,spacedim,space,total_spacedim,bandpp) result(arraymem)
346 :
347 : !Arguments ------------------------------------
348 : integer, intent(in ) :: bandpp
349 : integer, intent(in ) :: neigenpairs
350 : integer, intent(in ) :: space
351 : integer, intent(in ) :: spacedim
352 : integer, intent(in ) :: total_spacedim
353 :
354 : !Local variables-------------------------------
355 : !scalars
356 : real(dp) :: memX
357 : real(dp) :: memX_next
358 : real(dp) :: memX_prev
359 : real(dp) :: memAX
360 : real(dp) :: memBX
361 : !Transposer variables
362 : real(dp) :: memX_CR
363 : real(dp) :: memAX_CR
364 : real(dp) :: memBX_CR
365 : !chebfi_rayleighRitz function variables
366 : real(dp) :: memA_und_X
367 : real(dp) :: memB_und_X
368 : real(dp) :: memEigenvalues
369 : real(dp) :: cplx
370 : !arrays
371 : real(dp) :: arraymem(2)
372 :
373 : ! *********************************************************************
374 0 : cplx = 1
375 0 : if ( space == SPACE_C ) cplx = 2 !for now only complex
376 :
377 : !Permanent in chebfi
378 0 : memX = cplx * kind(1.d0) * spacedim * neigenpairs
379 :
380 : ! if (paral_kgb == 0) then
381 : ! memX_next = cplx * kind(1.d0) * spacedim * neigenpairs
382 : ! memX_prev = cplx * kind(1.d0) * spacedim * neigenpairs
383 : ! else
384 0 : memX_next = cplx * kind(1.d0) * total_spacedim * bandpp
385 0 : memX_prev = cplx * kind(1.d0) * total_spacedim * bandpp
386 : ! end if
387 :
388 0 : memAX = cplx * kind(1.d0) * spacedim * neigenpairs
389 0 : memBX = cplx * kind(1.d0) * spacedim * neigenpairs
390 :
391 : !Transposer colrow array
392 : ! if (paral_kgb == 1) then
393 0 : memX_CR = cplx * kind(1.d0) * total_spacedim * bandpp
394 0 : memAX_CR = cplx * kind(1.d0) * total_spacedim * bandpp
395 0 : memBX_CR = cplx * kind(1.d0) * total_spacedim * bandpp
396 : ! else
397 : ! memX_CR = 0
398 : ! memAX_CR = 0
399 : ! memBX_CR = 0
400 : ! end if
401 :
402 : !chebfi_rayleighRitz function variables
403 0 : memA_und_X = cplx * kind(1.d0) * neigenpairs * neigenpairs
404 0 : memB_und_X = cplx * kind(1.d0) * neigenpairs * neigenpairs
405 0 : memEigenvalues = kind(1.d0) * neigenpairs
406 :
407 : arraymem(1) = memX + memX_next + memX_prev + &
408 0 : memAX + memBX + memX_CR + memAX_CR + memBX_CR
409 0 : arraymem(2) = memA_und_X + memB_und_X + memEigenvalues
410 :
411 0 : end function chebfi_memInfo
412 : !!***
413 :
414 : !----------------------------------------------------------------------
415 :
416 : !!****f* m_chebfi2_cprj/chebfi_run
417 : !! NAME
418 : !! chebfi_run
419 : !!
420 : !! FUNCTION
421 : !! Apply the Chebyshev Filtering algorithm on a set of vectors.
422 : !!
423 : !! INPUTS
424 : !! mpi_enreg = information about MPI parallelization
425 : !! getAX_BX= pointer to the function giving A|X> and B|X>
426 : !! A is typically the Hamiltonian H, and B the overlap operator S
427 : !! getBm1X= pointer to the function giving B^-1|X>
428 : !! B is typically the overlap operator S
429 : !!
430 : !! OUTPUT
431 : !!
432 : !! SIDE EFFECTS
433 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
434 : !! eigen= Full eigenvalues (initial values on entry)
435 : !! residu= residuals, i.e. norm of (A-lambdaB)|X>
436 : !! X0= Full set of vectors (initial values on entry)
437 : !!
438 : !! SOURCE
439 :
440 15368 : subroutine chebfi_run_cprj(chebfi,X0,cprjX0,getAX,kin,eigen,occ,residu,enl,nspinor)
441 :
442 : !Arguments ------------------------------------
443 : type(chebfi_t) , intent(inout) :: chebfi
444 : integer, intent(in) :: nspinor
445 : type(xgBlock_t), intent(inout) :: X0
446 : type(xgBlock_t), intent(inout) :: cprjX0
447 : type(xgBlock_t), intent(inout) :: eigen
448 : type(xgBlock_t), intent(in) :: occ
449 : type(xgBlock_t), intent(inout) :: residu
450 : type(xgBlock_t), intent(inout) :: enl
451 : type(xgBlock_t), intent(in ) :: kin
452 : interface
453 : subroutine getAX(X,AX)
454 : use m_xg, only : xgBlock_t
455 : type(xgBlock_t), intent(inout) :: X
456 : type(xgBlock_t), intent(inout) :: AX
457 : end subroutine getAX
458 : end interface
459 :
460 : !Local variables-------------------------------
461 : !scalars
462 : integer :: spacedim
463 : integer :: space_res
464 : integer :: neigenpairs
465 : integer :: ndeg_filter,ndeg_filter_max
466 : integer :: ideg, ierr
467 : real(dp) :: tolerance
468 : real(dp) :: maxeig, maxeig_global
469 : real(dp) :: mineig, mineig_global
470 : real(dp) :: lambda_minus
471 : real(dp) :: lambda_plus
472 : real(dp) :: one_over_r
473 : real(dp) :: two_over_r
474 : real(dp) :: center
475 : real(dp) :: radius
476 : type(xg_t) :: DivResults
477 : !arrays
478 : real(dp) :: tsec(2)
479 : !Pointers similar to old Chebfi
480 15368 : integer,allocatable :: ndeg_filter_bands(:) !Oracle variable
481 15368 : type(xg_nonlop_t) :: xg_nonlop
482 :
483 : ! *********************************************************************
484 :
485 15368 : spacedim = chebfi%spacedim
486 15368 : neigenpairs = chebfi%neigenpairs
487 15368 : ndeg_filter = chebfi%ndeg_filter
488 15368 : xg_nonlop = chebfi%xg_nonlop
489 15368 : chebfi%eigenvalues = eigen
490 :
491 46104 : ABI_MALLOC(ndeg_filter_bands,(chebfi%bandpp))
492 15368 : if (chebfi%space==SPACE_C) then
493 9112 : space_res = SPACE_C
494 6256 : else if (chebfi%space==SPACE_CR) then
495 6256 : space_res = SPACE_R
496 : else
497 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
498 : end if
499 15368 : call xg_init(DivResults, space_res, chebfi%bandpp, 1)
500 :
501 15368 : tolerance = chebfi%tolerance
502 15368 : lambda_plus = chebfi%ecut
503 15368 : chebfi%X = X0
504 15368 : chebfi%cprjX = cprjX0
505 :
506 : ! Transpose
507 15368 : call timab(tim_transpose,1,tsec)
508 : call xgTransposer_constructor(chebfi%xgTransposerX,chebfi%X,chebfi%xXColsRows,nspinor,&
509 15368 : STATE_LINALG,TRANS_ALL2ALL,xmpi_comm_self,chebfi%spacecom,0,0,chebfi%me_g0_fft)
510 :
511 15368 : call xgTransposer_copyConstructor(chebfi%xgTransposerAX,chebfi%xgTransposerX,chebfi%AX%self,chebfi%xAXColsRows,STATE_LINALG)
512 :
513 15368 : call xgTransposer_transpose(chebfi%xgTransposerX,STATE_COLSROWS)
514 15368 : chebfi%xgTransposerAX%state = STATE_COLSROWS
515 15368 : call timab(tim_transpose,2,tsec)
516 :
517 15368 : call timab(tim_cprj,1,tsec)
518 15368 : call xg_nonlop_getcprj(xg_nonlop,chebfi%xXColsRows,chebfi%cprjX,chebfi%proj_work%self)
519 15368 : call timab(tim_cprj,2,tsec)
520 15368 : call timab(tim_AX_v,1,tsec)
521 15368 : call getAX(chebfi%xXColsRows,chebfi%xAXColsRows)
522 15368 : call timab(tim_AX_v,2,tsec)
523 15368 : call timab(tim_AX_k,1,tsec)
524 15368 : call xgBlock_add_diag(chebfi%xXColsRows,kin,nspinor,chebfi%xAXColsRows)
525 15368 : call timab(tim_AX_k,2,tsec)
526 15368 : call timab(tim_AX_nl,1,tsec)
527 15368 : call xg_nonlop_getHX(xg_nonlop,chebfi%xAXcolsRows,chebfi%cprjX,chebfi%cprj_work%self,chebfi%proj_work%self)
528 15368 : call timab(tim_AX_nl,2,tsec)
529 :
530 15368 : call timab(tim_barrier,1,tsec)
531 15368 : call xmpi_barrier(chebfi%spacecom)
532 15368 : call timab(tim_barrier,2,tsec)
533 :
534 : !********************* Compute Rayleigh quotients for every band, and set lambda equal to the largest one *****
535 15368 : call timab(tim_RR_q, 1, tsec)
536 15368 : call chebfi_rayleighRitzQuotients(chebfi, maxeig, mineig, DivResults%self)
537 :
538 15368 : call xmpi_max(maxeig,maxeig_global,chebfi%spacecom,ierr)
539 15368 : call xmpi_min(mineig,mineig_global,chebfi%spacecom,ierr)
540 15368 : call timab(tim_RR_q, 2, tsec)
541 :
542 15368 : lambda_minus = maxeig_global
543 :
544 15368 : call timab(tim_oracle,1,tsec)
545 :
546 : ! ndeg_filter_max limits the reduction of the residual of the smallest eigenvalue (i.e. the most amplified one by the filter) by a factor 1e8.
547 : ! Also, the maximal value of ndeg_filter_max is 40.
548 15368 : ndeg_filter_max = cheb_oracle1(mineig_global, lambda_minus, lambda_plus, 1D-16, 40)
549 15368 : ndeg_filter = MIN(ndeg_filter_max,chebfi%ndeg_filter)
550 15368 : if (chebfi%oracle>0) then
551 2312 : call chebfi_set_ndeg_from_residu(chebfi,lambda_minus,lambda_plus,occ,DivResults%self,ndeg_filter_max,ndeg_filter)
552 : end if
553 96712 : ndeg_filter_bands(:) = ndeg_filter
554 :
555 15368 : call timab(tim_oracle,2,tsec)
556 :
557 15368 : center = (lambda_plus + lambda_minus)*0.5
558 15368 : radius = (lambda_plus - lambda_minus)*0.5
559 :
560 15368 : one_over_r = 1/radius
561 15368 : two_over_r = 2/radius
562 :
563 129554 : do ideg = 0, ndeg_filter - 1
564 :
565 114186 : call timab(tim_cprj,1,tsec)
566 114186 : call xg_nonlop_getcprj(xg_nonlop,chebfi%xAXcolsrows,chebfi%cprjX,chebfi%proj_work%self)
567 114186 : call timab(tim_cprj,2,tsec)
568 114186 : call chebfi_computeNextOrderChebfiPolynom(chebfi, ideg, center, one_over_r, two_over_r)
569 :
570 114186 : call timab(tim_swap,1,tsec)
571 114186 : call chebfi_swapInnerBuffers(chebfi, chebfi%total_spacedim, chebfi%bandpp)
572 114186 : call timab(tim_swap,2,tsec)
573 :
574 : !A * Psi
575 114186 : call timab(tim_AX_v,1,tsec)
576 114186 : call getAX(chebfi%xXColsRows,chebfi%xAXColsRows)
577 114186 : call timab(tim_AX_v,2,tsec)
578 114186 : call timab(tim_AX_k,1,tsec)
579 114186 : call xgBlock_add_diag(chebfi%xXColsRows,kin,nspinor,chebfi%xAXColsRows)
580 114186 : call timab(tim_AX_k,2,tsec)
581 114186 : call timab(tim_cprj,1,tsec)
582 114186 : call xg_nonlop_getcprj(xg_nonlop,chebfi%xXColsRows,chebfi%cprjX,chebfi%proj_work%self)
583 114186 : call timab(tim_cprj,2,tsec)
584 114186 : call timab(tim_AX_nl,1,tsec)
585 114186 : call xg_nonlop_getHX(xg_nonlop,chebfi%xAXcolsRows,chebfi%cprjX,chebfi%cprj_work%self,chebfi%proj_work%self)
586 129554 : call timab(tim_AX_nl,2,tsec)
587 :
588 : end do
589 :
590 15368 : call timab(tim_barrier,1,tsec)
591 15368 : call xmpi_barrier(chebfi%spacecom)
592 15368 : call timab(tim_barrier,2,tsec)
593 :
594 15368 : call timab(tim_amp_f,1,tsec)
595 15368 : call chebfi_ampfactor(chebfi, DivResults%self, lambda_minus, lambda_plus, ndeg_filter_bands)
596 : ! this results in higher condition number so avoid
597 : !call chebfi_ampfactorMax(chebfi, DivResults%self, lambda_minus, lambda_plus, ndeg_filter_bands)
598 15368 : call timab(tim_amp_f,2,tsec)
599 :
600 15368 : call xg_free(DivResults)
601 15368 : ABI_FREE(ndeg_filter_bands)
602 :
603 15368 : call timab(tim_transpose,1,tsec)
604 15368 : call xmpi_barrier(chebfi%spacecom)
605 :
606 15368 : call xgTransposer_transpose(chebfi%xgTransposerX,STATE_LINALG)
607 15368 : call xgTransposer_transpose(chebfi%xgTransposerAX,STATE_LINALG)
608 :
609 15368 : if (xmpi_comm_size(chebfi%spacecom) == 1) then !only one MPI proc reset buffers to right addresses (because of X-Xcolwise swaps)
610 4216 : call xgBlock_setBlock(chebfi%xXColsRows , chebfi%X , spacedim, neigenpairs)
611 4216 : call xgBlock_setBlock(chebfi%xAXColsRows, chebfi%AX%self, spacedim, neigenpairs)
612 : end if
613 15368 : call timab(tim_transpose,2,tsec)
614 :
615 15368 : call timab(tim_cprj,1,tsec)
616 15368 : call xg_nonlop_getcprj(xg_nonlop,chebfi%X,chebfi%cprjX,chebfi%cprj_work%self)
617 15368 : call timab(tim_cprj,2,tsec)
618 : call xg_RayleighRitz_cprj(chebfi%xg_nonlop,chebfi%X,chebfi%cprjX,chebfi%AX%self,chebfi%eigenvalues,ierr,0,&
619 15368 : tim_RR,ABI_GPU_DISABLED,solve_ax_bx=.true.)
620 :
621 15368 : if (chebfi%paw) then
622 12920 : call timab(tim_AX_nl,1,tsec)
623 : call xg_nonlop_getHmeSX(xg_nonlop,chebfi%X,chebfi%cprjX,chebfi%AX%self,chebfi%eigenvalues,chebfi%cprj_work%self,&
624 12920 : & chebfi%cprj_work2%self,no_H=.True.)
625 12920 : call timab(tim_AX_nl,2,tsec)
626 : end if
627 :
628 15368 : call timab(tim_residu, 1, tsec)
629 :
630 15368 : if (.not.chebfi%paw) then
631 2448 : call xgBlock_yxmax(chebfi%AX%self,chebfi%eigenvalues,chebfi%X)
632 : end if
633 :
634 15368 : call xgBlock_colwiseNorm2(chebfi%AX%self, residu)
635 15368 : call timab(tim_residu, 2, tsec)
636 :
637 15368 : call timab(tim_copy, 1, tsec)
638 15368 : call xgBlock_copy(chebfi%X,X0)
639 15368 : call timab(tim_copy, 2, tsec)
640 :
641 15368 : call xgTransposer_free(chebfi%xgTransposerX)
642 15368 : call xgTransposer_free(chebfi%xgTransposerAX)
643 :
644 15368 : if (.not.chebfi%paw) then
645 2448 : call timab(tim_enl,1,tsec)
646 2448 : call xg_nonlop_colwiseXHX(chebfi%xg_nonlop,chebfi%cprjX,chebfi%cprj_work%self,enl)
647 2448 : call timab(tim_enl,2,tsec)
648 : end if
649 :
650 46104 : end subroutine chebfi_run_cprj
651 : !!***
652 :
653 : !----------------------------------------------------------------------
654 :
655 : !!****f* m_chebfi2_cprj/chebfi_rayleighRitzQuotients
656 : !! NAME
657 : !! chebfi_rayleighRitzQuotients
658 : !!
659 : !! FUNCTION
660 : !! Compute the Rayleigh-Ritz quotients.
661 : !!
662 : !! INPUTS
663 : !!
664 : !! OUTPUT
665 : !!
666 : !! SIDE EFFECTS
667 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
668 : !! maxeig= highest eigenvalue
669 : !! mineig= lowest eigenvalue
670 : !! DivResults= Rayleigh-Ritz quotients
671 : !!
672 : !! SOURCE
673 :
674 15368 : subroutine chebfi_rayleighRitzQuotients(chebfi,maxeig,mineig,DivResults)
675 :
676 : !Arguments ------------------------------------
677 : real(dp), intent(inout) :: maxeig
678 : real(dp), intent(inout) :: mineig
679 : type(chebfi_t), intent(inout) :: chebfi
680 : type(xgBlock_t), intent(inout) :: DivResults
681 :
682 : !Local variables-------------------------------
683 : !scalars
684 : type(xg_t)::Results1
685 : type(xg_t)::Results2
686 : type(xg_t)::Results_work
687 : !arrays
688 : integer :: maxeig_pos(2)
689 : integer :: mineig_pos(2)
690 : integer :: space_res
691 :
692 : ! *********************************************************************
693 :
694 15368 : if (space(chebfi%xXcolsRows)==SPACE_C) then
695 9112 : space_res = SPACE_C
696 6256 : else if (space(chebfi%xXcolsRows)==SPACE_CR) then
697 6256 : space_res = SPACE_R
698 : else
699 0 : ABI_ERROR('space(X) should be SPACE_C or SPACE_CR')
700 : end if
701 15368 : call xg_init(Results1, space_res, chebfi%bandpp, 1)
702 15368 : call xg_init(Results2, space_res, chebfi%bandpp, 1)
703 :
704 15368 : call xgBlock_colwiseDotProduct(chebfi%xXColsRows,chebfi%xAXColsRows,Results1%self,comm_loc=xmpi_comm_null)
705 :
706 15368 : call xgBlock_colwiseDotProduct(chebfi%xXColsRows,chebfi%xXColsRows,Results2%self,comm_loc=xmpi_comm_null)
707 15368 : if (chebfi%xg_nonlop%paw) then
708 12920 : call xg_init(Results_work, space_res, chebfi%bandpp, 1)
709 12920 : call xg_nonlop_colwiseXAX(chebfi%xg_nonlop,chebfi%xg_nonlop%Sij%self,chebfi%cprjX,chebfi%cprj_work%self,Results_work%self)
710 12920 : call xgBlock_add(Results2%self,Results_work%self)
711 12920 : call xg_free(Results_work)
712 : end if
713 :
714 15368 : call xgBlock_colwiseDivision(Results1%self, Results2%self, DivResults, maxeig, maxeig_pos, mineig, mineig_pos)
715 :
716 15368 : call xg_free(Results1)
717 15368 : call xg_free(Results2)
718 :
719 15368 : end subroutine chebfi_rayleighRitzQuotients
720 : !!***
721 :
722 : !----------------------------------------------------------------------
723 :
724 : !!****f* m_chebfi2_cprj/chebfi_computeNextOrderChebfiPolynom
725 : !! NAME
726 : !! chebfi_computeNextOrderChebfiPolynom
727 : !!
728 : !! FUNCTION
729 : !! From P_n(B-^1.A)|X> (where P_n is the Chebyshev polynom of order n),
730 : !! computes P_n+1(B-^1.A)|X>
731 : !!
732 : !! INPUTS
733 : !! ideg=current degree of polynom
734 : !! center=filter center
735 : !! one_over_r,two_over_r=1/R, 2/R, R being the radius of the filter
736 : !! getBm1X= pointer to the function giving B^-1|X>
737 : !! B is typically the overlap operator S
738 : !!
739 : !! OUTPUT
740 : !!
741 : !! SIDE EFFECTS
742 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
743 : !!
744 : !! SOURCE
745 :
746 114186 : subroutine chebfi_computeNextOrderChebfiPolynom(chebfi,ideg,center,one_over_r,two_over_r)
747 :
748 : !Arguments ------------------------------------
749 : real(dp) , intent(in) :: center
750 : integer , intent(in) :: ideg
751 : real(dp) , intent(in) :: one_over_r
752 : real(dp) , intent(in) :: two_over_r
753 : type(chebfi_t) , intent(inout) :: chebfi
754 :
755 : !Local variables-------------------------------
756 : real(dp) :: tsec(2)
757 :
758 : ! *********************************************************************
759 :
760 114186 : call timab(tim_copy, 1, tsec)
761 114186 : call xgBlock_copy(chebfi%xAXColsRows,chebfi%X_next)
762 114186 : call timab(tim_copy, 2, tsec)
763 :
764 114186 : if (chebfi%paw) then
765 94950 : call timab(tim_invovl, 1, tsec)
766 : call xg_nonlop_getSm1X(chebfi%xg_nonlop,chebfi%X_next,chebfi%cprjX,&
767 94950 : & chebfi%cprj_work%self,chebfi%cprj_work2%self,chebfi%proj_work%self)
768 94950 : call timab(tim_invovl, 2, tsec)
769 : else
770 19236 : call timab(tim_copy, 1, tsec)
771 19236 : call xgBlock_copy(chebfi%xAXColsRows,chebfi%X_next)
772 19236 : call timab(tim_copy, 2, tsec)
773 : end if
774 :
775 114186 : call timab(tim_postinvovl, 1, tsec)
776 114186 : call xgBlock_scale(chebfi%xXColsRows, center, 1) !scale by center
777 :
778 : !(B-1 * A * Psi^i-1 - c * Psi^i-1)
779 114186 : call xgBlock_saxpy(chebfi%X_next, dble(-1.0), chebfi%xXColsRows)
780 :
781 : !Psi^i-1 = 1/c * Psi^i-1
782 114186 : call xgBlock_scale(chebfi%xXColsRows, 1/center, 1) !counter scale by 1/center
783 :
784 114186 : if (ideg == 0) then
785 14984 : call xgBlock_scale(chebfi%X_next, one_over_r, 1)
786 : else
787 99202 : call xgBlock_scale(chebfi%X_next, two_over_r, 1)
788 :
789 99202 : call xgBlock_saxpy(chebfi%X_next, dble(-1.0), chebfi%X_prev)
790 : end if
791 :
792 114186 : call timab(tim_postinvovl, 2, tsec)
793 :
794 114186 : end subroutine chebfi_computeNextOrderChebfiPolynom
795 : !!***
796 :
797 : !----------------------------------------------------------------------
798 :
799 : !!****f* m_chebfi2_cprj/chebfi_swapInnerBuffers
800 : !! NAME
801 : !! chebfi_swapInnerBuffers
802 : !!
803 : !! FUNCTION
804 : !! Swap buffers inside a 'chebfi' datastructure.
805 : !!
806 : !! INPUTS
807 : !! neigenpairs= number of requested eigenvectors/eigenvalues
808 : !! spacedim= space dimension for one vector
809 : !!
810 : !! OUTPUT
811 : !!
812 : !! SIDE EFFECTS
813 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
814 : !!
815 : !! SOURCE
816 :
817 114186 : subroutine chebfi_swapInnerBuffers(chebfi,spacedim,neigenpairs)
818 :
819 : ! Arguments ------------------------------------
820 : integer , intent(in ) :: spacedim
821 : integer , intent(in ) :: neigenpairs
822 : type(chebfi_t) , intent(inout) :: chebfi
823 :
824 : ! *********************************************************************
825 :
826 114186 : call xgBlock_setBlock(chebfi%X_prev, chebfi%X_swap, spacedim, neigenpairs) !X_swap = X_prev
827 114186 : call xgBlock_setBlock(chebfi%xXColsRows, chebfi%X_prev, spacedim, neigenpairs) !X_prev = xXColsRows
828 114186 : call xgBlock_setBlock(chebfi%X_next, chebfi%xXColsRows, spacedim, neigenpairs) !xXColsRows = X_next
829 114186 : call xgBlock_setBlock(chebfi%X_swap, chebfi%X_next, spacedim, neigenpairs) !X_next = X_swap
830 :
831 114186 : end subroutine chebfi_swapInnerBuffers
832 : !!***
833 :
834 : !----------------------------------------------------------------------
835 :
836 : !!****f* m_chebfi2_cprj/chebfi_ampfactor
837 : !! NAME
838 : !! chebfi_ampfactor
839 : !!
840 : !! FUNCTION
841 : !! Compute amplification factor
842 : !!
843 : !! INPUTS
844 : !! eig (:,:)= eigenvalues
845 : !! lambda_minus,lambda_plus=
846 : !! ndeg_filter_bands(:)= degree of Chebyshev polynomial filter for each band
847 : !!
848 : !! OUTPUT
849 : !!
850 : !! SIDE EFFECTS
851 : !! residu<type(xgBlock_t)>= vector of residuals
852 : !! chebfi <type(chebfi_t)>=all data used to apply Chebyshev Filtering algorithm
853 : !!
854 : !! SOURCE
855 :
856 15368 : subroutine chebfi_ampfactor(chebfi,DivResults,lambda_minus,lambda_plus,ndeg_filter_bands)
857 :
858 : ! Arguments ------------------------------------
859 : integer, intent(in ) :: ndeg_filter_bands(:)
860 : type(xgBlock_t), intent(in ) :: DivResults
861 : real(dp), intent(in ) :: lambda_minus
862 : real(dp), intent(in ) :: lambda_plus
863 : type(chebfi_t), intent(inout) :: chebfi
864 :
865 : ! Local variables-------------------------------
866 : ! scalars
867 : integer :: iband
868 : real(dp) :: ampfactor
869 : real(dp) :: eig_per_band
870 : type(xgBlock_t) :: X_part
871 : type(xgBlock_t) :: AX_part
872 15368 : real(dp),pointer :: eig(:,:)
873 :
874 : ! *********************************************************************
875 :
876 15368 : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=chebfi%bandpp)
877 :
878 96712 : do iband = 1, chebfi%bandpp
879 :
880 81344 : eig_per_band = eig(1,iband)
881 :
882 : !cheb_poly1(x, n, a, b)
883 81344 : ampfactor = cheb_poly1(eig_per_band, ndeg_filter_bands(iband), lambda_minus, lambda_plus)
884 :
885 81344 : if(abs(ampfactor) < 1e-3) ampfactor = 1e-3 !just in case, avoid amplifying too much
886 :
887 81344 : call xgBlock_setBlock(chebfi%xXColsRows, X_part, chebfi%total_spacedim, 1, fcol=iband)
888 81344 : call xgBlock_setBlock(chebfi%xAXColsRows, AX_part, chebfi%total_spacedim, 1, fcol=iband)
889 :
890 81344 : call xgBlock_scale(X_part, 1/ampfactor, 1)
891 96712 : call xgBlock_scale(AX_part, 1/ampfactor, 1)
892 :
893 : end do
894 :
895 15368 : end subroutine chebfi_ampfactor
896 : !!***
897 :
898 : subroutine chebfi_ampfactorMax(chebfi,DivResults,lambda_minus,lambda_plus,ndeg_filter_bands)
899 :
900 : ! Arguments ------------------------------------
901 : integer, intent(in ) :: ndeg_filter_bands(:)
902 : type(xgBlock_t), intent(in ) :: DivResults
903 : real(dp), intent(in ) :: lambda_minus
904 : real(dp), intent(in ) :: lambda_plus
905 : type(chebfi_t), intent(inout) :: chebfi
906 :
907 : ! Local variables-------------------------------
908 : ! scalars
909 : integer :: iband
910 : real(dp) :: ampfactor
911 : !type(xgBlock_t) :: X_part
912 : !type(xgBlock_t) :: AX_part
913 : real(dp),pointer :: eig(:,:)
914 :
915 : ! *********************************************************************
916 :
917 : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=cols(DivResults))
918 :
919 : !cheb_poly1(x, n, a, b)
920 : ampfactor = maxval( (/ (cheb_poly1(eig(1,iband), ndeg_filter_bands(iband), lambda_minus, lambda_plus),&
921 : iband=1,cols(DivResults)) /) )
922 :
923 : call xgBlock_scale(chebfi%xXColsRows, 1/ampfactor, 1)
924 : call xgBlock_scale(chebfi%xAXColsRows, 1/ampfactor, 1)
925 :
926 : end subroutine chebfi_ampfactorMax
927 : !!***
928 :
929 : !----------------------------------------------------------------------
930 :
931 : !!****f* m_chebfi2_cprj/chebfi_oracle1
932 : !! NAME
933 : !! chebfi_oracle1
934 : !!
935 : !! FUNCTION
936 : !! Compute order of Chebyshev polynom necessary to converge to a given tol
937 : !!
938 : !! INPUTS
939 : !! xx= input variable
940 : !! aa= left bound of the interval
941 : !! bb= right bound of the interval
942 : !! tol= needed precision
943 : !! nmax= max number of iterations
944 : !!
945 : !! OUTPUT
946 : !!
947 : !! SIDE EFFECTS
948 : !!
949 : !! SOURCE
950 :
951 22620 : function cheb_oracle1(xx,aa,bb,tol,nmax) result(nn)
952 :
953 : ! Arguments ------------------------------------
954 : integer :: nn
955 : integer, intent(in) :: nmax
956 : real(dp), intent(in) :: xx,aa,bb
957 : real(dp), intent(in) :: tol
958 :
959 : ! Local variables-------------------------------
960 : integer :: ii
961 : real(dp) :: yy,yim1,xred,temp
962 :
963 : ! *************************************************************************
964 :
965 22620 : xred = (xx-(aa+bb)/2)/(bb-aa)*2
966 22620 : yy = xred
967 22620 : yim1 = 1 !ONE
968 :
969 22620 : nn = nmax
970 22620 : if(1/(yy**2) < tol) then
971 : nn = 1
972 : else
973 422342 : do ii=2, nmax-1
974 421826 : temp = yy
975 421826 : yy = 2*xred*yy - yim1
976 421826 : yim1 = temp
977 422342 : if(1/(yy**2) < tol) then
978 : nn = ii
979 : exit
980 : end if
981 : end do
982 : end if
983 :
984 22620 : end function cheb_oracle1
985 : !!***
986 :
987 : !----------------------------------------------------------------------
988 :
989 : !!****f* m_chebfi2_cprj/chebfi_poly1
990 : !! NAME
991 : !! chebfi_poly1
992 : !!
993 : !! FUNCTION
994 : !! Compute Chebyshev polynomial???
995 : !!
996 : !! INPUTS
997 : !! xx= input variable
998 : !! aa= left bound of the interval
999 : !! bb= right bound of the interval
1000 : !! nn=
1001 : !!
1002 : !! OUTPUT
1003 : !!
1004 : !! SIDE EFFECTS
1005 : !!
1006 : !! SOURCE
1007 :
1008 81344 : function cheb_poly1(xx,nn,aa,bb) result(yy)
1009 :
1010 : ! Arguments ------------------------------------
1011 : integer, intent(in) :: nn
1012 : real(dp), intent(in) :: xx, aa, bb
1013 : real(dp) :: yy
1014 :
1015 : ! Local variables-------------------------------
1016 : integer :: ii
1017 : real(dp) :: xred,yim1,temp
1018 :
1019 : ! *************************************************************************
1020 :
1021 81344 : xred = (xx-(aa+bb)/2)/(bb-aa)*2
1022 81344 : yy = xred
1023 81344 : yim1 = 1
1024 611328 : do ii= 2, nn
1025 529984 : temp = yy
1026 529984 : yy = 2*xred*yy - yim1
1027 611328 : yim1 = temp
1028 : end do
1029 :
1030 81344 : end function cheb_poly1
1031 : !!***
1032 :
1033 : !!****f* m_chebfi2_cprj/chebfi_set_ndeg_from_residu
1034 : !! NAME
1035 : !! chebfi_set_ndeg_from_residu
1036 : !!
1037 : !! FUNCTION
1038 : !! Compute ndeg_filter using the oracle and residuals.
1039 : !!
1040 : !! INPUTS
1041 : !!
1042 : !! OUTPUT
1043 : !!
1044 : !! SIDE EFFECTS
1045 : !!
1046 : !! SOURCE
1047 :
1048 2312 : subroutine chebfi_set_ndeg_from_residu(chebfi,lambda_minus,lambda_plus,occ,DivResults,ndeg_filter_max,ndeg_filter)
1049 :
1050 : integer,intent(in) :: ndeg_filter_max
1051 : integer,intent(out) :: ndeg_filter
1052 : type(chebfi_t), intent(inout) :: chebfi
1053 : type(xgBlock_t), intent(in) :: occ
1054 : type(xgBlock_t), intent(in) :: DivResults
1055 : real(dp), intent(in) :: lambda_minus, lambda_plus
1056 :
1057 : logical :: test1,test2,test3
1058 : integer :: iband_tot,iband
1059 : integer :: bandpp,ierr,ndeg_filter_tolwfr,ndeg_filter_decrease,nbdbuf,ndeg_filter_all,shift
1060 2312 : integer,allocatable :: ndeg_filter_bands(:)
1061 : type(xgBlock_t) :: occBlock,occ_reshaped
1062 : type(xg_t) :: residu
1063 2312 : real(dp),pointer :: residu_(:,:),occ_(:,:)
1064 : real(dp) :: eig_iband,res_iband,occ_iband
1065 2312 : real(dp),pointer :: eig(:,:)
1066 :
1067 2312 : bandpp = chebfi%bandpp
1068 :
1069 : !Compute residu here for oracle, use X_next as a work space
1070 : ! X_next = S|Psi>
1071 2312 : call xgBlock_copy(chebfi%xXColsRows,chebfi%X_next)
1072 2312 : if (chebfi%paw) then
1073 2312 : call xg_nonlop_getSX(chebfi%xg_nonlop,chebfi%X_next,chebfi%cprjX,chebfi%cprj_work%self,chebfi%proj_work%self)
1074 : end if
1075 : ! X_next = - eig * S|Psi>
1076 2312 : call xgBlock_ymax(chebfi%X_next,DivResults,0,1)
1077 : ! X_next = H|Psi> - eig * S|Psi>
1078 2312 : call xgBlock_add(chebfi%X_next,chebfi%xAXColsRows)
1079 : ! resid = |X_next|^2
1080 2312 : call xg_init(residu,SPACE_R,bandpp,1)
1081 2312 : call xgBlock_colwiseNorm2(chebfi%X_next, residu%self,comm_loc=xmpi_comm_null)
1082 :
1083 2312 : occ_reshaped = occ
1084 2312 : shift=xmpi_comm_rank(chebfi%spacecom)*bandpp
1085 2312 : call xgBlock_reshape(occ_reshaped,1,chebfi%neigenpairs)
1086 2312 : call xgBlock_setBlock(occ_reshaped,occBlock,1,bandpp,fcol=1+shift)
1087 2312 : call xgBlock_reshape(occBlock,bandpp,1)
1088 2312 : if (chebfi%nbdbuf==-101) then
1089 952 : call xgBlock_apply_diag(residu%self,occBlock,1)
1090 : end if
1091 :
1092 6936 : ABI_MALLOC(ndeg_filter_bands,(bandpp))
1093 :
1094 : ! DivResults could be complex (with null imaginary part), so bandpp has to be in cols, not rows
1095 2312 : call xgBlock_reverseMap(DivResults,eig,rows=1,cols=bandpp)
1096 2312 : call xgBlock_reverseMap(residu%self,residu_,rows=1,cols=bandpp)
1097 2312 : call xgBlock_reverseMap(occBlock,occ_,rows=1,cols=bandpp)
1098 :
1099 2312 : if (chebfi%nbdbuf>0) then
1100 : nbdbuf = chebfi%nbdbuf
1101 952 : else if (chebfi%nbdbuf==-101) then
1102 952 : nbdbuf = 0
1103 : end if
1104 :
1105 12104 : do iband=1, bandpp
1106 9792 : eig_iband = eig(1,iband)
1107 9792 : res_iband = residu_(1,iband)
1108 9792 : occ_iband = occ_(1,iband)
1109 9792 : iband_tot = iband + shift
1110 9792 : test1 = res_iband<chebfi%tolerance ! band already converged
1111 9792 : test2 = iband_tot>chebfi%neigenpairs-nbdbuf ! band in the buffer
1112 9792 : test3 = chebfi%nbdbuf==-101.and.occ_iband<chebfi%oracle_min_occ ! occupancy is too low
1113 12104 : if (test1.or.test2.or.test3) then
1114 5750 : ndeg_filter_bands(iband) = 0
1115 : else
1116 : !ndeg_filter necessary to converge to tolerance
1117 4042 : ndeg_filter_tolwfr = cheb_oracle1(eig_iband, lambda_minus, lambda_plus, chebfi%tolerance / res_iband, 1000)
1118 4042 : if (chebfi%oracle==1) then
1119 832 : ndeg_filter_bands(iband) = MIN(ndeg_filter_max, ndeg_filter_tolwfr, chebfi%ndeg_filter)
1120 3210 : else if (chebfi%oracle==2) then
1121 : !ndeg_filter necessary to decrease residual by a constant factor
1122 3210 : ndeg_filter_decrease = cheb_oracle1(eig_iband, lambda_minus, lambda_plus, chebfi%oracle_factor, 15)
1123 3210 : ndeg_filter_bands(iband) = MIN(ndeg_filter_max, ndeg_filter_tolwfr, ndeg_filter_decrease)
1124 : else
1125 0 : ABI_ERROR('Wrong value for chebfi%oracle')
1126 : end if
1127 : end if
1128 : end do
1129 12104 : ndeg_filter = MAXVAL(ndeg_filter_bands)
1130 2312 : call xmpi_max(ndeg_filter,ndeg_filter_all,chebfi%spacecom,ierr)
1131 2312 : ndeg_filter=ndeg_filter_all
1132 :
1133 2312 : call xg_free(residu)
1134 2312 : ABI_FREE(ndeg_filter_bands)
1135 :
1136 4624 : end subroutine chebfi_set_ndeg_from_residu
1137 : !!***
1138 :
1139 : !!****f* m_chebfi2/jackson_lowpass_coeffs
1140 : !! NAME
1141 : !! jackson_lowpass_coeffs
1142 : !!
1143 : !! FUNCTION
1144 : !! Compute Jackson-damped Chebyshev coefficients for a lowpass step
1145 : !! on interval [lambda_min, lambda_max], degree M
1146 :
1147 : subroutine jackson_lowpass_coeffs(M, c)
1148 : implicit none
1149 : integer, intent(in) :: M
1150 : real(dp), intent(out) :: c(1:M+1)
1151 : integer :: k
1152 : real(dp) :: theta, g
1153 :
1154 : ! Chebyshev coefficients for step at left end (-1 in scaled coords)
1155 : theta = acos(-1.0d0) ! step at left edge
1156 : c(1) = theta / Pi ! k = 0
1157 :
1158 : do k = 1, M
1159 : c(k+1) = -2.0d0 / Pi * sin(k*theta) / k
1160 : end do
1161 :
1162 : ! Apply Jackson damping
1163 : do k = 0, M
1164 : g = ((M-k+1)*cos(Pi*k/(M+1)) + sin(Pi*k/(M+1))/tan(Pi/(M+1)))/(M+1)
1165 : c(k+1) = c(k+1) * g
1166 : end do
1167 :
1168 : end subroutine jackson_lowpass_coeffs
1169 : !!***
1170 :
1171 0 : end module m_chebfi2_cprj
1172 : !!***
|
