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