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