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