Branch data Line data Source code
1 : : !!****m* ABINIT/m_chebfi
2 : : !! NAME
3 : : !! m_chebfi
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2014-2026 ABINIT group (AL)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_chebfi
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_xmpi
27 : : use m_abicore
28 : : use m_abi_linalg
29 : : use m_rayleigh_ritz
30 : : use m_invovl
31 : : use m_dtset
32 : :
33 : : use defs_abitypes, only : mpi_type
34 : : use m_time, only : timab
35 : : use m_cgtools, only : dotprod_g
36 : : use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_get_ikpt
37 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_axpby, pawcprj_copy
38 : : use m_hamiltonian, only : gs_hamiltonian_type
39 : : use m_getghc, only : getghc
40 : : use m_prep_kgb, only : prep_getghc, prep_index_wavef_bandpp
41 : :
42 : : implicit none
43 : :
44 : : private
45 : : !!***
46 : :
47 : : public :: chebfi
48 : : !!***
49 : :
50 : : contains
51 : : !!***
52 : :
53 : :
54 : : !!****f* ABINIT/chebfi
55 : : !! NAME
56 : : !! chebfi
57 : : !!
58 : : !! FUNCTION
59 : : !! this routine updates the wave functions at a given k-point,
60 : : !! using the ChebFi method (see paper by A. Levitt and M. Torrent)
61 : : !!
62 : : !! INPUTS
63 : : !! dtset <type(dataset_type)>=all input variales for this dataset
64 : : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k
65 : : !! kinpw(npw)=(modified) kinetic energy for each plane wave (hartree)
66 : : !! mpi_enreg=information about MPI parallelization
67 : : !! nband=number of bands at this k point for that spin polarization
68 : : !! npw=number of plane waves at this k point
69 : : !! nspinor=number of plane waves at this k point
70 : : !! prtvol=control print volume and debugging output
71 : : !!
72 : : !! OUTPUT
73 : : !! eig(nband)=array for holding eigenvalues (hartree)
74 : : !! resid(nband)=residuals for each states
75 : : !! If gs_hamk%usepaw==1:
76 : : !! gsc(2,*)=<g|s|c> matrix elements (s=overlap)
77 : : !! If gs_hamk%usepaw==0
78 : : !! enlx(nband)=contribution from each band to nonlocal psp + potential Fock ACE part of total energy, at this k-point
79 : : !!
80 : : !! SIDE EFFECTS
81 : : !! cg(2,*)=updated wavefunctions
82 : : !!
83 : : !! NOTES
84 : : !! -- TODO --
85 : : !! Normev?
86 : : !! Ecutsm
87 : : !! nspinor 2
88 : : !! spinors parallelisation
89 : : !! fock
90 : : !! -- Performance --
91 : : !! Improve load balancing
92 : : !! Don't diagonalize converged eigenvectors, just orthogonalize
93 : : !! Maybe don't diagonalize so often (once every two outer iterations?)
94 : : !! Benchmark diagonalizations, choose np_slk
95 : : !! How to chose npfft?
96 : : !! Implement MINRES for invovl
97 : : !! -- LOBPCG --
98 : : !! Improve stability (see paper by Lehoucq Sorensen, maybe use bunch-kaufman factorizations?)
99 : : !!
100 : : !! SOURCE
101 : :
102 : 48 : subroutine chebfi(cg,dtset,eig,enlx,gs_hamk,gsc,kinpw,mpi_enreg,nband,npw,nspinor,prtvol,resid)
103 : :
104 : : !Arguments ------------------------------------
105 : : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
106 : : type(dataset_type),intent(in) :: dtset
107 : : type(mpi_type),intent(inout) :: mpi_enreg
108 : : integer,intent(in) :: nband,npw,prtvol,nspinor
109 : : real(dp),intent(inout), target :: cg(2,npw*nspinor*nband),gsc(2,npw*nspinor*nband)
110 : : real(dp),intent(in) :: kinpw(npw)
111 : : real(dp),intent(out) :: resid(nband)
112 : : real(dp),intent(out) :: enlx(nband)
113 : : real(dp),intent(out) :: eig(nband)
114 : :
115 : : !Local variables-------------------------------
116 : 96 : real(dp) :: pcon(npw)
117 : : real(dp) :: filter_low
118 : : real(dp) :: filter_center, filter_radius
119 : 144 : real(dp), dimension(2, npw*nspinor*nband), target :: ghc, gvnlxc
120 : 48 : real(dp), allocatable, dimension(:,:) :: cg_filter_next, cg_filter_prev, gsm1hc_filter, gsc_filter_prev, gsc_filter_next
121 : 48 : real(dp), allocatable, dimension(:,:), target :: cg_alltoall1,gsc_alltoall1,ghc_alltoall1,gvnlxc_alltoall1
122 : 48 : real(dp), allocatable, dimension(:,:), target :: cg_alltoall2,gsc_alltoall2,ghc_alltoall2,gvnlxc_alltoall2
123 : 48 : real(dp), pointer, dimension(:,:) :: cg_filter, gsc_filter, ghc_filter, gvnlxc_filter
124 : 96 : real(dp) :: resid_vec(2, npw*nspinor)
125 : : logical :: has_fock,paw
126 : : integer :: shift, shift_cg_loadbalanced
127 : : integer :: iband, ideg, ispinor
128 : : integer :: sij_opt, cpopt
129 : : real(dp) :: eval, tsec(2)
130 : : integer :: tim_getghc = 5, ierr
131 : : integer :: i
132 : 48 : integer, allocatable :: index_wavef_band(:)
133 : : real(dp) :: maxeig, mineig
134 : 48 : real(dp), allocatable :: resids_filter(:), residvec_filter(:,:)
135 : 48 : integer, allocatable :: ndeg_filter_bands(:)
136 : : integer :: iactive, nactive
137 : : real(dp) :: ampfactor
138 : : integer :: ndeg_filter_max, ndeg_filter_decrease, ndeg_filter_tolwfr
139 : : ! real(dp) :: load_imbalance
140 : : integer :: mcg
141 : : real(dp) :: dprod_r, dprod_i
142 : : character(len=500) :: message
143 : 96 : integer :: rdisplsloc(mpi_enreg%nproc_band), recvcountsloc(mpi_enreg%nproc_band)
144 : 96 : integer :: sdisplsloc(mpi_enreg%nproc_band), sendcountsloc(mpi_enreg%nproc_band)
145 : : integer :: ikpt_this_proc, npw_filter, nband_filter
146 : 48 : type(pawcprj_type), allocatable :: cwaveprj(:,:), cwaveprj_next(:,:), cwaveprj_prev(:,:)
147 : : ! integer :: ndeg_filter_total
148 : :
149 : : ! timers
150 : : integer, parameter :: timer_chebfi = 1600, timer_alltoall = 1601, timer_apply_inv_ovl = 1602, timer_rotation = 1603
151 : : integer, parameter :: timer_subdiago = 1604, timer_subham = 1605, timer_ortho = 1606, timer_getghc = 1607
152 : : integer, parameter :: timer_residuals = 1608, timer_update_eigen = 1609, timer_sync = 1610
153 : :
154 : : ! *************************************************************************
155 : :
156 : : !======================================================================================================
157 : : ! Initialize, transpose input cg if paral_kgb
158 : : !======================================================================================================
159 : 48 : call timab(timer_chebfi,1,tsec)
160 : :
161 : : !Initializations
162 : 48 : paw = gs_hamk%usepaw == 1
163 : 48 : has_fock=(associated(gs_hamk%fockcommon))
164 : 48 : mcg = npw*nspinor*nband
165 : :
166 : : ! Init pcon
167 [ + + ]: 3156 : pcon = (27+kinpw*(18+kinpw*(12+8*kinpw))) / (27+kinpw*(18+kinpw*(12+8*kinpw)) + 16*kinpw**4)
168 : :
169 [ + + + + : 223824 : ghc=zero; gvnlxc=zero
+ + + + ]
170 : :
171 : : ! Initialize the _filter pointers. Depending on paral_kgb, they might point to the actual arrays or to _alltoall variables
172 [ + - ]: 48 : if (dtset%paral_kgb == 1) then
173 : 48 : ikpt_this_proc = bandfft_kpt_get_ikpt()
174 : 48 : npw_filter = bandfft_kpt(ikpt_this_proc)%ndatarecv
175 : 48 : nband_filter = mpi_enreg%bandpp
176 : :
177 [ - + + - : 144 : ABI_MALLOC(cg_alltoall1, (2, npw_filter*nspinor*nband_filter))
+ - + - -
+ - + ]
178 [ - + + - : 96 : ABI_MALLOC(gsc_alltoall1, (2, npw_filter*nspinor*nband_filter))
- + - + ]
179 [ - + + - : 96 : ABI_MALLOC(ghc_alltoall1, (2, npw_filter*nspinor*nband_filter))
- + - + ]
180 [ - + + - : 96 : ABI_MALLOC(gvnlxc_alltoall1, (2, npw_filter*nspinor*nband_filter))
- + - + ]
181 [ - + + - : 96 : ABI_MALLOC(cg_alltoall2, (2, npw_filter*nspinor*nband_filter))
- + - + ]
182 [ - + + - : 96 : ABI_MALLOC(gsc_alltoall2, (2, npw_filter*nspinor*nband_filter))
- + - + ]
183 [ - + + - : 96 : ABI_MALLOC(ghc_alltoall2, (2, npw_filter*nspinor*nband_filter))
- + - + ]
184 [ - + + - : 96 : ABI_MALLOC(gvnlxc_alltoall2, (2, npw_filter*nspinor*nband_filter))
- + - + ]
185 : :
186 : : ! Init tranpose variables
187 [ + + ]: 144 : recvcountsloc = bandfft_kpt(ikpt_this_proc)%recvcounts * 2 * nspinor * mpi_enreg%bandpp
188 [ + + ]: 144 : rdisplsloc = bandfft_kpt(ikpt_this_proc)%rdispls * 2 * nspinor * mpi_enreg%bandpp
189 [ + + ]: 144 : sendcountsloc = bandfft_kpt(ikpt_this_proc)%sendcounts * 2 * nspinor
190 [ + + ]: 144 : sdisplsloc = bandfft_kpt(ikpt_this_proc)%sdispls * 2 * nspinor
191 : :
192 : : ! Load balancing, so that each processor has approximately the same number of converged and non-converged bands
193 : : ! for two procs, rearrange 1 2 3 4 5 6 as 1 4 2 5 3 6
194 : : !
195 : : ! trick to save memory: ghc has the necessary size, and will be overwritten afterwards anyway
196 : : #define cg_loadbalanced ghc
197 : 48 : shift = 0
198 [ + + ]: 144 : do i=1, mpi_enreg%nproc_band
199 [ + + ]: 720 : do iband=1, mpi_enreg%bandpp
200 : 576 : shift_cg_loadbalanced = (i-1 + (iband-1)*mpi_enreg%nproc_band)*npw*nspinor
201 [ + + + + ]: 112464 : cg_loadbalanced(:, shift+1:shift+npw*nspinor) = cg(:, shift_cg_loadbalanced+1:shift_cg_loadbalanced+npw*nspinor)
202 : 672 : shift = shift + npw*nspinor
203 : : end do
204 : : end do
205 : :
206 : : ! Transpose input cg into cg_alloall1. cg_alltoall1 is now (npw_filter, nband_filter)
207 : 48 : call timab(timer_alltoall, 1, tsec)
208 : : call xmpi_alltoallv(cg_loadbalanced,sendcountsloc,sdisplsloc,cg_alltoall1,&
209 : 48 : & recvcountsloc,rdisplsloc,mpi_enreg%comm_band,ierr)
210 : 48 : call timab(timer_alltoall, 2, tsec)
211 : : #undef cg_loadbalanced
212 : :
213 : : ! sort according to bandpp (from lobpcg, I don't fully understand what's going on but it works and it's fast)
214 : : call prep_index_wavef_bandpp(mpi_enreg%nproc_band,mpi_enreg%bandpp,&
215 : : & nspinor,bandfft_kpt(ikpt_this_proc)%ndatarecv,&
216 : : & bandfft_kpt(ikpt_this_proc)%recvcounts,bandfft_kpt(ikpt_this_proc)%rdispls,&
217 : 48 : & index_wavef_band)
218 : :
219 [ + + + + ]: 111936 : cg_alltoall2(:,:) = cg_alltoall1(:,index_wavef_band)
220 : :
221 : 48 : cg_filter => cg_alltoall2
222 : 48 : gsc_filter => gsc_alltoall2
223 : 48 : ghc_filter => ghc_alltoall2
224 : 48 : gvnlxc_filter => gvnlxc_alltoall2
225 : : else
226 : 0 : npw_filter = npw
227 : 0 : nband_filter = nband
228 : :
229 : 0 : cg_filter => cg
230 : 0 : gsc_filter => gsc
231 : 0 : ghc_filter => ghc
232 : 0 : gvnlxc_filter => gvnlxc
233 : : end if
234 : : ! from here to the next alltoall, all computation is done on _filter variables, agnostic
235 : : ! to whether it's nband x npw (paral_kgb == 0) or ndatarecv*bandpp (paral_kgb = 1)
236 : :
237 : : ! Allocate filter variables for the application of the Chebyshev polynomial
238 [ - + + - : 144 : ABI_MALLOC(cg_filter_next, (2, npw_filter*nspinor*nband_filter))
+ - + - -
+ - + ]
239 [ - + + - : 96 : ABI_MALLOC(cg_filter_prev, (2, npw_filter*nspinor*nband_filter))
- + - + ]
240 [ - + + - : 96 : ABI_MALLOC(gsc_filter_prev, (2, npw_filter*nspinor*nband_filter))
- + - + ]
241 [ - + + - : 96 : ABI_MALLOC(gsc_filter_next, (2, npw_filter*nspinor*nband_filter))
- + - + ]
242 [ - + + - : 96 : ABI_MALLOC(gsm1hc_filter, (2, npw_filter*nspinor*nband_filter))
- + - + ]
243 : :
244 : : ! PAW init
245 [ + - ]: 48 : if(paw) then
246 [ - + + - : 1056 : ABI_MALLOC(cwaveprj, (gs_hamk%natom,nspinor*nband_filter))
- + + - +
- + - - +
- + - + +
+ + + +
- ]
247 [ - + + - : 1056 : ABI_MALLOC(cwaveprj_next, (gs_hamk%natom,nspinor*nband_filter))
- + + - +
- + - - +
- + + + +
+ ]
248 [ - + + - : 1056 : ABI_MALLOC(cwaveprj_prev, (gs_hamk%natom,nspinor*nband_filter))
- + + - +
- + - - +
- + + + +
+ ]
249 : 48 : call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
250 : 48 : call pawcprj_alloc(cwaveprj_next,0,gs_hamk%dimcprj)
251 : 48 : call pawcprj_alloc(cwaveprj_prev,0,gs_hamk%dimcprj)
252 : :
253 : 48 : sij_opt = 1 ! recompute S
254 : 48 : cpopt = 0 ! save cprojs
255 : : else
256 : 0 : sij_opt = 0
257 : 0 : cpopt = -1
258 : : end if
259 : :
260 : :
261 : :
262 : : !======================================================================================================
263 : : ! Data in npfft x npband distribution. First getghc, update eigenvalues and residuals
264 : : !======================================================================================================
265 : 48 : write(message, *) 'First getghc'
266 : 48 : call wrtout(std_out,message,'COLL')
267 : :
268 : : ! get_ghc on cg
269 : 48 : call timab(timer_getghc, 1, tsec)
270 [ - + ]: 48 : if (dtset%paral_kgb == 0) then
271 : : call getghc(cpopt,cg_filter,cwaveprj,ghc_filter,gsc_filter,gs_hamk,gvnlxc_filter,&
272 : 0 : & eval,mpi_enreg,nband,prtvol,sij_opt,tim_getghc,0)
273 : : else
274 : : call prep_getghc(cg_filter,gs_hamk,gvnlxc_filter,ghc_filter,gsc_filter,eval,nband,mpi_enreg,&
275 : 48 : & prtvol,sij_opt,cpopt,cwaveprj,already_transposed=.true.)
276 : : end if
277 : 48 : call timab(timer_getghc, 2, tsec)
278 : :
279 : : ! Debug barrier: should be invisible
280 : 48 : call timab(timer_sync, 1, tsec)
281 : 48 : call xmpi_barrier(mpi_enreg%comm_band)
282 : 48 : call timab(timer_sync, 2, tsec)
283 : :
284 : 48 : write(message, *) 'Computing residuals'
285 : 48 : call wrtout(std_out,message,'COLL')
286 : : ! update eigenvalues and residuals
287 : 48 : call timab(timer_update_eigen, 1, tsec)
288 [ - + + - : 144 : ABI_MALLOC(resids_filter, (nband_filter))
+ - + - -
+ - + ]
289 [ - + + - : 144 : ABI_MALLOC(residvec_filter, (2, npw_filter*nspinor))
+ - + - -
+ - + ]
290 [ - + + - : 144 : ABI_MALLOC(ndeg_filter_bands, (nband_filter))
+ - + - -
+ - + ]
291 [ + + ]: 336 : do iband=1, nband_filter
292 : 288 : shift = npw_filter*nspinor*(iband-1)
293 : : call dotprod_g(eig(iband),dprod_i,gs_hamk%istwf_k,npw_filter*nspinor,1,ghc_filter(:, shift+1:shift+npw_filter*nspinor),&
294 [ - + - + ]: 288 : & cg_filter(:, shift+1:shift+npw_filter*nspinor),mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
295 [ - + ]: 288 : if(paw) then
296 : : call dotprod_g(dprod_r,dprod_i,gs_hamk%istwf_k,npw_filter*nspinor,1,gsc_filter(:, shift+1:shift+npw_filter*nspinor),&
297 [ - + - + ]: 288 : & cg_filter(:, shift+1:shift+npw_filter*nspinor),mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
298 : 288 : eig(iband) = eig(iband)/dprod_r
299 : : end if
300 : :
301 : : if(paw) then
302 : : residvec_filter = ghc_filter(:, shift+1 : shift+npw_filter*nspinor) &
303 [ + - - + : 112464 : & - eig(iband)*gsc_filter(:, shift+1 : shift+npw_filter*nspinor)
- - + + +
+ ]
304 : : else
305 : : residvec_filter = ghc_filter(:, shift+1 : shift+npw_filter*nspinor) &
306 [ # # # # : 0 : & - eig(iband)*cg_filter(:, shift+1 : shift+npw_filter*nspinor)
# # # # #
# ]
307 : : end if
308 [ + + + + ]: 112224 : resids_filter(iband) = SUM(residvec_filter**2)
309 : : end do
310 : 48 : call xmpi_sum(resids_filter,mpi_enreg%comm_fft,ierr)
311 [ + - - + : 384 : call xmpi_max(MAXVAL(eig(1:nband_filter)),maxeig,mpi_enreg%comm_band,ierr)
- - + + ]
312 [ + - - + : 384 : call xmpi_min(MINVAL(eig(1:nband_filter)),mineig,mpi_enreg%comm_band,ierr)
- - + + ]
313 : 48 : filter_low = maxeig
314 : 48 : call timab(timer_update_eigen, 2, tsec)
315 : :
316 : : ! Decide which polynomial filter degree per band is needed
317 : : ! don't go above this, or face bad conditioning of the Gram matrix.
318 : 48 : ndeg_filter_max = cheb_oracle(mineig, filter_low, dtset%ecut, 1e-16_dp, 40)
319 [ + + ]: 336 : do iband=1, nband_filter
320 : : ! Filter degree necessary to converge to tolwfr
321 : 288 : ndeg_filter_tolwfr = cheb_oracle(eig(iband), filter_low, dtset%ecut, dtset%tolwfr_diago / resids_filter(iband), dtset%mdeg_filter)
322 : : ! Filter degree necessary to decrease residual by a constant factor
323 : : ndeg_filter_decrease = cheb_oracle(eig(iband), filter_low, dtset%ecut, 0.1_dp, dtset%mdeg_filter)
324 : :
325 : 288 : ndeg_filter_bands(iband) = MAX(MIN(ndeg_filter_tolwfr, ndeg_filter_decrease, ndeg_filter_max, dtset%mdeg_filter), 1)
326 : 336 : ndeg_filter_bands(iband) = dtset%mdeg_filter ! fiddle with this to use locking
327 : : end do
328 : :
329 : :
330 : : !!!!! Uncomment for diagnostics
331 : : ! ndeg_filter_total = SUM(ndeg_filter_bands)
332 : : ! call xmpi_sum(ndeg_filter_total, mpi_enreg%comm_band, ierr)
333 : : ! load_imbalance = (SUM(ndeg_filter_bands) - REAL(ndeg_filter_total)/REAL(mpi_enreg%nproc_band)) / &
334 : : ! & (REAL(ndeg_filter_total)/REAL(mpi_enreg%nproc_band))
335 : : ! call xmax_mpi(load_imbalance, mpi_enreg%comm_band, ierr)
336 : :
337 : : ! write(message, *) 'Mean ndeg_filter', REAL(ndeg_filter_total)/REAL(nband), 'max imbalance (%)', load_imbalance*100
338 : : ! call wrtout(std_out,message,'COLL')
339 : :
340 [ - + ]: 48 : ABI_FREE(resids_filter)
341 [ - + ]: 48 : ABI_FREE(residvec_filter)
342 : :
343 : : !======================================================================================================
344 : : ! Chebyshev polynomial application
345 : : !======================================================================================================
346 : : ! Filter by a chebyshev polynomial of max. degree mdeg_filter
347 [ + + ]: 336 : do ideg=1,dtset%mdeg_filter
348 : : ! Filter only on [iactive, iactive+nactive-1]
349 : 288 : iactive = nband_filter
350 [ + - ]: 288 : do iband = 1, nband_filter
351 : : ! does iband need a higher degree for polynomila filter?
352 [ - + ]: 288 : if (ndeg_filter_bands(iband) >= ideg) then
353 : : iactive = iband
354 : : exit
355 : : end if
356 : : end do
357 : 288 : nactive = nband_filter - iactive + 1
358 : 288 : shift = npw_filter*nspinor*(iactive-1) + 1
359 : : ! trick the legacy prep_getghc
360 : 288 : mpi_enreg%bandpp = nactive
361 : :
362 : : ! Define the filter position
363 : 288 : filter_center = (dtset%ecut+filter_low)/2
364 : 288 : filter_radius = (dtset%ecut-filter_low)/2
365 : :
366 : : ! write(message, *) 'Applying invovl, iteration', ideg
367 : : ! call wrtout(std_out,message,'COLL')
368 : :
369 : : ! If paw, have to apply S^-1
370 [ + - ]: 288 : if(paw) then
371 : 288 : call timab(timer_apply_inv_ovl, 1, tsec)
372 : : call apply_invovl(gs_hamk, ghc_filter(:,shift:), gsm1hc_filter(:,shift:), cwaveprj_next(:,iactive:), &
373 [ - + ]: 288 : & npw_filter, nactive, mpi_enreg, nspinor, dtset%invovl_blksliced, dtset%invovl_precision)
374 : 288 : call timab(timer_apply_inv_ovl, 2, tsec)
375 : : else
376 [ # # # # ]: 0 : gsm1hc_filter(:,shift:) = ghc_filter(:,shift:)
377 : : end if
378 : :
379 : : ! Chebyshev iteration: UPDATE cg
380 [ + + ]: 288 : if(ideg == 1) then
381 [ + + + + ]: 111936 : cg_filter_next(:,shift:) = one/filter_radius * (gsm1hc_filter(:,shift:) - filter_center*cg_filter(:,shift:))
382 : : else
383 : : cg_filter_next(:,shift:) = two/filter_radius * (gsm1hc_filter(:,shift:) - filter_center*cg_filter(:,shift:)) &
384 [ + + + + ]: 559680 : & - cg_filter_prev(:,shift:)
385 : : end if
386 : : ! Update gsc and cwaveprj
387 [ + - ]: 288 : if(paw) then
388 [ + + ]: 288 : if(ideg == 1) then
389 [ + + + + ]: 111936 : gsc_filter_next(:,shift:) = one/filter_radius * (ghc_filter(:,shift:) - filter_center*gsc_filter(:,shift:))
390 : : !cwaveprj_next = one/filter_radius * (cwaveprj_next - filter_center*cwaveprj)
391 : 48 : call pawcprj_axpby(-filter_center/filter_radius, one/filter_radius,cwaveprj(:,iactive:),cwaveprj_next(:,iactive:))
392 : : else
393 : : gsc_filter_next(:,shift:) = two/filter_radius * (ghc_filter(:,shift:) - filter_center*gsc_filter(:,shift:))&
394 [ + + + + ]: 559680 : & - gsc_filter_prev(:,shift:)
395 : : !cwaveprj_next = two/filter_radius * (cwaveprj_next - filter_center*cwaveprj) - cwaveprj_prev
396 : 240 : call pawcprj_axpby(-two*filter_center/filter_radius, two/filter_radius,cwaveprj(:,iactive:),cwaveprj_next(:,iactive:))
397 : 240 : call pawcprj_axpby(-one, one,cwaveprj_prev(:,iactive:),cwaveprj_next(:,iactive:))
398 : : end if
399 : : end if
400 : :
401 : : ! Bookkeeping of the _prev variables
402 [ + + + + ]: 671616 : cg_filter_prev(:,shift:) = cg_filter(:,shift:)
403 [ + + + + ]: 671616 : cg_filter(:,shift:) = cg_filter_next(:,shift:)
404 [ - + ]: 288 : if(paw) then
405 [ + + + + ]: 671616 : gsc_filter_prev(:,shift:) = gsc_filter(:,shift:)
406 [ + + + + ]: 671616 : gsc_filter(:,shift:) = gsc_filter_next(:,shift:)
407 : :
408 : 288 : call pawcprj_copy(cwaveprj(:,iactive:),cwaveprj_prev(:,iactive:))
409 : 288 : call pawcprj_copy(cwaveprj_next(:,iactive:),cwaveprj(:,iactive:))
410 : : end if
411 : :
412 : : ! Update ghc
413 : : if(paw) then
414 : : !! DEBUG use this to remove the optimization and recompute gsc/cprojs
415 : : ! sij_opt = 1
416 : : ! cpopt = 0
417 : :
418 : 288 : sij_opt = 0 ! gsc is already computed
419 : 288 : cpopt = 2 ! reuse cprojs
420 : : else
421 : 0 : sij_opt = 0
422 : 0 : cpopt = -1
423 : : end if
424 : :
425 : 288 : write(message, *) 'Getghc, degree of polynom', ideg
426 : 288 : call wrtout(std_out,message,'COLL')
427 : :
428 : 288 : call timab(timer_getghc, 1, tsec)
429 [ - + ]: 288 : if (dtset%paral_kgb == 0) then
430 : : call getghc(cpopt,cg_filter(:,shift:),cwaveprj(:,iactive:),ghc_filter(:,shift:),&
431 : : & gsc_filter(:,shift:),gs_hamk,gvnlxc_filter(:,shift:),eval,mpi_enreg,&
432 : 0 : & nband,prtvol,sij_opt,tim_getghc,0)
433 : : else
434 : : call prep_getghc(cg_filter(:,shift:),gs_hamk,gvnlxc_filter(:,shift:),ghc_filter(:,shift:),&
435 : : & gsc_filter(:,shift:),eval,nband,mpi_enreg,prtvol,sij_opt,cpopt,&
436 : 288 : & cwaveprj(:,iactive:),already_transposed=.true.)
437 : : end if
438 : :
439 : : ! end of the trick
440 : 288 : mpi_enreg%bandpp = nband_filter
441 : :
442 : 336 : call timab(timer_getghc, 2, tsec)
443 : : end do ! end loop on ideg
444 : :
445 : : ! normalize according to the previously computed rayleigh quotients (inaccurate, but cheap)
446 [ + + ]: 336 : do iband = 1, nband_filter
447 : 288 : ampfactor = cheb_poly(eig(iband), ndeg_filter_bands(iband), filter_low, dtset%ecut)
448 [ - + ]: 288 : if(abs(ampfactor) < 1e-3) ampfactor = 1e-3 ! just in case, avoid amplifying too much
449 : 288 : shift = npw_filter*nspinor*(iband-1)
450 [ + + + + ]: 112176 : cg_filter(:, shift+1:shift+npw_filter*nspinor) = cg_filter(:, shift+1:shift+npw_filter*nspinor) / ampfactor
451 [ + + + + ]: 112176 : ghc_filter(:, shift+1:shift+npw_filter*nspinor) = ghc_filter(:, shift+1:shift+npw_filter*nspinor) / ampfactor
452 [ + - ]: 288 : if(paw) then
453 [ + + + + ]: 112176 : gsc_filter(:, shift+1:shift+npw_filter*nspinor) = gsc_filter(:, shift+1:shift+npw_filter*nspinor) / ampfactor
454 : : endif
455 [ - + ]: 336 : if(.not.paw .or. has_fock)then
456 [ # # # # ]: 0 : gvnlxc_filter(:, shift+1:shift+npw_filter*nspinor) = gvnlxc_filter(:, shift+1:shift+npw_filter*nspinor) / ampfactor
457 : : end if
458 : : end do
459 : :
460 : : ! Cleanup
461 [ + - ]: 48 : if(paw) then
462 : 48 : call pawcprj_free(cwaveprj)
463 : 48 : call pawcprj_free(cwaveprj_next)
464 : 48 : call pawcprj_free(cwaveprj_prev)
465 [ + - + + : 624 : ABI_FREE(cwaveprj)
- + - + -
+ ]
466 [ + - + + : 624 : ABI_FREE(cwaveprj_next)
- + - + -
+ ]
467 [ + - + + : 624 : ABI_FREE(cwaveprj_prev)
- + - + -
+ ]
468 : : end if
469 [ - + ]: 48 : ABI_FREE(ndeg_filter_bands)
470 [ - + ]: 48 : ABI_FREE(cg_filter_next)
471 [ - + ]: 48 : ABI_FREE(cg_filter_prev)
472 [ - + ]: 48 : ABI_FREE(gsc_filter_prev)
473 [ - + ]: 48 : ABI_FREE(gsc_filter_next)
474 [ - + ]: 48 : ABI_FREE(gsm1hc_filter)
475 : :
476 : : !======================================================================================================
477 : : ! Filtering done, tranpose back
478 : : !======================================================================================================
479 : :
480 : 48 : write(message, *) 'Filtering done, transposing back'
481 : 48 : call wrtout(std_out,message,'COLL')
482 : :
483 : : ! transpose back
484 [ - + ]: 48 : if(dtset%paral_kgb == 1) then
485 [ + + + + ]: 111936 : cg_alltoall1(:,index_wavef_band) = cg_alltoall2(:,:)
486 [ + + + + ]: 111936 : ghc_alltoall1(:,index_wavef_band) = ghc_alltoall2(:,:)
487 [ + - ]: 48 : if(paw) then
488 [ + + + + ]: 111936 : gsc_alltoall1(:,index_wavef_band) = gsc_alltoall2(:,:)
489 : : else
490 [ # # # # ]: 0 : gvnlxc_alltoall1(:,index_wavef_band) = gvnlxc_alltoall2(:,:)
491 : : end if
492 : :
493 [ - + ]: 48 : ABI_FREE(index_wavef_band)
494 : :
495 : 48 : call timab(timer_sync, 1, tsec)
496 : 48 : call xmpi_barrier(mpi_enreg%comm_band)
497 : 48 : call timab(timer_sync, 2, tsec)
498 : :
499 : 48 : call timab(timer_alltoall, 1, tsec)
500 : :
501 : : ! Do we pack the arrays in the alltoall, saving latency, or do we do it separately, saving memory and copies?
502 : : call xmpi_alltoallv(cg_alltoall1,recvcountsloc,rdisplsloc,cg,&
503 : 48 : & sendcountsloc,sdisplsloc,mpi_enreg%comm_band,ierr)
504 : : call xmpi_alltoallv(ghc_alltoall1,recvcountsloc,rdisplsloc,ghc,&
505 : 48 : & sendcountsloc,sdisplsloc,mpi_enreg%comm_band,ierr)
506 [ + - ]: 48 : if(paw) then
507 : : call xmpi_alltoallv(gsc_alltoall1,recvcountsloc,rdisplsloc,gsc,&
508 : 48 : & sendcountsloc,sdisplsloc,mpi_enreg%comm_band,ierr)
509 : : else
510 : : call xmpi_alltoallv(gvnlxc_alltoall1,recvcountsloc,rdisplsloc,gvnlxc,&
511 : 0 : & sendcountsloc,sdisplsloc,mpi_enreg%comm_band,ierr)
512 : : end if
513 : 48 : call timab(timer_alltoall, 2, tsec)
514 : :
515 [ - + ]: 48 : if(mpi_enreg%paral_kgb == 1) then
516 [ - + ]: 48 : ABI_FREE(cg_alltoall1)
517 [ - + ]: 48 : ABI_FREE(gsc_alltoall1)
518 [ - + ]: 48 : ABI_FREE(ghc_alltoall1)
519 [ - + ]: 48 : ABI_FREE(gvnlxc_alltoall1)
520 [ - + ]: 48 : ABI_FREE(cg_alltoall2)
521 [ - + ]: 48 : ABI_FREE(gsc_alltoall2)
522 [ - + ]: 48 : ABI_FREE(ghc_alltoall2)
523 [ - + ]: 48 : ABI_FREE(gvnlxc_alltoall2)
524 : : end if
525 : : else
526 : : ! nothing to do, the _filter variables already point to the right ones
527 : : end if
528 : :
529 : :
530 : :
531 : : !======================================================================================================
532 : : ! Data in (npfft*npband) x 1 distribution. Rayleigh-Ritz step
533 : : !======================================================================================================
534 : :
535 : : ! _subdiago might use less memory when using only one proc, should maybe call it, or just remove it
536 : : ! and always call _distributed
537 : : #if defined HAVE_LINALG_SCALAPACK
538 : : call rayleigh_ritz_distributed(cg,ghc,gsc,gvnlxc,eig,has_fock,gs_hamk%istwf_k,mpi_enreg,nband,npw,nspinor,gs_hamk%usepaw)
539 : : #else
540 : 48 : call rayleigh_ritz_subdiago(cg,ghc,gsc,gvnlxc,eig,has_fock,gs_hamk%istwf_k,mpi_enreg,nband,npw,nspinor,gs_hamk%usepaw)
541 : : #endif
542 : :
543 : : ! Build residuals
544 : 48 : call timab(timer_residuals, 1, tsec)
545 [ + + ]: 624 : do iband=1,nband
546 : 576 : shift = npw*nspinor*(iband-1)
547 [ + - ]: 576 : if(paw) then
548 [ + + + + ]: 112464 : resid_vec = ghc(:, shift+1 : shift+npw*nspinor) - eig(iband)*gsc(:, shift+1 : shift+npw*nspinor)
549 : : else
550 [ # # # # ]: 0 : resid_vec = ghc(:, shift+1 : shift+npw*nspinor) - eig(iband)*cg (:, shift+1 : shift+npw*nspinor)
551 : : end if
552 : :
553 : : ! precondition resid_vec
554 [ + + ]: 1152 : do ispinor = 1,nspinor
555 [ + + ]: 37872 : resid_vec(1, npw*(ispinor-1)+1:npw*ispinor) = resid_vec(1, npw*(ispinor-1)+1:npw*ispinor) * pcon
556 [ + + ]: 38448 : resid_vec(2, npw*(ispinor-1)+1:npw*ispinor) = resid_vec(2, npw*(ispinor-1)+1:npw*ispinor) * pcon
557 : : end do
558 : :
559 : : call dotprod_g(resid(iband),dprod_i,gs_hamk%istwf_k,npw*nspinor,1,resid_vec,&
560 : 576 : & resid_vec,mpi_enreg%me_g0,mpi_enreg%comm_bandspinorfft)
561 : :
562 [ + - - + ]: 624 : if(.not. paw .or. has_fock) then
563 : : call dotprod_g(enlx(iband),dprod_i,gs_hamk%istwf_k,npw*nspinor,1,cg(:, shift+1:shift+npw*nspinor),&
564 : 0 : & gvnlxc(:, shift+1:shift+npw_filter*nspinor),mpi_enreg%me_g0,mpi_enreg%comm_bandspinorfft)
565 : : end if
566 : : end do
567 : 48 : call timab(timer_residuals, 2, tsec)
568 : :
569 : : ! write(message, '(a,4e10.2)') 'Resids (1, N, min, max) ', resid(1), resid(nband), MINVAL(resid), MAXVAL(resid)
570 : : ! call wrtout(std_out,message,'COLL')
571 : :
572 : : ! write(message,*)'Eigens(1,nocc,nband) ',eig(1), eig(ilastocc),eig(nband)
573 : : ! call wrtout(std_out,message,'COLL')
574 : : ! write(message,*)'Resids(1,nocc,nband) ',resid(1), resid(ilastocc),resid(nband)
575 : : ! call wrtout(std_out,message,'COLL')
576 : :
577 : 48 : call timab(timer_chebfi,2,tsec)
578 : :
579 [ - + - + : 144 : end subroutine chebfi
- + - + -
+ - + - +
- + - - -
- - - - +
- - - - -
- - + - -
- - - - -
+ - + ]
580 : : !!***
581 : :
582 : : !!****f* ABINIT/cheb_poly
583 : : !! NAME
584 : : !! cheb_poly
585 : : !!
586 : : !! FUNCTION
587 : : !! Computes the value of the Chebyshev polynomial of degree n on the interval [a,b] at x
588 : : !!
589 : : !! INPUTS
590 : : !! x= input variable
591 : : !! n= degree
592 : : !! a= left bound of the interval
593 : : !! b= right bound of the interval
594 : : !!
595 : : !! OUTPUT
596 : : !! y= Tn(x)
597 : : !!
598 : : !! NOTES
599 : : !!
600 : : !! SOURCE
601 : :
602 : 288 : function cheb_poly(x, n, a, b) result(y)
603 : :
604 : : integer, intent(in) :: n
605 : : integer :: i
606 : : real(dp), intent(in) :: x, a, b
607 : : real(dp) :: y, xred, temp
608 : : real(dp) :: yim1
609 : :
610 : : ! *************************************************************************
611 : :
612 : 288 : xred = (x-(a+b)/2)/(b-a)*2
613 : 288 : y = xred
614 : 288 : yim1 = one
615 [ + + ]: 1728 : do i=2, n
616 : 1440 : temp = y
617 : 1440 : y = 2*xred*y - yim1
618 : 1728 : yim1 = temp
619 : : end do
620 : :
621 : 288 : end function cheb_poly
622 : : !!***
623 : :
624 : : !!****f* ABINIT/cheb_oracle
625 : : !! NAME
626 : : !! cheb_oracle
627 : : !!
628 : : !! FUNCTION
629 : : !! Returns the number of necessary iterations to decrease residual by at least tol
630 : : !! Here as in the rest of the code, the convention is that residuals are squared (||Ax-lx||^2)
631 : : !!
632 : : !! INPUTS
633 : : !! x= input variable
634 : : !! a= left bound of the interval
635 : : !! b= right bound of the interval
636 : : !! tol= needed precision
637 : : !! nmax= max number of iterations
638 : : !!
639 : : !! OUTPUT
640 : : !! n= number of iterations needed to decrease residual by tol
641 : : !!
642 : : !! NOTES
643 : : !!
644 : : !! SOURCE
645 : :
646 : : function cheb_oracle(x, a, b, tol, nmax) result(n)
647 : :
648 : : real(dp) :: tol
649 : :
650 : : integer :: nmax
651 : : integer :: n, i
652 : : real(dp), intent(in) :: x, a, b
653 : : real(dp) :: y, xred, temp
654 : : real(dp) :: yim1
655 : :
656 : : ! *************************************************************************
657 : :
658 : : xred = (x-(a+b)/2)/(b-a)*2
659 : : y = xred
660 : : yim1 = one
661 : :
662 : : n = nmax
663 : : if(1/(y**2) < tol) then
664 : : n = 1
665 : : else
666 : : do i=2, nmax-1
667 : : temp = y
668 : : y = 2*xred*y - yim1
669 : : yim1 = temp
670 : : if(1/(y**2) < tol) then
671 : : n = i
672 : : exit
673 : : end if
674 : : end do
675 : : end if
676 : :
677 : : end function cheb_oracle
678 : : !!***
679 : :
680 : : end module m_chebfi
681 : : !!***
|
