Branch data Line data Source code
1 : : !!****m* ABINIT/m_exc_itdiago
2 : : !! NAME
3 : : !! m_exc_itdiago
4 : : !!
5 : : !! FUNCTION
6 : : !! Iterative diagonalization of the BSE Hamiltonian with band-by-band conjugate gradient method
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (MG)
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_exc_itdiago
23 : :
24 : : use defs_basis
25 : : use m_bs_defs
26 : : use m_errors
27 : : use m_abicore
28 : : use m_linalg_interfaces
29 : : use m_hdr
30 : : USE_MPI
31 : : use m_xmpi
32 : :
33 : : use m_io_tools, only : open_file
34 : : use m_time, only : cwtime
35 : : use m_numeric_tools, only : stats_t, stats_eval
36 : : use m_hide_lapack, only : xhpev !xheev,
37 : : use m_bse_io, only : exc_read_rcblock
38 : :
39 : : implicit none
40 : :
41 : : private
42 : :
43 : : #ifdef HAVE_MPI1
44 : : include 'mpif.h'
45 : : #endif
46 : :
47 : : public :: exc_iterative_diago ! Calculates eigenvalues and eigenvectors of the Resonant BSE Hamiltonian
48 : : !!***
49 : :
50 : : CONTAINS !===============================================================
51 : : !!***
52 : :
53 : : !!****f* m_exc_itdiago/exc_iterative_diago
54 : : !!
55 : : !! NAME
56 : : !! exc_iterative_diago
57 : : !!
58 : : !! FUNCTION
59 : : !! Calculates eigenvalues and eigenvectors of the Hermitian excitonic Hamiltonian (coupling is neglected).
60 : : !!
61 : : !! INPUTS
62 : : !! Bsp
63 : : !! %nreh=Rank of the resonant block of the Hamiltonian.
64 : : !! %nstates=Number of eigenstates required.
65 : : !! %nline=Max number of line minimizations.
66 : : !! %tolwfr=Tolerance on the residuals.
67 : : !! %nbdbuf
68 : : !! %nstep
69 : : !! BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code.
70 : : !! %exh=Name of the file storing the excitonin resonant part.
71 : : !! %out_eig_out=Name of the file where final results are store.
72 : : !! %in_eig=Name of the file used to initialize the calculation.
73 : : !! comm=MPI communicator.
74 : : !!
75 : : !! OUTPUT
76 : : !! Eigenvalues and eigenvectors are written on file %out_eig
77 : : !!
78 : : !! NOTES
79 : : !! Concernig the approach followed to parallelize this routine: the most important
80 : : !! bottleneck is represent by the storage of the excitonic Hamiltonian since a
81 : : !! large number of k-points is needed to obtain converged exciton energies.
82 : : !! The number of eigenstates is usually much smaller than the rank of the full matrix,
83 : : !! this is especially true if we are only interested in the binding energy of the
84 : : !! exciton or in the excitonic states close to the single-particle gap.
85 : : !! Therefore good scaling and good performance should be obtained by distributing
86 : : !! the row of the excitonic Hamiltonian among the nodes while the required
87 : : !! eigenvectors are duplicated on each node. The memory needed to stores the eigenvalues
88 : : !! scales like nreh*nstates where nreh is the rank of the Hamiltonian and this might
89 : : !! render the calculation unfeasible when nstates is large.
90 : : !! On the other hand, having the complex set of trial eigenvectors on each node permits to parallelize
91 : : !! tasks such as the application of the Hamiltonian as well as the orthogonalization or the sub-space rotation
92 : : !! the later two algorithms represent the most CPU demanding part in standard KS calculations
93 : : !! as they scale with the third power of the number of atoms.
94 : : !! The conjugate direction and the gradient as well as Hphi are not distributed as the line minimization
95 : : !! requires the evaluation of <cg_dir_|H_exc|cg_dir>.
96 : : !! Note that this routine has been written having in mind an homogeneous network of machines.
97 : : !! A network made of different CPU will lead to unpredictable results as each node has
98 : : !! to check for the converge of the calculation.
99 : : !!
100 : : !! SOURCE
101 : :
102 : 1 : subroutine exc_iterative_diago(BSp,BS_files,Hdr_bse,prtvol,comm)
103 : :
104 : : !Arguments ------------------------------------
105 : : !scalars
106 : : integer,intent(in) :: comm,prtvol
107 : : type(excparam),intent(in) :: BSp
108 : : type(excfiles),intent(in) :: BS_files
109 : : type(Hdr_type),intent(in) :: Hdr_bse
110 : :
111 : : !Local variables ------------------------------
112 : : !scalars
113 : : integer,parameter :: STRICT=2,MEDIUM=1,WORST=0
114 : : integer,parameter :: master=0
115 : : integer(i8b) :: bsize_hmat,bsize_phi_block
116 : : integer :: hexc_size,nstates,nline,nbdbuf,nstep
117 : : integer :: cg_nsteps,my_t1,my_t2,my_nt
118 : : integer :: ii,jj,state,line,max_nline
119 : : integer :: cg_step,nbdbuf_,nsppol,ierr,nproc,my_rank
120 : : real(dp) :: exc_gap,exc_maxene,norm,etrial,etrial_old,deltae,tolwfr
121 : : ! real(dp) :: deold
122 : : real(dp) :: dhd,dhc,den,fac,poly,xx,root,swap,tan2th,diff,tolwfr_
123 : : !complex(dp) :: cg_gamma,dotgg,old_dotgg
124 : : real(dp) :: cg_gamma,dotgg,old_dotgg
125 : : real(dp) :: max_resid,costh,sinth
126 : : complex(dp) :: zz,kprc
127 : : logical,parameter :: DEBUGME=.False.
128 : : logical :: use_mpio,is_resonant,diago_is_real
129 : : character(len=500) :: msg
130 : : character(len=fnlen) :: hexc_fname,ihexc_fname,oeig_fname
131 : : type(stats_t) :: stats
132 : : !arrays
133 : 3 : integer :: nline_for(Bsp%nstates),convergence_of(Bsp%nstates)
134 : 3 : real(dp) :: resid(Bsp%nstates),exc_energy(Bsp%nstates),rbuf2(2)
135 : : ! real(dp),allocatable :: gsc(:,:),cg(:,:)
136 : : !complex,allocatable :: hexc(:,:)
137 : 1 : complex(dp),allocatable :: hexc(:,:),hji(:),vec_tmp(:)
138 : 1 : complex(dp),contiguous, pointer :: my_phi(:)
139 : 1 : real(dp),allocatable :: hexc_diagonal(:)
140 : 1 : complex(dp),target,allocatable :: phi_block(:,:)
141 : 1 : complex(dp),allocatable :: hphi(:) !,buffer_dpc(:)
142 : 1 : complex(dp),allocatable :: cg_dir(:),grad_dir(:),prc_dir(:)
143 : 1 : complex(dp),allocatable :: old_cg_dir(:)
144 : : !************************************************************************
145 : :
146 : : DBG_ENTER("COLL")
147 : :
148 : 1 : if (Bsp%use_coupling>0) then
149 : 0 : ABI_ERROR("CG Method does not support coupling")
150 : : end if
151 : :
152 : 1 : nsppol = Hdr_bse%nsppol
153 [ - + ]: 1 : if (Hdr_bse%nsppol == 2) then
154 : 0 : ABI_WARNING("nsppol==2 with cg method is still under development")
155 : : end if
156 : :
157 : 1 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
158 : :
159 : 1 : use_mpio=.FALSE.
160 : : #ifdef HAVE_MPI_IO
161 : : use_mpio = (nproc > 1)
162 : : #endif
163 : : use_mpio=.FALSE.
164 : : !use_mpio = .TRUE.
165 : :
166 [ + + ]: 2 : hexc_size = SUM(Bsp%nreh)
167 : 1 : nstates= Bsp%nstates
168 : 1 : nline = Bsp%nline
169 : 1 : nbdbuf = Bsp%nbdbuf
170 : 1 : nstep = Bsp%niter
171 : 1 : tolwfr = Bsp%cg_tolwfr
172 : :
173 : 1 : write(msg,'(a,i0)')' Iterative diagonalization of the resonant excitonic Hamiltonian, Matrix size= ',hexc_size
174 : 1 : call wrtout(std_out,msg,"COLL")
175 : 1 : call wrtout(ab_out,msg,"COLL")
176 : :
177 [ - + ]: 1 : ABI_CHECK(hexc_size>=nproc,"hexc_size<nproc!")
178 [ - + ]: 1 : ABI_CHECK(nstates <= hexc_size,"nstates cannot be greater that hexc size!")
179 : :
180 : : ! Divide the columns of the Hamiltonian among the nodes.
181 : : call xmpi_split_work(hexc_size,comm,my_t1,my_t2)
182 : :
183 : 1 : my_nt = my_t2-my_t1+1
184 : 1 : write(msg,'(a,i0,a)')" Will handle ",my_nt," columns of the excitonic Hamiltonian. "
185 : 1 : call wrtout(std_out,msg,"PERS")
186 : :
187 : 1 : tolwfr_ = tolwfr
188 [ - + ]: 1 : if (tolwfr < 10**(-30)) then
189 : 0 : tolwfr_ = tol12
190 : 0 : write(msg,'(2(a,es12.4))')" Input tolwfr= ",tolwfr," Using tolwfr= ",tolwfr_
191 : 0 : ABI_WARNING(msg)
192 : : end if
193 : :
194 : 1 : cg_nsteps = nstep
195 [ - + ]: 1 : if (cg_nsteps<=0) then
196 : 0 : cg_nsteps = 30
197 : 0 : write(msg,'(2(a,es12.4))')" Input nstep= ",nstep," Using cg_nsteps= ",cg_nsteps
198 : 0 : ABI_WARNING(msg)
199 : : end if
200 : :
201 : 1 : nbdbuf_ = nbdbuf
202 [ - + ]: 1 : if (nbdbuf<=0) then
203 : 0 : nbdbuf_ = 4
204 : 0 : write(msg,'(2(a,i0))')" Input nbdbuf= ",nbdbuf," Using nbdbuf= ",nbdbuf_
205 : 0 : ABI_WARNING(msg)
206 : : end if
207 : :
208 : : write(msg,"(4(a,i0,a),a,es12.4)")&
209 : 1 : & " cg_nsteps: ",cg_nsteps,ch10,&
210 : 1 : & " nstates: ",nstates,ch10,&
211 : 1 : & " nline: ",nline,ch10,&
212 : 1 : & " nbdbuf: ",nbdbuf_,ch10,&
213 : 2 : & " tolwfr: ",tolwfr_
214 : 1 : call wrtout(std_out,msg,"COLL")
215 : 1 : call wrtout(ab_out,msg,"COLL")
216 : :
217 : 1 : bsize_hmat = 2*dp*hexc_size*my_nt
218 : 1 : write(msg,'(a,f8.1,a)')' Allocating excitonic Hamiltonian. Memory requested: ',bsize_hmat*b2Mb,' Mb.'
219 : 1 : call wrtout(std_out,msg,"COLL")
220 : :
221 [ - + + - : 3 : ABI_MALLOC(hexc_diagonal,(my_t1:my_t2))
+ - + - -
+ - + ]
222 [ - + + - : 4 : ABI_MALLOC_OR_DIE(hexc,(hexc_size,my_t1:my_t2), ierr)
- + + - +
- + - + -
+ - - + +
- - + ]
223 : : !
224 : : ! Read and construct full excitonic Hamiltonian using Hermiticity.
225 [ - + ]: 1 : if (BS_files%in_hreso /= BSE_NOFILE) then
226 : 0 : hexc_fname = BS_files%in_hreso
227 : : else
228 : 1 : hexc_fname = BS_files%out_hreso
229 : : end if
230 : : !
231 : : ! Read the resonant block from file.
232 : 1 : is_resonant=.TRUE.
233 : 1 : diago_is_real=(.not.BSp%have_complex_ene)
234 : 1 : call exc_read_rcblock(hexc_fname,Bsp,is_resonant,diago_is_real,nsppol,Bsp%nreh,hexc_size,my_t1,my_t2,hexc,use_mpio,comm)
235 : : !
236 : : ! Save diagonal part for preconditioning.
237 [ + + ]: 385 : do jj=my_t1,my_t2
238 : 385 : hexc_diagonal(jj) = REAL(hexc(jj,jj),kind=dp)
239 : : end do
240 : : !
241 : : ! === Initialisation of the excitonic wavefunctions ===
242 : : ! Two cases are possible.
243 : : ! 1) Fill trial eigenvectors with random numbers
244 : : ! One needs to initialize wfs in such a way to avoid symmetry traps,
245 : : ! and to avoid linear dependencies between wavefunctions
246 : : ! 2) Read eigenstates generated by a previous calculation.
247 : :
248 : 1 : bsize_phi_block = 2*sp*my_nt*nstates
249 : 1 : write(msg,'(a,f8.1,a)')' Allocating BSE eigenvectors. Memory requested: ',bsize_phi_block*b2Mb,' Mb.'
250 : 1 : call wrtout(std_out,msg,"COLL",do_flush=.True.)
251 : :
252 [ - + + - : 4 : ABI_MALLOC_OR_DIE(phi_block,(my_t1:my_t2,nstates), ierr)
- + + - +
- + - + -
+ - - + +
- - + ]
253 : :
254 : 1 : ihexc_fname = ""
255 [ - + ]: 1 : if (BS_files%in_eig /= BSE_NOFILE) ihexc_fname = BS_files%in_eig
256 : :
257 : 1 : call exc_init_phi_block(ihexc_fname,use_mpio,comm)
258 : : !
259 : : ! =========================
260 : : ! === Orthogonalization ===
261 : : ! =========================
262 : 1 : call exc_cholesky_ortho()
263 : : call exc_check_phi_block("First cholesky ortho")
264 : :
265 : : ! * Sub-space rotation.
266 : 2 : call exc_subspace_rotation()
267 : : call exc_check_phi_block("First subspace_rotation")
268 : : !
269 : : ! ===========================
270 : : ! ==== Conjugate gradient ===
271 : : ! ===========================
272 [ + - + - : 3 : ABI_MALLOC(hphi,(hexc_size))
+ - - + -
+ ]
273 [ - + + - : 2 : ABI_MALLOC(cg_dir,(hexc_size))
- + - + ]
274 [ - + + - : 2 : ABI_MALLOC(old_cg_dir,(hexc_size))
- + - + ]
275 [ - + + - : 2 : ABI_MALLOC(grad_dir,(hexc_size))
- + - + ]
276 [ - + + - : 2 : ABI_MALLOC(prc_dir,(hexc_size))
- + - + ]
277 : :
278 : : max_nline=nline
279 [ + + + + : 31 : resid(:)=HUGE(one); nline_for(1:nstates)=max_nline; convergence_of(1:nstates)=WORST
+ + ]
280 : :
281 [ + + ]: 21 : do cg_step=1,cg_nsteps
282 [ + + ]: 220 : do state=1,nstates
283 : :
284 [ - + ]: 200 : if (prtvol>=10) then ! Tell us what is going on:
285 : 0 : write(msg,'(a,i6,2x,a,i3,a)')' --- exc_iterative_diago is called for state ',state,'for',nline_for(state),' lines'
286 : 0 : call wrtout(std_out,msg,'PERS')
287 : : end if
288 : :
289 : : ! Extraction of the vector that is iteratively updated.
290 : 200 : my_phi => phi_block(my_t1:my_t2,state)
291 : :
292 [ + + ]: 933 : do line=1,nline_for(state)
293 : : ! Compute etrial=<phi|H|phi> and the residual [H-etrial]|phi>.
294 [ + + ]: 283360 : hphi = czero
295 : : #ifdef FC_NVHPC
296 : : !Buggy NVHPC compiler
297 : : do jj=1,my_t2-my_t1+1;do ii=1,hexc_size
298 : : hphi(ii)=hphi(ii)+hexc(ii,jj)*my_phi(jj)
299 : : enddo ; enddo
300 : : #else
301 [ + - - + : 326148832 : hphi = MATMUL(hexc, my_phi)
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + + + +
+ + + - +
- + - ]
302 : : #endif
303 : 736 : call xmpi_sum(hphi,comm,ierr)
304 : :
305 [ + + ]: 283360 : etrial = DOT_PRODUCT(my_phi, hphi(my_t1:my_t2))
306 : 736 : call xmpi_sum(etrial,comm,ierr)
307 : 736 : exc_energy(state) = etrial
308 : :
309 : : ! Compute residual (squared) norm.
310 [ + + ]: 283360 : grad_dir(my_t1:my_t2) = hphi(my_t1:my_t2) - etrial*my_phi
311 [ + + ]: 283360 : resid(state) = DOT_PRODUCT(grad_dir(my_t1:my_t2), grad_dir(my_t1:my_t2))
312 : 736 : call xmpi_sum(resid(state),comm,ierr)
313 : 1472 : convergence_of(state) = convergence_degree(resid(state))
314 : : !
315 : : ! Check that etrial is decreasing on succeeding lines:
316 [ + + - + ]: 736 : if (line>1 .and. (etrial > etrial_old+tol12)) then
317 : : write(msg,'(a,i8,a,1p,e14.6,a1,3x,a,1p,e14.6,a1)')&
318 : 0 : & 'New trial exc_energy at line ',line,' = ',etrial,ch10,&
319 : 0 : & 'is higher than former:',etrial_old,ch10
320 : 0 : ABI_WARNING(msg)
321 : : end if
322 : 736 : etrial_old = etrial
323 : : !
324 : : ! If residual sufficiently small stop line minimization.
325 [ + + ]: 736 : if (convergence_of(state)==STRICT) then
326 [ - + ]: 23 : if (prtvol>=10) then
327 : : write(msg,'(a,i4,a,i2,a,es12.4)')&
328 : 0 : & ' exc_iterative_diago: state ',state,' converged after ',line,&
329 : 0 : & ' line minimizations : resid =',resid(state)
330 : 0 : call wrtout(std_out,msg,'PERS')
331 : : end if
332 : : EXIT !line
333 : : end if
334 : :
335 : : ! === PROJECT THE STEEPEST DESCENT DIRECTION OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ===
336 : : ! The following projection over the subspace orthogonal to occupied bands
337 : : ! is optional. It is a bit more accurate, but doubles the number of N^3 ops.
338 : : ! It is done only if ortalg>=0.
339 : :
340 : : ! Project the steepest descent direction: direc(2,npw)=<G|H|Cnk> - \sum_{(i<=n)} <G|H|Cik> , normalized.
341 : :
342 : : ! Grad_dir is already orthogonal to this band
343 [ - + + - : 2139 : ABI_MALLOC(hji,(nstates))
+ - + - -
+ - + -
+ ]
344 [ + + ]: 7843 : hji=czero
345 : :
346 : : ! MG TODO Don't know why here we sum over i=<=n!!!!!!!!
347 [ + + ]: 7843 : do jj=1,nstates
348 [ + + + + ]: 2471971 : if (jj/=state) hji(jj) = DOT_PRODUCT(phi_block(:,jj), hphi(my_t1:my_t2) )
349 : : end do
350 : 713 : call xmpi_sum(hji,comm,ierr)
351 : :
352 [ + + ]: 7843 : do jj=1,nstates
353 [ + + + + ]: 2471971 : if (jj/=state) grad_dir(my_t1:my_t2) = grad_dir(my_t1:my_t2) - hji(jj)*phi_block(:,jj)
354 : : end do
355 [ - + ]: 713 : ABI_FREE(hji)
356 : : !
357 : : ! === PRECONDITION THE STEEPEST DESCENT DIRECTION ===
358 [ + + ]: 274505 : den = DOT_PRODUCT(grad_dir(my_t1:my_t2), hexc_diagonal(my_t1:my_t2)*grad_dir(my_t1:my_t2) )
359 : 713 : call xmpi_sum(den,comm,ierr)
360 : :
361 [ + + ]: 274505 : do ii=my_t1,my_t2
362 : : ! Teter polynomial ratio, modified according to Kresse, Furthmuller, PRB 54, 11169 (1996) [[cite:Kresse1996]]
363 : 273792 : xx = hexc_diagonal(ii)/den
364 : 273792 : poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
365 : 273792 : fac=poly/(poly+16._dp*xx**4)
366 : 273792 : kprc = fac*four/(three*den)
367 : 274505 : prc_dir(ii) = kprc * grad_dir(ii)
368 : : end do
369 : : !
370 : : ! * PROJECT THE PRECOND. STEEPEST DESCENT DIRECTION OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS.
371 [ - + + - : 2139 : ABI_MALLOC(hji,(nstates))
+ - + - -
+ - + ]
372 [ + + ]: 7843 : hji=czero
373 [ + + ]: 7843 : do jj=1,nstates
374 [ + + ]: 2745763 : hji(jj) = DOT_PRODUCT(phi_block(:,jj), prc_dir(my_t1:my_t2) )
375 : : end do
376 : 713 : call xmpi_sum(hji,comm,ierr)
377 : :
378 [ + + ]: 7843 : do jj=1,nstates
379 [ + + ]: 2745763 : prc_dir(my_t1:my_t2) = prc_dir(my_t1:my_t2) - hji(jj)*phi_block(:,jj)
380 : : end do
381 [ - + ]: 713 : ABI_FREE(hji)
382 : : !
383 : : ! === COMPUTE THE CONJUGATE-GRADIENT ===
384 [ + + ]: 274505 : dotgg = DOT_PRODUCT(prc_dir(my_t1:my_t2),grad_dir(my_t1:my_t2))
385 : 713 : call xmpi_sum(dotgg,comm,ierr)
386 : :
387 [ + + ]: 713 : if (line==1) then ! At first iteration, cg_gamma is set to zero
388 : 179 : cg_gamma=zero
389 : 179 : old_dotgg=dotgg
390 [ - + + + ]: 69094 : cg_dir = prc_dir
391 [ - + + + ]: 69094 : old_cg_dir = cg_dir
392 : : else
393 : 534 : cg_gamma=dotgg/old_dotgg
394 : 534 : old_dotgg=dotgg
395 : : !write(std_out,*)"cg_gamma= ",cg_gamma
396 : : !cg_dir = prc_dir + cg_gamma*cg_dir
397 [ - + + + ]: 206124 : cg_dir = prc_dir + cg_gamma*old_cg_dir !TODO check this, anyhow it is much faster.
398 [ - + + + ]: 206124 : old_cg_dir =cg_dir ! old_cg_dir is used to store the previsou CG direction, cg_dir will be orthonormalized to the band
399 : : end if
400 : : !
401 : : ! === PROJECTION OF THE CONJUGATED GRADIENT ===
402 [ + + ]: 274505 : zz = DOT_PRODUCT(my_phi, cg_dir(my_t1:my_t2))
403 : 713 : call xmpi_sum(zz,comm,ierr)
404 [ + + ]: 274505 : cg_dir(my_t1:my_t2) = cg_dir(my_t1:my_t2) -zz*my_phi(:)
405 : :
406 [ + + ]: 274505 : norm = DOT_PRODUCT(cg_dir(my_t1:my_t2), cg_dir(my_t1:my_t2) )
407 : 713 : call xmpi_sum(norm,comm,ierr)
408 : 713 : norm = SQRT(norm)
409 [ + + ]: 274505 : cg_dir = cg_dir/norm ! Have to normalize it.
410 : :
411 : : ! Line minimization of the Raileigh functional.
412 [ - + + - : 2139 : ABI_MALLOC(vec_tmp,(hexc_size))
+ - + - -
+ - + -
+ ]
413 [ + + ]: 274505 : vec_tmp=czero
414 : : #ifdef FC_NVHPC
415 : : !Buggy NVHPC compiler
416 : : do jj=my_t1,my_t2;do ii=1,hexc_size
417 : : vec_tmp(ii)=vec_tmp(ii)+hexc(ii,jj-my_t1+1)*cg_dir(jj)
418 : : enddo ; enddo
419 : : #else
420 [ + - - + : 210820553 : vec_tmp = MATMUL(hexc, cg_dir(my_t1:my_t2))
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + + + +
+ + + - +
- ]
421 : : #endif
422 : 713 : call xmpi_sum(vec_tmp,comm,ierr)
423 : :
424 : : !if (my_rank==master) then
425 : : ! write(777,*)"cg_step, state, line",cg_step, state, line
426 : : ! write(777,*)vec_tmp
427 : : !end if
428 : :
429 [ + + ]: 274505 : dhd = DOT_PRODUCT( cg_dir(my_t1:my_t2), vec_tmp(my_t1:my_t2)) ! is this always real?
430 [ + + ]: 274505 : dhc = REAL( DOT_PRODUCT( my_phi, vec_tmp(my_t1:my_t2) ))
431 [ - + ]: 713 : ABI_FREE(vec_tmp)
432 : :
433 [ + + ]: 2139 : rbuf2 = (/dhd,dhc/)
434 : 713 : call xmpi_sum(rbuf2,comm,ierr)
435 : 713 : dhd = rbuf2(1)
436 : 713 : dhc = rbuf2(2)
437 : :
438 : : !write(201*(my_rank+1),*)"cg_step, state, line dotgg dhd dhc",cg_step,state,line,dotgg,dhd,dhc
439 : :
440 : : ! Taken from cgwf
441 : : ! Compute tan(2 theta),sin(theta) and cos(theta)
442 : 713 : tan2th=2.0_dp*dhc/(etrial-dhd)
443 : :
444 [ + + ]: 713 : if (abs(tan2th)<1.d-05) then
445 : 67 : costh=1.0_dp-0.125_dp*tan2th**2
446 : 67 : sinth=0.5_dp*tan2th*(1.0_dp-0.375_dp*tan2th**2)
447 : : ! Check that result is above machine precision
448 : : ! FIXME This part is not safe on clusters made of different machines or different compiation options.
449 [ - + ]: 67 : if (abs(sinth)<epsilon(0._dp)) then
450 : 0 : write(msg, '(a,es16.4)' ) ' exc_iterative_diago: converged with tan2th= ',tan2th
451 : 0 : call wrtout(std_out,msg,'PERS')
452 : 0 : EXIT !Exit from the loop on line
453 : : end if
454 : :
455 : : else
456 : 646 : root =sqrt(1.0_dp+tan2th**2)
457 : 646 : costh=sqrt(0.5_dp+0.5_dp/root)
458 : 646 : sinth=sign(sqrt(0.5_dp-0.5_dp/root),tan2th)
459 : : end if
460 : : !
461 : : ! Check for lower of two possible roots (same sign as curvature at theta where slope is zero)
462 : 713 : diff=(etrial-dhd)
463 [ + + ]: 713 : if (diff>zero) then ! Swap c and d if value of diff is positive
464 : 10 : swap=costh
465 : 10 : costh=-sinth
466 : 10 : sinth=swap
467 [ - + ]: 10 : if (prtvol<0 .or. prtvol>=10) then
468 : 0 : write(msg,'(a,i4,es16.6)')' exc_iterative_diago: swap roots, line,diff= ',line,diff
469 : 0 : call wrtout(std_out,msg,'PERS')
470 : : end if
471 : : end if
472 : : !
473 : : ! === GENERATE NEW |wf>, H|wf> =============
474 [ + + ]: 274505 : my_phi = costh*my_phi + sinth*cg_dir(my_t1:my_t2)
475 : : !write(100*(my_rank+1),*)"cg_step state, line costh sinth etrial",cg_step,state,line,costh,sinth,etrial
476 : : !end taken from cgwf
477 : :
478 : : !norm = SQRT( DOT_PRODUCT(my_phi,my_phi) )
479 : : !my_phi = my_phi /norm
480 : : !write(std_out,*)norm
481 : : !write(std_out,*)DOT_PRODUCT(hphi,my_phi),cos(theta_min)
482 : :
483 : : ! ======================================================================
484 : : ! =========== CHECK CONVERGENCE AGAINST TRIAL ENERGY ===================
485 : : ! ======================================================================
486 : : ! Compute delta(E)
487 : : !deltae=chc*(costh**2-1._dp)+dhd*sinth**2+2._dp*costh*sinth*dhc
488 : 8089 : deltae=etrial*(costh**2-1._dp)+dhd*sinth**2+2._dp*costh*sinth*dhc
489 : :
490 : : ! Check convergence and eventually exit
491 : : ! if (line==1) then
492 : : ! deold=deltae
493 : : ! else if (abs(deltae)<0.005_dp*abs(deold) .and. line/=nline_for(state))then
494 : : ! if (prtvol>=10)then
495 : : ! write(msg, '(a,i4,1x,a,1p,e12.4,a,e12.4,a)' ) &
496 : : !& ' cgwf: line',line,&
497 : : !& ' deltae=',deltae,' < 0.005*',deold,' =>skip lines'
498 : : ! call wrtout(std_out,msg,'PERS')
499 : : ! end if
500 : : ! exc_energy(state) = exc_energy(state) + deltae
501 : : ! EXIT
502 : : ! end if
503 : : end do ! LOOP FOR A GIVEN BAND. Note that there are three "exit" instructions inside
504 : : ! Modify nline_for(state) according to converge degree.
505 : : !if (convergence_of(state) == STRICT) nline_for(state) = MAX(max_nline-2,2)
506 : : !if (convergence_of(state) == MEDIUM) nline_for(state) = MAX(max_nline-1,2)
507 : : !if (convergence_of(state) == WORST ) nline_for(state) = max_nline
508 : : end do !state
509 : :
510 [ - + ]: 20 : if (prtvol>2) then
511 [ # # ]: 0 : do ii=0,(nstates-1)/8
512 : 0 : write(msg,'(a,8es10.2)')' res:',(resid(state),state=1+ii*8,MIN(nstates,8+ii*8))
513 : 0 : call wrtout(std_out,msg,'COLL')
514 : : end do
515 [ # # ]: 0 : do ii=0,(nstates-1)/8
516 : 0 : write(msg,'(a,8es10.2)')' ene:',(exc_energy(state),state=1+ii*8,MIN(nstates,8+ii*8))
517 : 0 : call wrtout(std_out,msg,'COLL')
518 : : end do
519 : : end if
520 : :
521 : 20 : write(msg,'(a,i0)')"After cg_step: ",cg_step
522 : : call exc_check_phi_block(msg)
523 : :
524 : : ! Find largest residual over bands and Print residuals
525 [ - + + - : 160 : max_resid=MAXVAL( resid(:MAX(1,nstates-nbdbuf_)) )
+ + ]
526 : :
527 [ - + ]: 20 : if (max_resid < tolwfr_) then
528 : : write(msg,'(a,i0,2(a,es10.2),a,i0,a)')&
529 : 0 : & " After ",cg_step," iterations, max_resid= ",max_resid," < tolwfr= ",tolwfr_," ( Excluding nbdbuf= ",nbdbuf_,")"
530 : 0 : call wrtout(std_out,msg,'COLL')
531 : 0 : EXIT ! cg_step
532 : : end if
533 : :
534 : : if (cg_step==1.or.MOD(cg_step,1)==0) then
535 : 20 : call wrtout(std_out," Subspace rotation + exc_cholesky_ortho ","COLL")
536 : :
537 : 20 : call exc_subspace_rotation()
538 : 20 : call exc_cholesky_ortho()
539 : :
540 : : !mcg=hexc_size; mgsc=hexc_size; useoverlap=0
541 : : !allocate(cg(2,mcg),gsc(2,mgsc*useoverlap))
542 : : !do ii=1,nstates
543 : : ! cg(1,:) = REAL (phi_block(:,ii))
544 : : ! cg(2,:) = AIMAG(phi_block(:,ii))
545 : : ! call fxphas(cg,gsc,0,0,1,mcg,mgsc,MPI_enreg_seq,1,hexc_size,useoverlap)
546 : : ! phi_block(:,ii)=CMPLX(cg(1,:),cg(2,:))
547 : : !end do
548 : : !deallocate(cg,gsc)
549 : :
550 : : end if
551 : :
552 : : !XG20141126 : At present, keep this fake test. Seems that there is some bug on shiva
553 : : !(compiler ?) so that the result is erroneous without it. Even suppressing the
554 : : !"stop" instruction triggers the bug ...
555 : 1 : if(cg_step==cg_nsteps+1)then
556 : : write(std_out,*)' One should not be here : cg_step==cg_nsteps+1 '
557 : : stop
558 : : endif
559 : :
560 : : end do !cg_step
561 : :
562 : : ! Release some memory before entering RMM-DIIS
563 [ - + ]: 1 : ABI_FREE(hphi)
564 [ - + ]: 1 : ABI_FREE(cg_dir)
565 [ - + ]: 1 : ABI_FREE(old_cg_dir)
566 [ - + ]: 1 : ABI_FREE(grad_dir)
567 [ - + ]: 1 : ABI_FREE(prc_dir)
568 : :
569 [ + + ]: 3 : do ii=0,(nstates-1)/8
570 : 2 : write(msg,'(a,8es10.2)')' res:',(resid(state),state=1+ii*8,min(nstates,8+ii*8))
571 : 3 : call wrtout(std_out,msg,'COLL')
572 : : end do
573 [ + + ]: 3 : do ii=0,(nstates-1)/8
574 : 2 : write(msg,'(a,8es10.2)')' ene:',(exc_energy(state),state=1+ii*8,min(nstates,8+ii*8))
575 : 3 : call wrtout(std_out,msg,'COLL')
576 : : end do
577 : :
578 : 1 : stats = stats_eval(resid)
579 : :
580 : : write(msg,"(4(a,es10.2))")&
581 : 1 : & ". Residuals: min value: ",stats%min,", Max value: ",stats%max,", mean: ",stats%mean,", stdev: ",stats%stdev
582 : 1 : call wrtout(std_out,msg)
583 : 1 : call wrtout(ab_out,msg)
584 : :
585 [ + - ]: 1 : if (max_resid > tolwfr_) then
586 : 1 : write(msg,'(2a,i5,2a,2(a,es10.2),a,i3,a)')ch10,&
587 : 1 : & " WARNING: conjugate-gradient not converged after ",cg_step," iterations.",ch10,&
588 : 2 : & " max_resid= ",max_resid," > tolwfr= ",tolwfr_," ( Excluding nbdbuf= ",nbdbuf_,")"
589 : 1 : call wrtout(ab_out,msg,'COLL')
590 : 1 : call wrtout(std_out,msg,'COLL')
591 : : end if
592 : :
593 [ + - - + : 12 : exc_gap = MINVAL(exc_energy)
- - + + ]
594 [ + - - + : 12 : exc_maxene = MAXVAL(exc_energy)
- - + + ]
595 : :
596 : 1 : write(msg,'(a,2(a,f7.2,2a))')ch10,&
597 : 1 : & " First excitonic eigenvalue= ",exc_gap*Ha_eV, " [eV]. ",ch10,&
598 : 2 : & " Last excitonic eigenvalue= ",exc_maxene*Ha_eV," [eV]. ",ch10
599 : 1 : call wrtout(std_out,msg,"COLL")
600 : 1 : call wrtout(ab_out,msg,"COLL")
601 : :
602 : : call exc_check_phi_block("END OF CONJUGATE-GRADIENT")
603 : :
604 : 0 : ABI_FREE(hexc)
605 [ - + ]: 1 : ABI_FREE(hexc_diagonal)
606 : : !
607 : : ! =====================================
608 : : ! ==== Write final results on file ====
609 : : ! =====================================
610 : 1 : oeig_fname = BS_files%out_eig
611 [ - + ]: 1 : if (oeig_fname== BSE_NOFILE) then
612 : 0 : ABI_WARNING("oeig_fname was set to "//TRIM(BSE_NOFILE))
613 [ # # ]: 0 : oeig_fname = TRIM(BS_files%out_basename)//"_BSEIG"
614 [ # # ]: 0 : ABI_WARNING("using oeig_fname : "//TRIM(oeig_fname))
615 : : end if
616 : :
617 : 1 : call exc_write_phi_block(oeig_fname,use_mpio)
618 : :
619 [ - + ]: 1 : ABI_FREE(phi_block)
620 : :
621 [ - + - + : 3 : call xmpi_barrier(comm)
- + - + -
+ ]
622 : :
623 : : DBG_EXIT("COLL")
624 : :
625 : : CONTAINS !===========================================================
626 : : !!***
627 : :
628 : : !----------------------------------------------------------------------
629 : :
630 : : !!****f* m_exc_itdiago/exc_init_phi_block
631 : : !! NAME
632 : : !! exc_init_phi_block
633 : : !!
634 : : !! FUNCTION
635 : : !! Initialize the eigenstates either from file or fill them with random number
636 : : !! if restart file is not available
637 : : !!
638 : : !! INPUTS
639 : : !! ihexc_fname=
640 : : !! Name of the file from which the eigenvectors will be read.
641 : : !! Empty string to initialize trial eigenvectors with random numbers.
642 : : !!
643 : : !! SIDE EFFECTS
644 : : !! phi_block(my_t1:my_t2,nstates)=Contains the trial eigenstates.
645 : : !!
646 : : !! SOURCE
647 : :
648 : 1 : subroutine exc_init_phi_block(ihexc_fname,use_mpio,comm)
649 : :
650 : : !Arguments ------------------------------------
651 : : !scalars
652 : : integer,intent(in) :: comm
653 : : logical,intent(in) :: use_mpio
654 : : character(len=*),intent(in) :: ihexc_fname
655 : :
656 : : !Local variables ------------------------------
657 : : integer :: eig_unt,hexc_size_restart,ii,state,seed
658 : : integer :: fold1,fold2,foldim,foldre
659 : : real(dp) :: cputime,walltime,gflops
660 : : character(len=500) :: errmsg
661 : : !arrays
662 : 1 : complex(dp),allocatable :: buffer_dpc(:)
663 : : #ifdef HAVE_MPI_IO
664 : : integer:: amode,mpi_fh,mpi_err,old_type,etype,eig_type,my_nel,ierr,my_nrows
665 : : integer(XMPI_OFFSET_KIND) :: ehdr_offset,my_offset,my_offpad,fmarker
666 : 1 : integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecord(:)
667 : : integer :: array_of_sizes(2),array_of_subsizes(2),array_of_starts(2)
668 : : #endif
669 : : !************************************************************************
670 : :
671 [ + - ]: 1 : if (LEN_TRIM(ihexc_fname) == 0) then
672 : 1 : call wrtout(std_out," Initializing eigenvectors with random numbers","COLL")
673 : : !
674 : : ! Use random number generator. For portability, use only integer numbers
675 : : ! The series of couples (fold1,fold2) is periodic with a period of
676 : : ! 3x5x7x11x13x17x19x23x29x31, that is, larger than 2**32, the largest integer*4
677 : : ! fold1 is between 0 and 34, fold2 is between 0 and 114. As sums of five
678 : : ! uniform random variables, their distribution is close to a gaussian
679 : : ! the gaussian distributions are folded, in order to be back to a uniform distribution
680 : : ! foldre is between 0 and 20, foldim is between 0 and 18.
681 : : !
682 [ + + ]: 11 : do state=1,nstates
683 [ + + ]: 3851 : do ii=my_t1,my_t2
684 : 3840 : seed=ii+(state-1)*hexc_size ! Different seed for different transitions and bands
685 : 3840 : fold1 =mod(seed,3)+mod(seed,5)+mod(seed,7)+mod(seed,11)+mod(seed,13)
686 : 3840 : fold2 =mod(seed,17)+mod(seed,19)+mod(seed,23)+mod(seed,29)+mod(seed,31)
687 : 3840 : foldre=mod(fold1+fold2,21)
688 : 3840 : foldim=mod(3*fold1+2*fold2,19)
689 : :
690 : 3850 : phi_block(ii,state) = DCMPLX(foldre,foldim)
691 : : end do
692 : : end do
693 : :
694 : : else
695 : :
696 : 0 : call cwtime(cputime, walltime, gflops, "start")
697 : :
698 [ # # ]: 0 : if (.not.use_mpio) then
699 [ # # ]: 0 : call wrtout(std_out," Initializing eigenvectors from file: "//TRIM(ihexc_fname)//" using Fortran IO.","COLL")
700 : :
701 [ # # ]: 0 : if (open_file(ihexc_fname,msg,newunit=eig_unt,form='unformatted',status="old") /=0 ) then
702 : 0 : ABI_ERROR(msg)
703 : : end if
704 : :
705 [ # # ]: 0 : read(eig_unt, err=10, iomsg=errmsg) hexc_size_restart
706 [ # # ]: 0 : ABI_CHECK(hexc_size_restart==hexc_size,"hexc_size_restart /= hexc_size")
707 [ # # ]: 0 : read(eig_unt, err=10, iomsg=errmsg) !skip DCMPLX(exevl(1:hexc_size))
708 : :
709 [ # # # # : 0 : ABI_MALLOC(buffer_dpc,(hexc_size))
# # # # #
# # # #
# ]
710 [ # # ]: 0 : do ii=1,nstates
711 [ # # ]: 0 : read(eig_unt, err=10, iomsg=errmsg) buffer_dpc
712 [ # # ]: 0 : phi_block(my_t1:my_t2,ii) = buffer_dpc(my_t1:my_t2)
713 : : end do
714 [ # # ]: 0 : ABI_FREE(buffer_dpc)
715 : :
716 [ # # ]: 0 : close(eig_unt, err=10, iomsg=errmsg)
717 : : else
718 [ # # ]: 0 : call wrtout(std_out," Initializing eigenvectors from file: "//TRIM(ihexc_fname)//" using MPI-IO.","COLL")
719 : : #ifdef HAVE_MPI_IO
720 : : !
721 : : ! Open the file with MPI-IO
722 : 0 : amode=MPI_MODE_RDONLY
723 : :
724 : 0 : call MPI_FILE_OPEN(comm, ihexc_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
725 [ # # ]: 0 : msg = " MPI_IO error opening file: "//TRIM(ihexc_fname)
726 : 0 : ABI_CHECK_MPI(mpi_err,msg)
727 : :
728 : : ! Move the file pointer to skip the first two records.
729 : 0 : ehdr_offset=0
730 : 0 : call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err)
731 : 0 : write(std_out,*)"fmarker first record ",fmarker
732 : 0 : call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err)
733 : 0 : write(std_out,*)"fmarker first record ",fmarker
734 : : !%call hdr_mpio_skip(mpi_fh,fform,ehdr_offset)
735 : : !%ehdr_offset = 4*xmpio_bsize_frm + xmpio_bsize_int + nstates*xmpio_bsize_dpc
736 : :
737 : 0 : etype=MPI_BYTE; old_type=MPI_DOUBLE_COMPLEX
738 : :
739 : 0 : my_nrows=my_t2-my_t1+1; old_type=MPI_DOUBLE_COMPLEX
740 [ # # ]: 0 : array_of_sizes = (/hexc_size,nstates/)
741 [ # # ]: 0 : array_of_subsizes = (/my_nrows,nstates/)
742 [ # # ]: 0 : array_of_starts = (/my_t1,1/)
743 : 0 : call xmpio_create_fsubarray_2D(array_of_sizes,array_of_subsizes,array_of_starts,old_type,eig_type,my_offpad,mpi_err)
744 : 0 : ABI_CHECK_MPI(mpi_err,"fsubarray_2D")
745 : : !
746 : : ! Each node uses a different offset to skip the header and the blocks written by the other CPUs.
747 : 0 : my_offset = ehdr_offset + my_offpad
748 : :
749 : 0 : call MPI_FILE_SET_VIEW(mpi_fh, my_offset, etype, eig_type, 'native', MPI_INFO_NULL, mpi_err)
750 : 0 : ABI_CHECK_MPI(mpi_err,"SET_VIEW")
751 : :
752 : 0 : call MPI_TYPE_FREE(eig_type,mpi_err)
753 : 0 : ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE")
754 : :
755 : 0 : my_nel = my_nrows*nstates
756 : 0 : call MPI_FILE_READ_ALL(mpi_fh, phi_block, my_nel, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err)
757 : 0 : ABI_CHECK_MPI(mpi_err,"FILE_READ")
758 : :
759 : : ! It seems that personal calls make the code stuck
760 : : ! check the fortran markers.
761 [ # # # # : 0 : ABI_MALLOC(bsize_frecord,(nstates))
# # # # #
# # # ]
762 [ # # ]: 0 : bsize_frecord = hexc_size * xmpi_bsize_dpc
763 : : ! ehdr_offset points to the end of the header.
764 : 0 : call xmpio_check_frmarkers(mpi_fh,ehdr_offset,xmpio_collective,nstates,bsize_frecord,ierr)
765 [ # # ]: 0 : ABI_CHECK(ierr==0,"Error in Fortran markers")
766 [ # # ]: 0 : ABI_FREE(bsize_frecord)
767 : : !
768 : : ! Close the file.
769 : 0 : call MPI_FILE_CLOSE(mpi_fh, mpi_err)
770 : 0 : ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
771 : : #else
772 : : ABI_ERROR("You should not be here")
773 : : #endif
774 : : end if
775 : :
776 : 0 : call cwtime(cputime, walltime, gflops, "stop")
777 : 0 : write(msg,'(2(a,f9.1),a)')" IO operation completed. cpu_time: ",cputime,"[s], walltime: ",walltime," [s]"
778 : 0 : call wrtout(std_out, msg, "COLL", do_flush=.True.)
779 : : end if
780 : :
781 : 1 : return
782 : :
783 : : ! Handle IO-error
784 : : 10 continue
785 : 0 : ABI_ERROR(errmsg)
786 : :
787 [ - - + - ]: 1 : end subroutine exc_init_phi_block
788 : : !!***
789 : :
790 : : !----------------------------------------------------------------------
791 : :
792 : : !!****f* m_exc_itdiago/exc_write_phi_block
793 : : !! NAME
794 : : !! exc_write_phi_block
795 : : !!
796 : : !! FUNCTION
797 : : !! Write phi_block on the Fortran file oeig_fname.
798 : : !!
799 : : !! SOURCE
800 : :
801 : 1 : subroutine exc_write_phi_block(oeig_fname,use_mpio)
802 : :
803 : : !Arguments ------------------------------------
804 : : !scalars
805 : : character(len=*),intent(in) :: oeig_fname
806 : : logical,intent(in) :: use_mpio
807 : :
808 : : !Local variables ------------------------------
809 : : integer :: eig_unt,state,mpi_err !,fform
810 : : real(dp) :: cputime,walltime,gflops
811 : : character(len=500) :: msg,errmsg
812 : : ! type(Hdr_type) :: hexc_Hdr
813 : : logical :: do_ep_lifetime
814 : : !!arrays
815 : 1 : complex(dp),allocatable :: buffer_dpc(:)
816 : : #ifdef HAVE_MPI_IO
817 : : integer:: amode,mpi_fh,old_type,etype,eig_type,my_nel,ierr,my_nrows
818 : : integer(XMPI_OFFSET_KIND) :: ehdr_offset,my_offset,my_offpad,fmarker
819 : 1 : integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecord(:)
820 : : integer :: array_of_sizes(2),array_of_subsizes(2),array_of_starts(2)
821 : : #endif
822 : : !************************************************************************
823 : :
824 : 1 : do_ep_lifetime = .FALSE.
825 : :
826 : 1 : call cwtime(cputime, walltime, gflops, "start")
827 : :
828 [ + - ]: 1 : if (.not.use_mpio) then
829 : :
830 : : ! * Master writes the header.
831 [ + - ]: 1 : if (my_rank==master) then
832 [ + - ]: 1 : call wrtout(std_out," Writing eigenstates on file "//TRIM(oeig_fname)//" via Fortran-IO","COLL")
833 [ - + ]: 1 : if (open_file(oeig_fname,msg, newunit=eig_unt, form='unformatted') /= 0) then
834 : 0 : ABI_ERROR(msg)
835 : : end if
836 [ + - ]: 1 : write(eig_unt, err=10, iomsg=errmsg) do_ep_lifetime
837 [ + - ]: 1 : write(eig_unt, err=10, iomsg=errmsg) hexc_size, nstates
838 [ + + + - ]: 11 : write(eig_unt, err=10, iomsg=errmsg) CMPLX(exc_energy(1:nstates),kind=dp)
839 : : end if
840 : :
841 : : ! Wavefunctions are gathered on the master node band-by-band.
842 : : ! TODO bands should be treated in blocks to minimize the number of MPI calls.
843 [ - + + - : 3 : ABI_MALLOC_OR_DIE(buffer_dpc,(hexc_size), ierr)
+ - + - +
- + - - +
+ - - + ]
844 : :
845 [ + + ]: 11 : do state=1,nstates
846 [ + + ]: 3850 : buffer_dpc=czero
847 [ + + ]: 3850 : buffer_dpc(my_t1:my_t2) = phi_block(:,state)
848 : 10 : call xmpi_sum_master(buffer_dpc,master,comm,mpi_err)
849 [ + - + - ]: 11 : if (my_rank==master) write(eig_unt, err=10, iomsg=errmsg) buffer_dpc(1:hexc_size)
850 : : end do
851 [ - + ]: 1 : ABI_FREE(buffer_dpc)
852 : :
853 [ + - + - ]: 1 : if (my_rank==master) close(eig_unt, err=10, iomsg=errmsg)
854 : :
855 : : else
856 : : #ifdef HAVE_MPI_IO
857 [ # # ]: 0 : call wrtout(std_out," Writing eigenstates on file "//TRIM(oeig_fname)//" with MPI-IO","COLL")
858 : :
859 : : ! Write the header.
860 [ # # ]: 0 : if (my_rank==master) then
861 : : ! Write header using Fortran primitives.
862 [ # # ]: 0 : if (open_file(oeig_fname,msg,newunit=eig_unt,form='unformatted') /= 0) then
863 : 0 : ABI_ERROR(msg)
864 : : end if
865 [ # # ]: 0 : write(eig_unt, err=10, iomsg=errmsg) nstates
866 [ # # # # ]: 0 : write(eig_unt, err=10, iomsg=errmsg) CMPLX(exc_energy(1:nstates),kind=dp)
867 : : ! TODO: change setup_bse so that Hdr_bse reflects the parameters of the run.
868 [ # # ]: 0 : close(eig_unt, err=10, iomsg=errmsg)
869 : : end if
870 : :
871 : 0 : call xmpi_barrier(comm)
872 : :
873 : : ! Open the file with MPI-IO
874 : 0 : amode=MPI_MODE_RDWR
875 : :
876 : 0 : call MPI_FILE_OPEN(comm, oeig_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
877 [ # # ]: 0 : msg = " MPI_IO error opening file: "//TRIM(oeig_fname)
878 : 0 : ABI_CHECK_MPI(mpi_err,msg)
879 : :
880 : : ! Move the file pointer to skip the first two records.
881 : 0 : ehdr_offset=0
882 : 0 : call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err)
883 : : !write(std_out,*)"fmarker first record ",fmarker
884 : 0 : call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err)
885 : : !write(std_out,*)"fmarker first record ",fmarker
886 : : !$call hdr_mpio_skip(mpi_fh,fform,ehdr_offset)
887 : : !$ehdr_offset = 4*xmpio_bsize_frm + xmpio_bsize_int + nstates*xmpio_bsize_dpc
888 : :
889 : 0 : etype=MPI_BYTE; old_type=MPI_DOUBLE_COMPLEX
890 : :
891 : 0 : my_nrows=my_t2-my_t1+1; old_type=MPI_DOUBLE_COMPLEX
892 [ # # ]: 0 : array_of_sizes = (/hexc_size,nstates/)
893 [ # # ]: 0 : array_of_subsizes = (/my_nrows,nstates/)
894 [ # # ]: 0 : array_of_starts = (/my_t1,1/)
895 : 0 : call xmpio_create_fsubarray_2D(array_of_sizes,array_of_subsizes,array_of_starts,old_type,eig_type,my_offpad,mpi_err)
896 : 0 : ABI_CHECK_MPI(mpi_err,"fsubarray_2D")
897 : : !
898 : : ! Each node uses a different offset to skip the header and the blocks written by the other CPUs.
899 : 0 : my_offset = ehdr_offset + my_offpad
900 : :
901 : 0 : call MPI_FILE_SET_VIEW(mpi_fh, my_offset, etype, eig_type, 'native', MPI_INFO_NULL, mpi_err)
902 : 0 : ABI_CHECK_MPI(mpi_err,"SET_VIEW")
903 : :
904 : 0 : call MPI_TYPE_FREE(eig_type,mpi_err)
905 : 0 : ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE")
906 : :
907 : 0 : my_nel = my_nrows*nstates
908 : 0 : call MPI_FILE_WRITE_ALL(mpi_fh, phi_block, my_nel, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err)
909 : 0 : ABI_CHECK_MPI(mpi_err,"FILE_WRITE")
910 : :
911 : : ! It seems that personal calls make the code stuck
912 [ # # # # : 0 : ABI_MALLOC(bsize_frecord,(nstates))
# # # # #
# # # ]
913 [ # # ]: 0 : bsize_frecord = hexc_size * xmpi_bsize_dpc
914 : : ! ehdr_offset points to the end of the header.
915 : 0 : call xmpio_write_frmarkers(mpi_fh,ehdr_offset,xmpio_collective,nstates,bsize_frecord,ierr)
916 [ # # ]: 0 : ABI_CHECK(ierr==0,"Error while writing Fortran markers")
917 [ # # ]: 0 : ABI_FREE(bsize_frecord)
918 : : !
919 : : ! Close the file.
920 : 0 : call MPI_FILE_CLOSE(mpi_fh, mpi_err)
921 : 0 : ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
922 : : #else
923 : : ABI_ERROR("MPI-IO support not enabled")
924 : : #endif
925 : : end if
926 : :
927 : 1 : call cwtime(cputime, walltime, gflops, "stop")
928 : 1 : write(msg,'(2(a,f9.1),a)')" IO operation completed. cpu_time: ",cputime,"[s], walltime: ",walltime," [s]"
929 : 1 : call wrtout(std_out, msg, "COLL", do_flush=.True.)
930 : :
931 : 1 : return
932 : :
933 : : ! Handle IO-error
934 : : 10 continue
935 : 0 : ABI_ERROR(errmsg)
936 : :
937 [ + - - - ]: 1 : end subroutine exc_write_phi_block
938 : : !!***
939 : :
940 : : !----------------------------------------------------------------------
941 : :
942 : : !!****f* m_exc_itdiago/exc_subspace_rotation
943 : : !! NAME
944 : : !! exc_subspace_rotation
945 : : !!
946 : : !! FUNCTION
947 : : !! This routine performs the subspace rotation.
948 : : !!
949 : : !! SIDE EFFECTS
950 : : !! phi_block(my_t1:my_t2,nstates)=Contains the trial eigenstates.
951 : : !!
952 : : !! SOURCE
953 : :
954 : 21 : subroutine exc_subspace_rotation()
955 : :
956 : : !Local variables ------------------------------
957 : : integer :: ii,jj,ipack,ierr
958 : : !arrays
959 : 21 : real(dp),allocatable :: sub_ene(:)
960 : : ! real(dp),allocatable :: evec(:,:)
961 : 21 : complex(dp),allocatable :: sub_ham(:,:),sub_pham(:),hphi_tot(:)
962 : : ! complex(dp),allocatable :: phi_tmp(:,:)
963 : : !************************************************************************
964 : :
965 : : ! * Sub-space rotation. Calculate <phi_i|H|phi_j> in packed form.
966 : : ! TODO: this part can be rewritten using BLAS3 routines.
967 : :
968 [ - + + - : 63 : ABI_MALLOC(hphi_tot,(hexc_size))
+ - + - -
+ - + ]
969 [ + - - + : 63 : ABI_MALLOC(sub_pham,(nstates*(nstates+1)/2))
+ - + - +
- - + -
+ ]
970 [ + + ]: 1176 : sub_pham=czero; ipack=0
971 : :
972 [ + + ]: 231 : do jj=1,nstates
973 [ + + ]: 80850 : hphi_tot = czero
974 [ + + + + : 124024530 : hphi_tot(:) = MATMUL(hexc, phi_block(:,jj))
+ + + - +
- + - +
- ]
975 : 210 : call xmpi_sum(hphi_tot,comm,ierr)
976 : :
977 [ + + ]: 1386 : do ii=1,jj
978 : 1155 : ipack=ipack+1
979 [ + + ]: 444675 : sub_pham(ipack) = DOT_PRODUCT(phi_block(my_t1:my_t2,ii), hphi_tot(my_t1:my_t2) )
980 [ + + ]: 1365 : if (ii==jj) sub_pham(ipack) = REAL(sub_pham(ipack),kind=dp)
981 : : end do
982 : : end do
983 : 21 : call xmpi_sum(sub_pham,comm,ierr)
984 : :
985 [ - + + - : 84 : ABI_MALLOC(sub_ham,(nstates,nstates))
- + + - +
- + - - +
- + ]
986 [ + + + + ]: 2331 : sub_ham=czero
987 [ - + + - : 63 : ABI_MALLOC(sub_ene,(nstates))
+ - + - -
+ - + ]
988 : :
989 : 21 : call xhpev("Vectors","Upper",nstates,sub_pham,sub_ene,sub_ham,nstates) !,comm)
990 : :
991 [ - + ]: 21 : ABI_FREE(hphi_tot)
992 [ - + ]: 21 : ABI_FREE(sub_pham)
993 [ - + ]: 21 : ABI_FREE(sub_ene)
994 : :
995 : : !do ii=1,nstates
996 : : ! norm = DOT_PRODUCT(sub_ham(:,ii),sub_ham(:,ii))
997 : : ! write(std_out,*)"norm subspac",norm
998 : : ! sub_ham(:,ii) = sub_ham(:,ii)/norm
999 : : !end do
1000 : :
1001 : : !allocate(evec(2*nstates,nstates))
1002 : :
1003 : : !do ii=1,nstates
1004 : : ! do jj=1,nstates
1005 : : ! evec(jj, ii) = REAL (sub_ham(jj,ii))
1006 : : ! evec(jj+1,ii) = AIMAG(sub_ham(jj,ii))
1007 : : ! end do
1008 : : !end do
1009 : :
1010 : : !call normev(evec,nstates,nstates)
1011 : :
1012 : : !do ii=1,nstates
1013 : : ! do jj=1,nstates
1014 : : ! sub_ham(jj,ii) = CMPLX( evec(jj,ii),evec(jj+1,ii) )
1015 : : ! end do
1016 : : !end do
1017 : : !deallocate(evec)
1018 : :
1019 : : #if 0
1020 : : ABI_MALLOC(phi_tmp,(my_nt,nstates))
1021 : : phi_tmp = phi_block
1022 : :
1023 : : call ZGEMM('N','N',my_nt,nstates,nstates,cone,phi_tmp,my_nt,sub_ham,nstates,czero,phi_block,my_nt)
1024 : :
1025 : : ABI_FREE(phi_tmp)
1026 : : #else
1027 [ + - + - : 161784 : phi_block = MATMUL(phi_block,sub_ham)
+ + + + -
+ - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - +
- + - - +
- - + + +
+ + - ]
1028 : : #endif
1029 : :
1030 [ - + ]: 21 : ABI_FREE(sub_ham)
1031 : :
1032 : 21 : end subroutine exc_subspace_rotation
1033 : : !!***
1034 : :
1035 : : !----------------------------------------------------------------------
1036 : :
1037 : : !!****f* m_exc_itdiago/exc_cholesky_ortho
1038 : : !! NAME
1039 : : !! exc_cholesky_ortho
1040 : : !!
1041 : : !! FUNCTION
1042 : : !! This routine performs the orthogonalization of the trial eigenstates using the
1043 : : !! Cholesky Algorithm.
1044 : : !!
1045 : : !! SIDE EFFECTS
1046 : : !! phi_block(my_t1:my_t2,nstates)=Contains the trial eigenstates.
1047 : : !!
1048 : : !! SOURCE
1049 : :
1050 : 21 : subroutine exc_cholesky_ortho()
1051 : :
1052 : : !Local variables ------------------------------
1053 : : integer :: my_info,ii,jj,ipack,ierr
1054 : : logical,parameter :: use_unpacked = .False.
1055 : : !arrays
1056 : 21 : complex(dp),allocatable :: overlap(:,:),povlp(:)
1057 : : !************************************************************************
1058 : :
1059 : : ! 1) overlap_ij = <phi_i|phi_j>
1060 [ - + + - : 84 : ABI_MALLOC(overlap, (nstates,nstates))
- + + - +
- + - - +
- + ]
1061 : :
1062 : : if (use_unpacked) then
1063 : : overlap = czero
1064 : :
1065 : : call ZGEMM('C','N',nstates,nstates,my_nt,cone,phi_block,my_nt,phi_block,my_nt,czero,overlap,nstates)
1066 : : call xmpi_sum(overlap,comm,ierr)
1067 : :
1068 : : do ii=1,nstates
1069 : : overlap(ii,ii)=REAL(overlap(ii,ii),kind=dp)
1070 : : end do
1071 : :
1072 : : ! 2) Cholesky factorization: overlap = U^H U with U upper triangle matrix.
1073 : : call ZPOTRF('U',nstates,overlap,nstates,my_info)
1074 : : if (my_info/=0) then
1075 : : write(msg,'(a,i3)')' ZPOTRF returned info= ',my_info
1076 : : ABI_ERROR(msg)
1077 : : end if
1078 : :
1079 : : else
1080 : : ! 1) Calculate overlap_ij = <phi_i|phi_j> in packed form.
1081 [ + - - + : 63 : ABI_MALLOC(povlp,(nstates*(nstates+1)/2))
+ - + - +
- - + -
+ ]
1082 [ + + ]: 1176 : povlp = czero; ipack=0
1083 [ + + ]: 231 : do jj=1,nstates
1084 [ + + ]: 1386 : do ii=1,jj
1085 : 1155 : ipack=ipack+1
1086 [ + + ]: 444675 : povlp(ipack) = DOT_PRODUCT( phi_block(my_t1:my_t2,ii), phi_block(my_t1:my_t2,jj) )
1087 [ + + ]: 1365 : if (ii==jj) povlp(ipack) = REAL(povlp(ipack),kind=dp)
1088 : : end do
1089 : : end do
1090 : 21 : call xmpi_sum(povlp,comm,ierr)
1091 : :
1092 : : ! 2) Cholesky factorization: overlap = U^H U with U upper triangle matrix.
1093 : 21 : call ZPPTRF("U",nstates,povlp,my_info)
1094 [ - + ]: 21 : if (my_info/=0) then
1095 : 0 : write(msg,'(a,i3)')' ZPPTRF returned info= ',my_info
1096 : 0 : ABI_ERROR(msg)
1097 : : end if
1098 : : !call xmpi_sum(povlp,comm,ierr)
1099 : : !povlp=povlp/nproc
1100 : :
1101 : : !unpack povlp to prepare call to ZTRSM.
1102 : 21 : ipack=0
1103 [ + + ]: 231 : do jj=1,nstates
1104 [ + + ]: 1386 : do ii=1,jj
1105 : 1155 : ipack=ipack+1
1106 [ + + ]: 1365 : if (ii/=jj) then
1107 : 945 : overlap(ii,jj)= povlp(ipack)
1108 : 945 : overlap(jj,ii)=CONJG(povlp(ipack))
1109 : : else
1110 : 210 : overlap(ii,ii)=REAL(povlp(ipack),kind=dp)
1111 : : end if
1112 : : end do
1113 : : end do
1114 [ - + ]: 42 : ABI_FREE(povlp)
1115 : : end if
1116 : :
1117 : : ! Check if this can be done with Scalapack. Direct PZTRSM is not provided
1118 : :
1119 : : ! 3) Solve X U = phi_block, on exit the phi_block treated by this node is orthonormalized.
1120 : : !call ZTRSM('R','U','N','N',hexc_size,nstates,cone,overlap,nstates,phi_block,hexc_size)
1121 : 21 : call ZTRSM('Right','Upper','Normal','Normal',my_nt,nstates,cone,overlap,nstates,phi_block,my_nt)
1122 [ - + ]: 21 : ABI_FREE(overlap)
1123 : :
1124 : 21 : end subroutine exc_cholesky_ortho
1125 : : !!***
1126 : :
1127 : : !----------------------------------------------------------------------
1128 : :
1129 : : !!****f* m_exc_itdiago/convergence_degree
1130 : : !! NAME
1131 : : !! convergence_degree
1132 : : !!
1133 : : !! FUNCTION
1134 : : !! Return the degree of convergence from the input residual.
1135 : : !!
1136 : : !! INPUTS
1137 : : !! resid=Residual.
1138 : : !!
1139 : : !! OUTPUT
1140 : : !!
1141 : : !! SOURCE
1142 : :
1143 : 736 : function convergence_degree(resid)
1144 : :
1145 : : !Arguments
1146 : : integer :: convergence_degree
1147 : : real(dp),intent(in) :: resid
1148 : :
1149 : : !************************************************************************
1150 : :
1151 [ + + ]: 736 : if (resid<tolwfr_) then
1152 : : convergence_degree = STRICT
1153 : : else
1154 : 713 : convergence_degree = WORST
1155 [ + + ]: 713 : if (resid<tolwfr_*10**5) convergence_degree = MEDIUM
1156 : : end if
1157 : :
1158 : : end function convergence_degree
1159 : : !!***
1160 : :
1161 : : !----------------------------------------------------------------------
1162 : :
1163 : : !!****f* m_exc_itdiago/exc_check_phi_block
1164 : : !! NAME
1165 : : !! exc_check_phi_block
1166 : : !!
1167 : : !! FUNCTION
1168 : : !! Debugging tools
1169 : : !!
1170 : : !! INPUTS
1171 : : !!
1172 : : !! OUTPUT
1173 : : !!
1174 : : !! SOURCE
1175 : :
1176 : : subroutine exc_check_phi_block(string)
1177 : :
1178 : : !Arguments ------------------------------------
1179 : : !scalars
1180 : : character(len=*),intent(in) :: string
1181 : :
1182 : : !Local variables ------------------------------
1183 : : !scalars
1184 : : integer :: ii,jj,ierr
1185 : : real(dp) :: err,rdum
1186 : : !arrays
1187 [ - + - + ]: 22 : complex(dp),allocatable :: lbuff(:,:)
1188 : : !************************************************************************
1189 : :
1190 : : if (.not. DEBUGME) return
1191 : :
1192 : : #if 0
1193 : : ABI_MALLOC(lbuff,(hexc_size,nstates))
1194 : : err = -one
1195 : : do irank=1,nproc-1
1196 : : call xmpi_exch(phi_block,hexc_size*nstates,irank,lbuff,master,comm,11,ierr)
1197 : : if (my_rank==master) then
1198 : : lbuff = lbuff-phi_block
1199 : : err = MAX(err,MAXVAL(MAXVAL(ABS(lbuff),DIM=1)))
1200 : : end if
1201 : : call xmpi_barrier(comm)
1202 : : end do
1203 : : ABI_FREE(lbuff)
1204 : : #else
1205 : :
1206 : : ABI_MALLOC(lbuff,(nstates,nstates))
1207 : : lbuff=czero
1208 : : do jj=1,nstates
1209 : : do ii=1,jj
1210 : : lbuff(ii,jj) = DOT_PRODUCT( phi_block(my_t1:my_t2,ii), phi_block(my_t1:my_t2,jj) )
1211 : : end do
1212 : : end do
1213 : : call xmpi_sum(lbuff,comm,ierr)
1214 : :
1215 : : err = -one
1216 : : do jj=1,nstates
1217 : : do ii=1,jj
1218 : : if (ii==jj) then
1219 : : rdum = ABS(lbuff(ii,jj)-one)
1220 : : else
1221 : : rdum = ABS(lbuff(ii,jj))
1222 : : end if
1223 : : err = MAX(err,rdum)
1224 : : end do
1225 : : end do
1226 : : ABI_FREE(lbuff)
1227 : : #endif
1228 : :
1229 : : if (my_rank==master) then
1230 : : write(std_out,*)"After ",TRIM(string),", MAX inconsistency error in phi_block= ",err
1231 : : end if
1232 : :
1233 : : !write(std_out,*)"master casts its own data"
1234 : : !call xmpi_bcast(phi_block,master,comm,ierr)
1235 : :
1236 : : end subroutine exc_check_phi_block
1237 : : !!***
1238 : :
1239 : : !----------------------------------------------------------------------
1240 : :
1241 : : end subroutine exc_iterative_diago
1242 : : !!***
1243 : :
1244 : : end module m_exc_itdiago
1245 : : !!***
|
