Branch data Line data Source code
1 : : !!****m* ABINIT/m_rayleigh_ritz
2 : : !! NAME
3 : : !! m_rayleigh_ritz
4 : : !!
5 : : !! FUNCTION
6 : : !! This file contains routines that perform the Rayleigh-Ritz,
7 : : !! either by forming the full matrix (_subdiago) or by forming the
8 : : !! distributed matrix in block-cyclic form (_distributed)
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2014-2026 ABINIT group (AL)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! SOURCE
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : :
22 : : #include "abi_common.h"
23 : :
24 : : module m_rayleigh_ritz
25 : :
26 : : use defs_basis
27 : : use m_errors
28 : : use m_cgtools
29 : : use m_xmpi
30 : : use m_abicore
31 : : use m_abi_linalg
32 : : use m_slk
33 : :
34 : : use defs_abitypes, only : mpi_type
35 : : use m_time, only : timab
36 : : use m_numeric_tools, only : pack_matrix
37 : :
38 : : implicit none
39 : :
40 : : private
41 : : !!***
42 : :
43 : : public :: rayleigh_ritz_subdiago
44 : : #if defined HAVE_LINALG_SCALAPACK
45 : : public :: rayleigh_ritz_distributed
46 : : #endif
47 : : !!***
48 : :
49 : : contains
50 : : !!***
51 : :
52 : : !!****f* ABINIT/rayleigh_ritz_subdiago
53 : : !! NAME
54 : : !! rayleigh_ritz_subdiago
55 : : !!
56 : : !! FUNCTION
57 : : !! Performs a rayleigh-ritz procedure (subspace rotation), building the
58 : : !! hamiltonian/overlap matrices in full and calling the subdiago method
59 : : !!
60 : : !! INPUTS
61 : : !! mpi_enreg=information about MPI parallelization
62 : : !! nband=number of bands at this k point for that spin polarization
63 : : !! npw=number of plane waves at this k point
64 : : !! nspinor=number of plane waves at this k point
65 : : !! usepaw=if 1 we use the PAW method
66 : : !!
67 : : !! OUTPUT
68 : : !! eig(nband)=array for holding eigenvalues (hartree)
69 : : !!
70 : : !! SIDE EFFECTS
71 : : !! cg(2,*)=updated wavefunctions
72 : : !! ghc(2,*)=updated ghc
73 : : !! gsc(2,*)=updated gsc
74 : : !! gvnlxc(2,*)=updated gvnlxc
75 : : !!
76 : : !! NOTES
77 : : !! TODO choose generalized eigenproblem or ortho + diago (see #if 1)
78 : : !!
79 : : !! SOURCE
80 : :
81 : 48 : subroutine rayleigh_ritz_subdiago(cg,ghc,gsc,gvnlxc,eig,has_fock,istwf_k,mpi_enreg,nband,npw,nspinor,usepaw)
82 : :
83 : : ! Arguments
84 : : type(mpi_type),intent(inout) :: mpi_enreg
85 : : integer,intent(in) :: nband,npw,nspinor,usepaw,istwf_k
86 : : real(dp),intent(inout) :: cg(2,npw*nspinor*nband),gsc(2,npw*nspinor*nband),ghc(2,npw*nspinor*nband),gvnlxc(2,npw*nspinor*nband)
87 : : real(dp),intent(out) :: eig(nband)
88 : : logical :: has_fock
89 : :
90 : : ! Locals
91 : 48 : real(dp), allocatable :: subham(:), totham(:,:)
92 : 48 : real(dp), allocatable :: subovl(:), totovl(:,:)
93 : 48 : real(dp), allocatable :: evec(:,:), edummy(:,:)
94 : : integer :: ierr, cplx, vectsize
95 : 96 : real(dp) :: gtempc(2,npw*nspinor*nband), tsec(2)
96 : : character :: blas_transpose
97 : :
98 : : integer, parameter :: timer_chebfi = 1600, timer_alltoall = 1601, timer_apply_inv_ovl = 1602, timer_rotation = 1603
99 : : integer, parameter :: timer_subdiago = 1604, timer_subham = 1605, timer_ortho = 1606, timer_getghc = 1607
100 : : integer, parameter :: timer_residuals = 1608, timer_update_eigen = 1609, timer_sync = 1610
101 : :
102 : : ! *************************************************************************
103 : :
104 [ + - ]: 48 : if(istwf_k == 1) then
105 : 48 : cplx = 2
106 : 48 : vectsize = npw*nspinor
107 : 48 : blas_transpose = 'c'
108 : : else
109 : 0 : cplx = 1
110 : 0 : vectsize = 2*npw*nspinor
111 : 0 : blas_transpose = 't'
112 : : end if
113 : :
114 : : #if 1
115 : 48 : call timab(timer_subham, 1, tsec)
116 : :
117 : : ! Transform cg, ghc and maybe gsc, according to istwf_k
118 [ - + ]: 48 : if(istwf_k == 2) then
119 [ # # # # ]: 0 : cg = cg * sqrt2
120 [ # # # # : 0 : if(mpi_enreg%me_g0 == 1) cg(:, 1:npw*nspinor*nband:npw) = cg(:, 1:npw*nspinor*nband:npw) / sqrt2
# # ]
121 [ # # # # ]: 0 : ghc = ghc * sqrt2
122 [ # # # # : 0 : if(mpi_enreg%me_g0 == 1) ghc(:, 1:npw*nspinor*nband:npw) = ghc(:, 1:npw*nspinor*nband:npw) / sqrt2
# # ]
123 [ # # ]: 0 : if(usepaw == 1) then
124 [ # # # # ]: 0 : gsc = gsc * sqrt2
125 [ # # # # : 0 : if(mpi_enreg%me_g0 == 1) gsc(:, 1:npw*nspinor*nband:npw) = gsc(:, 1:npw*nspinor*nband:npw) / sqrt2
# # ]
126 : : end if
127 : : end if
128 : :
129 : : ! Build, pack and sum suham
130 [ + - - + : 144 : ABI_MALLOC(subham, (cplx*nband*(nband+1)/2))
+ - + - +
- - + -
+ ]
131 [ - + + - : 192 : ABI_MALLOC(totham, (cplx, nband*nband))
- + + - +
- + - - +
- + ]
132 : : call abi_xgemm(blas_transpose,'n',nband,nband,vectsize,cone,ghc,vectsize,&
133 : 48 : & cg,vectsize,czero,totham,nband, x_cplx=cplx)
134 : 48 : call pack_matrix(totham, subham, nband, cplx)
135 [ - + ]: 48 : ABI_FREE(totham)
136 : 48 : call xmpi_sum(subham,mpi_enreg%comm_bandspinorfft,ierr)
137 : :
138 : :
139 : : ! Same for subovl
140 [ - + + - : 96 : ABI_MALLOC(subovl, (cplx*nband*(nband+1)/2))
- + - + ]
141 [ - + + - : 144 : ABI_MALLOC(totovl, (cplx, nband*nband))
- + + - -
+ - + ]
142 [ + - ]: 48 : if(usepaw == 1) then
143 : : call abi_xgemm(blas_transpose,'n',nband,nband,vectsize,cone,gsc,vectsize,&
144 : 48 : & cg,vectsize,czero,totovl,nband, x_cplx=cplx)
145 : : else
146 : : call abi_xgemm(blas_transpose,'n',nband,nband,vectsize,cone,cg,vectsize,&
147 : 0 : & cg,vectsize,czero,totovl,nband, x_cplx=cplx)
148 : : end if
149 : 48 : call pack_matrix(totovl, subovl, nband, cplx)
150 [ - + ]: 48 : ABI_FREE(totovl)
151 : 48 : call xmpi_sum(subovl,mpi_enreg%comm_bandspinorfft,ierr)
152 : :
153 : :
154 : : ! Transform back
155 [ - + ]: 48 : if(istwf_k == 2) then
156 [ # # # # ]: 0 : cg = cg / sqrt2
157 [ # # # # : 0 : if(mpi_enreg%me_g0 == 1) cg(:, 1:npw*nspinor*nband:npw) = cg(:, 1:npw*nspinor*nband:npw) * sqrt2
# # ]
158 [ # # # # ]: 0 : ghc = ghc / sqrt2
159 [ # # # # : 0 : if(mpi_enreg%me_g0 == 1) ghc(:, 1:npw*nspinor*nband:npw) = ghc(:, 1:npw*nspinor*nband:npw) * sqrt2
# # ]
160 [ # # ]: 0 : if(usepaw == 1) then
161 [ # # # # ]: 0 : gsc = gsc / sqrt2
162 [ # # # # : 0 : if(mpi_enreg%me_g0 == 1) gsc(:, 1:npw*nspinor*nband:npw) = gsc(:, 1:npw*nspinor*nband:npw) * sqrt2
# # ]
163 : : end if
164 : : end if
165 : 48 : call timab(timer_subham, 2, tsec)
166 : :
167 : :
168 : 48 : call timab(timer_subdiago, 1, tsec)
169 [ - + + - : 192 : ABI_MALLOC(evec, (cplx*nband, nband))
- + + - +
- + - - +
- + ]
170 : :
171 : 48 : call abi_xhpgv(1,'V','U',nband,subham,subovl,eig,evec,nband,istwf_k=istwf_k,use_slk=mpi_enreg%paral_kgb)
172 : :
173 [ - + ]: 48 : ABI_FREE(subham)
174 [ - + ]: 48 : ABI_FREE(subovl)
175 : :
176 : : ! Fix the phase (this is because of the simultaneous diagonalisation of this
177 : : ! matrix by different processors, allowing to get different unitary transforms, thus breaking the
178 : : ! coherency of parts of cg stored on different processors).
179 : : ! call cg_normev(evec,nband,nband) ! Unfortunately, for cg_normev to work, one needs the vectors to be normalized, so uses fxphas_seq
180 [ - + + - : 144 : ABI_MALLOC(edummy, (cplx*nband, nband))
- + + - -
+ - + ]
181 : 48 : call fxphas_seq(evec,edummy,0,0,1,nband*nband,nband*nband,nband,nband,0)
182 [ - + ]: 48 : ABI_FREE(edummy)
183 : :
184 : : ! Rotate
185 : 48 : call abi_xgemm('n','n',vectsize,nband, nband,cone,cg , vectsize, evec, nband, czero, gtempc, vectsize, x_cplx=cplx)
186 [ + + + + ]: 111936 : cg = gtempc
187 : 48 : call abi_xgemm('n','n',vectsize,nband, nband,cone,ghc, vectsize, evec, nband, czero, gtempc, vectsize, x_cplx=cplx)
188 [ + + + + ]: 111936 : ghc = gtempc
189 [ + - ]: 48 : if(usepaw == 1) then
190 : 48 : call abi_xgemm('n','n',vectsize,nband, nband,cone,gsc, vectsize, evec, nband, czero, gtempc, vectsize, x_cplx=cplx)
191 [ + + + + ]: 111936 : gsc = gtempc
192 : : endif
193 [ + - - + ]: 48 : if(usepaw==0 .or. has_fock)then
194 : 0 : call abi_xgemm('n','n',vectsize,nband, nband,cone,gvnlxc, vectsize, evec, nband, czero, gtempc, vectsize, x_cplx=cplx)
195 [ - - - - ]: 48 : gvnlxc = gtempc
196 : : end if
197 [ - + ]: 48 : ABI_FREE(evec)
198 : 48 : call timab(timer_subdiago, 2, tsec)
199 : :
200 : : #else
201 : : !! TODO non-functional, should be rewritten. Possibly faster (tests needed)
202 : : call wrtout(std_out, 'Transposed, orthogonalizing')
203 : :
204 : : ! orthonormalization
205 : : call timab(timer_ortho, 1, tsec)
206 : : if (usepaw==1) then
207 : : call abi_xorthonormalize(cg, gsc,nband, mpi_enreg%comm_bandspinorfft, sqgram, npw*nspinor, 2)
208 : : else
209 : : call abi_xorthonormalize(cg, cg, nband, mpi_enreg%comm_bandspinorfft, sqgram, npw*nspinor, 2)
210 : : end if
211 : : call timab(timer_ortho, 2, tsec)
212 : :
213 : : ! rotate ghc, gsc and gvnlxc
214 : : call timab(timer_rotation, 1, tsec)
215 : : call abi_xtrsm('r','u','n','n',npw*nspinor,nband,cone,sqgram,nband, ghc,npw*nspinor,x_cplx=2)
216 : : if(usepaw==1) then
217 : : call abi_xtrsm('r','u','n','n',npw*nspinor,nband,cone,sqgram,nband, gsc,npw*nspinor,x_cplx=2)
218 : : endif
219 : : if(usepaw==0 .or has_fock)then
220 : : call abi_xtrsm('r','u','n','n',npw*nspinor,nband,cone,sqgram,nband, gvnlxc,npw*nspinor,x_cplx=2)
221 : : end if
222 : : call timab(timer_rotation, 2, tsec)
223 : :
224 : : call wrtout(std_out, 'Orthogonalized, building subham')
225 : :
226 : : ! build hamiltonian in subspace
227 : : call timab(timer_subham, 1, tsec)
228 : : call abi_xgemm(blas_transpose,'n',nband,nband,npw*nspinor,cone,ghc,npw*nspinor,&
229 : : & cg,npw*nspinor,czero,totham,nband, x_cplx=2)
230 : : ! pack for compatibility with subdiago
231 : : call pack_matrix(totham, subham, nband)
232 : : call xmpi_sum(subham,mpi_enreg%comm_bandspinorfft,ierr)
233 : : call timab(timer_subham, 2, tsec)
234 : :
235 : : call wrtout(std_out, 'Subham built, diagonalizing')
236 : :
237 : : ! Rayleigh-Ritz
238 : : call timab(timer_subdiago,1,tsec)
239 : : call subdiago(cg,eig,evec,gsc,0,0,gs_hamk%istwf_k,&
240 : : & mcg,mcg,nband,npw,nspinor,dtset%paral_kgb,&
241 : : & subham,dummy,0,gs_hamk%usepaw,mpi_enreg%me_g0)
242 : : call timab(timer_subdiago,2,tsec)
243 : :
244 : : call wrtout(std_out, 'Diagonalization done')
245 : :
246 : : ! Rotate ghc and gvnlxc according to evecs
247 : : call timab(timer_rotation, 1, tsec)
248 : : call abi_xgemm('n','n',npw*nspinor,nband, nband,cone,ghc, npw*nspinor, evec, nband, czero, gtempc, npw*nspinor, x_cplx=2)
249 : : ghc = gtempc
250 : : if(usepaw==0 .or has_fock)then
251 : : call abi_xgemm('n','n',npw*nspinor,nband, nband,cone,gvnlxc, npw*nspinor, evec, nband, czero, gtempc, npw*nspinor, x_cplx=2)
252 : : gvnlxc = gtempc
253 : : end if
254 : : call timab(timer_rotation, 2, tsec)
255 : :
256 : : #endif
257 : :
258 : 96 : end subroutine rayleigh_ritz_subdiago
259 : : !!***
260 : :
261 : : #if defined HAVE_LINALG_SCALAPACK
262 : : !!****f* ABINIT/rayleigh_ritz_distributed
263 : : !! NAME
264 : : !! rayleigh_ritz_distributed
265 : : !!
266 : : !! FUNCTION
267 : : !! Performs a rayleigh-ritz procedure (subspace rotation), building the distributed
268 : : !! hamiltonian/overlap matrices directly, and calling the ScaLapack routines
269 : : !!
270 : : !! INPUTS
271 : : !! mpi_enreg=information about MPI parallelization
272 : : !! nband=number of bands at this k point for that spin polarization
273 : : !! npw=number of plane waves at this k point
274 : : !! nspinor=number of plane waves at this k point
275 : : !! usepaw=do we use the PAW method
276 : : !!
277 : : !! OUTPUT
278 : : !! eig(nband)=array for holding eigenvalues (hartree)
279 : : !!
280 : : !! SIDE EFFECTS
281 : : !! cg(2,*)=updated wavefunctions
282 : : !! ghc(2,*)=updated ghc
283 : : !! gsc(2,*)=updated gsc
284 : : !! gvnlxc(2,*)=updated gvnlxc
285 : : !!
286 : : !! NOTES
287 : : !! Should profile for large test cases and see where the bottleneck is.
288 : : !! Is it the copies? Should we do partial GEMMs?
289 : : !! Is it the latency? Should we buffer more?
290 : : !! Should we overlap computations and communications? (easy in theory, tedious in practice)
291 : : !!
292 : : !! SOURCE
293 : :
294 : : subroutine rayleigh_ritz_distributed(cg,ghc,gsc,gvnlxc,eig,has_fock,istwf_k,mpi_enreg,nband,npw,nspinor,usepaw)
295 : :
296 : : integer,external :: NUMROC
297 : :
298 : : ! Arguments
299 : : type(mpi_type),intent(inout) :: mpi_enreg
300 : : integer,intent(in) :: nband,npw,nspinor,usepaw,istwf_k
301 : : real(dp),intent(inout) :: cg(2,npw*nspinor*nband),gsc(2,npw*nspinor*nband),ghc(2,npw*nspinor*nband),gvnlxc(2,npw*nspinor*nband)
302 : : real(dp),intent(out) :: eig(nband)
303 : : logical :: has_fock
304 : :
305 : : ! Locals
306 : : integer :: blocksize,nbproc,iproc,ierr,cplx,vectsize
307 : : integer :: buffsize_iproc(2), coords_iproc(2), grid_dims(2)
308 : : real(dp) :: cg_new(2,npw*nspinor*nband),gsc_or_vnlxc_new(2,npw*nspinor*nband),ghc_new(2,npw*nspinor*nband)
309 : : real(dp), allocatable :: ham_iproc(:,:), ovl_iproc(:,:), evec_iproc(:,:), left_temp(:,:), right_temp(:,:)
310 : : real(dp) :: tsec(2)
311 : : type(slkmat_dp_t) :: sca_ham, sca_ovl, sca_evec
312 : : !character(len=500) :: message
313 : : character :: blas_transpose
314 : :
315 : : integer, parameter :: timer_chebfi = 1600, timer_alltoall = 1601, timer_apply_inv_ovl = 1602, timer_rotation = 1603
316 : : integer, parameter :: timer_subdiago = 1604, timer_subham = 1605, timer_ortho = 1606, timer_getghc = 1607
317 : : integer, parameter :: timer_residuals = 1608, timer_update_eigen = 1609, timer_sync = 1610
318 : :
319 : : ! *************************************************************************
320 : :
321 : : if(istwf_k == 1) then
322 : : cplx = 2
323 : : vectsize = npw*nspinor
324 : : blas_transpose = 'c'
325 : : else
326 : : cplx = 1
327 : : vectsize = 2*npw*nspinor
328 : : blas_transpose = 't'
329 : : end if
330 : :
331 : : !write(message, *) 'RR: init'
332 : : !call wrtout(std_out,message)
333 : : !======================================================================================================
334 : : ! Init Scalapack matrices
335 : : !======================================================================================================
336 : : call sca_ham%init(nband,nband,slk_processor,istwf_k)
337 : : call sca_ovl%init(nband,nband,slk_processor,istwf_k)
338 : : call sca_evec%init(nband,nband,slk_processor,istwf_k)
339 : :
340 : : ! Get info
341 : : blocksize = sca_ham%size_blocs(1) ! Assume square blocs
342 : : nbproc = slk_processor%grid%nprocs
343 : : grid_dims = slk_processor%grid%dims
344 : :
345 : : !======================================================================================================
346 : : ! Build hamiltonian and overlap matrices
347 : : !======================================================================================================
348 : : ! TODO maybe we should avoid copies at the price of less BLAS efficiency (when blocksize is small)? must profile.
349 : :
350 : : call timab(timer_subham, 1, tsec)
351 : :
352 : : ! Transform cg, ghc and maybe gsc, according to istwf_k
353 : : if(istwf_k == 2) then
354 : : cg = cg * sqrt2
355 : : if(mpi_enreg%me_g0 == 1) cg(:, 1:npw*nspinor*nband:npw) = cg(:, 1:npw*nspinor*nband:npw) / sqrt2
356 : : ghc = ghc * sqrt2
357 : : if(mpi_enreg%me_g0 == 1) ghc(:, 1:npw*nspinor*nband:npw) = ghc(:, 1:npw*nspinor*nband:npw) / sqrt2
358 : : if(usepaw == 1) then
359 : : gsc = gsc * sqrt2
360 : : if(mpi_enreg%me_g0 == 1) gsc(:, 1:npw*nspinor*nband:npw) = gsc(:, 1:npw*nspinor*nband:npw) / sqrt2
361 : : end if
362 : : end if
363 : :
364 : : do iproc=0,nbproc-1
365 : : ! Build the local matrix belonging to processor iproc
366 : : !write(message, *) 'RR: build', iproc
367 : : !call wrtout(std_out,message)
368 : :
369 : : ! Get coordinates of iproc
370 : : coords_iproc(1) = INT(iproc / grid_dims(2))
371 : : coords_iproc(2) = MOD(iproc, grid_dims(2))
372 : :
373 : : ! Get buffersize of iproc
374 : : buffsize_iproc(1) = NUMROC(nband,blocksize,coords_iproc(1),0,slk_processor%grid%dims(1))
375 : : buffsize_iproc(2) = NUMROC(nband,blocksize,coords_iproc(2),0,slk_processor%grid%dims(2))
376 : :
377 : : ! Allocate matrices_iproc, that will gather the contribution of this proc to the block owned by iproc
378 : : ABI_MALLOC(ham_iproc, (cplx*buffsize_iproc(1), buffsize_iproc(2)))
379 : : ABI_MALLOC(ovl_iproc, (cplx*buffsize_iproc(1), buffsize_iproc(2)))
380 : :
381 : : ! Build them
382 : : ABI_MALLOC(left_temp, (2, npw*nspinor*buffsize_iproc(1)))
383 : : ABI_MALLOC(right_temp, (2, npw*nspinor*buffsize_iproc(2)))
384 : :
385 : : ! ovl
386 : : call from_mat_to_block_cyclic(cg, npw*nspinor, nband, left_temp, &
387 : : & buffsize_iproc(1), blocksize, coords_iproc(1), grid_dims(1))
388 : : if(usepaw == 1) then
389 : : call from_mat_to_block_cyclic(gsc, npw*nspinor, nband, right_temp, &
390 : : & buffsize_iproc(2), blocksize, coords_iproc(2), grid_dims(2))
391 : : else
392 : : call from_mat_to_block_cyclic(cg, npw*nspinor, nband, right_temp, &
393 : : & buffsize_iproc(2), blocksize, coords_iproc(2), grid_dims(2))
394 : : end if
395 : : call abi_xgemm(blas_transpose,'n',buffsize_iproc(1),buffsize_iproc(2),vectsize,cone,left_temp,vectsize,&
396 : : & right_temp,vectsize,czero,ovl_iproc,buffsize_iproc(1), x_cplx=cplx)
397 : :
398 : : ! ham
399 : : call from_mat_to_block_cyclic(ghc, npw*nspinor, nband, right_temp, &
400 : : & buffsize_iproc(2), blocksize, coords_iproc(2), grid_dims(2))
401 : : call abi_xgemm(blas_transpose,'n',buffsize_iproc(1),buffsize_iproc(2),vectsize,cone,left_temp,vectsize,&
402 : : & right_temp,vectsize,czero,ham_iproc,buffsize_iproc(1), x_cplx=cplx)
403 : :
404 : : ! Sum to iproc, and fill sca_ matrices
405 : : call xmpi_sum_master(ham_iproc, iproc, slk_communicator, ierr)
406 : : call xmpi_sum_master(ovl_iproc, iproc, slk_communicator, ierr)
407 : : if(iproc == slk_processor%my_rank) then
408 : : ! DCOPY to bypass the real/complex issue
409 : : if(cplx == 2) then
410 : : call DCOPY(cplx*buffsize_iproc(1)*buffsize_iproc(2), ham_iproc, 1, sca_ham%buffer_cplx, 1)
411 : : call DCOPY(cplx*buffsize_iproc(1)*buffsize_iproc(2), ovl_iproc, 1, sca_ovl%buffer_cplx, 1)
412 : : else
413 : : call DCOPY(cplx*buffsize_iproc(1)*buffsize_iproc(2), ham_iproc, 1, sca_ham%buffer_real, 1)
414 : : call DCOPY(cplx*buffsize_iproc(1)*buffsize_iproc(2), ovl_iproc, 1, sca_ovl%buffer_real, 1)
415 : : end if
416 : : end if
417 : :
418 : : ABI_FREE(ham_iproc)
419 : : ABI_FREE(ovl_iproc)
420 : : ABI_FREE(left_temp)
421 : : ABI_FREE(right_temp)
422 : : end do
423 : :
424 : : ! Final sum
425 : : if(cplx == 2) then
426 : : call xmpi_sum(sca_ham%buffer_cplx, slk_complement_communicator, ierr)
427 : : call xmpi_sum(sca_ovl%buffer_cplx, slk_complement_communicator, ierr)
428 : : else
429 : : call xmpi_sum(sca_ham%buffer_real, slk_complement_communicator, ierr)
430 : : call xmpi_sum(sca_ovl%buffer_real, slk_complement_communicator, ierr)
431 : : end if
432 : :
433 : : ! Transform back
434 : : if(istwf_k == 2) then
435 : : cg = cg / sqrt2
436 : : if(mpi_enreg%me_g0 == 1) cg(:, 1:npw*nspinor*nband:npw) = cg(:, 1:npw*nspinor*nband:npw) * sqrt2
437 : : ghc = ghc / sqrt2
438 : : if(mpi_enreg%me_g0 == 1) ghc(:, 1:npw*nspinor*nband:npw) = ghc(:, 1:npw*nspinor*nband:npw) * sqrt2
439 : : if(usepaw == 1) then
440 : : gsc = gsc / sqrt2
441 : : if(mpi_enreg%me_g0 == 1) gsc(:, 1:npw*nspinor*nband:npw) = gsc(:, 1:npw*nspinor*nband:npw) * sqrt2
442 : : end if
443 : : end if
444 : : call timab(timer_subham, 2, tsec)
445 : :
446 : : !======================================================================================================
447 : : ! Do the diagonalization
448 : : !======================================================================================================
449 : : !write(message, *) 'RR: diag'
450 : : !call wrtout(std_out,message)
451 : : call timab(timer_subdiago, 1, tsec)
452 : : call compute_generalized_eigen_problem(slk_processor,sca_ham,sca_ovl,&
453 : : & sca_evec,eig,slk_communicator,istwf_k)
454 : : call timab(timer_subdiago, 2, tsec)
455 : :
456 : : !======================================================================================================
457 : : ! Perform rotation
458 : : !======================================================================================================
459 : : call timab(timer_rotation, 1, tsec)
460 : : cg_new = zero
461 : : gsc_or_vnlxc_new = zero
462 : : ghc_new = zero
463 : : do iproc=0,nbproc-1
464 : : ! Compute the contribution to the rotated matrices from this block
465 : : !write(message, *) 'RR: rot', iproc
466 : : !call wrtout(std_out,message)
467 : :
468 : : ! Get coordinates of iproc
469 : : coords_iproc(1) = INT(iproc / grid_dims(2))
470 : : coords_iproc(2) = MOD(iproc, grid_dims(2))
471 : :
472 : : ! Get buffersize of iproc
473 : : buffsize_iproc(1) = NUMROC(nband,blocksize,coords_iproc(1),0,slk_processor%grid%dims(1))
474 : : buffsize_iproc(2) = NUMROC(nband,blocksize,coords_iproc(2),0,slk_processor%grid%dims(2))
475 : :
476 : : ! Get data from iproc
477 : : ABI_MALLOC(evec_iproc, (cplx*buffsize_iproc(1), buffsize_iproc(2)))
478 : : if(iproc == slk_processor%my_rank) then
479 : : if(cplx == 2) then
480 : : call DCOPY(cplx*buffsize_iproc(1)*buffsize_iproc(2), sca_evec%buffer_cplx, 1, evec_iproc, 1)
481 : : else
482 : : call DCOPY(cplx*buffsize_iproc(1)*buffsize_iproc(2), sca_evec%buffer_real, 1, evec_iproc, 1)
483 : : end if
484 : : end if
485 : : call xmpi_bcast(evec_iproc,iproc,slk_communicator,ierr)
486 : :
487 : : ! Compute contribution to the rotated matrices from iproc
488 : : ABI_MALLOC(left_temp, (2,npw*nspinor*buffsize_iproc(1)))
489 : : ABI_MALLOC(right_temp, (2,npw*nspinor*buffsize_iproc(2)))
490 : :
491 : : ! cg
492 : : call from_mat_to_block_cyclic(cg, npw*nspinor, nband, left_temp, &
493 : : & buffsize_iproc(1), blocksize, coords_iproc(1), grid_dims(1))
494 : : call abi_xgemm('n','n',vectsize,buffsize_iproc(2),buffsize_iproc(1),cone,left_temp,vectsize,&
495 : : & evec_iproc, buffsize_iproc(1), czero, right_temp, vectsize, x_cplx=cplx)
496 : : call from_block_cyclic_to_mat(cg_new, npw*nspinor, nband, right_temp, &
497 : : & buffsize_iproc(2), blocksize, coords_iproc(2), grid_dims(2))
498 : :
499 : : ! ghc
500 : : call from_mat_to_block_cyclic(ghc, npw*nspinor, nband, left_temp, &
501 : : & buffsize_iproc(1), blocksize, coords_iproc(1), grid_dims(1))
502 : : call abi_xgemm('n','n',vectsize,buffsize_iproc(2),buffsize_iproc(1),cone,left_temp,vectsize,&
503 : : & evec_iproc, buffsize_iproc(1), czero, right_temp, vectsize,x_cplx=cplx)
504 : : call from_block_cyclic_to_mat(ghc_new, npw*nspinor, nband, right_temp,&
505 : : & buffsize_iproc(2), blocksize, coords_iproc(2), grid_dims(2))
506 : :
507 : : ! gsc or vnlc
508 : : if(usepaw == 1) then
509 : : call from_mat_to_block_cyclic(gsc, npw*nspinor, nband, left_temp, &
510 : : & buffsize_iproc(1), blocksize, coords_iproc(1), grid_dims(1))
511 : : endif
512 : : if(usepaw==0 .or. has_fock)then
513 : : call from_mat_to_block_cyclic(gvnlxc, npw*nspinor, nband, left_temp, &
514 : : & buffsize_iproc(1), blocksize, coords_iproc(1), grid_dims(1))
515 : : end if
516 : : call abi_xgemm('n','n',vectsize,buffsize_iproc(2),buffsize_iproc(1),cone,left_temp,vectsize,&
517 : : & evec_iproc, buffsize_iproc(1), czero, right_temp, vectsize, x_cplx=cplx)
518 : : call from_block_cyclic_to_mat(gsc_or_vnlxc_new, npw*nspinor, nband, right_temp, &
519 : : & buffsize_iproc(2), blocksize, coords_iproc(2), grid_dims(2))
520 : :
521 : : ABI_FREE(evec_iproc)
522 : : ABI_FREE(left_temp)
523 : : ABI_FREE(right_temp)
524 : : end do
525 : :
526 : : ! Overwrite with _new
527 : : cg = cg_new
528 : : ghc = ghc_new
529 : : if(usepaw == 1) then
530 : : gsc = gsc_or_vnlxc_new
531 : : endif
532 : : if(usepaw==0 .or. has_fock)then
533 : : gvnlxc = gsc_or_vnlxc_new
534 : : end if
535 : : call timab(timer_rotation, 2, tsec)
536 : :
537 : : call sca_ham%free()
538 : : call sca_ovl%free()
539 : : call sca_evec%free()
540 : :
541 : : end subroutine rayleigh_ritz_distributed
542 : : !!***
543 : :
544 : : !!****f* ABINIT/from_mat_to_block_cyclic
545 : : !! NAME
546 : : !! from_mat_to_block_cyclic
547 : : !!
548 : : !! FUNCTION
549 : : !! Fills block_cyclic_mat with the columns of full_mat owned by iproc, using a 1D block-cyclic distribution
550 : : !!
551 : : !! INPUTS
552 : : !! full_mat(2,vectsize*nband)=full input matrix
553 : : !! vectsize=number of rows
554 : : !! nband=number of columns of the full matrix
555 : : !! buffsize=size of the local part of the matrix
556 : : !! blocksize=size of block for the block-cyclic distribution
557 : : !! iproc=local processor number
558 : : !! nprocs=total number of processors
559 : : !!
560 : : !! OUTPUT
561 : : !! block_cyclic_mat(2,vectsize*buffsize)=local part of full_mat
562 : : !!
563 : : !! SIDE EFFECTS
564 : : !!
565 : : !! SOURCE
566 : : subroutine from_mat_to_block_cyclic(full_mat, vectsize, nband, block_cyclic_mat, buffsize, blocksize, iproc, nprocs)
567 : :
568 : : integer, intent(in) :: vectsize, nband, buffsize, blocksize, iproc, nprocs
569 : : real(dp), intent(in) :: full_mat(2, vectsize*nband)
570 : : real(dp), intent(inout) :: block_cyclic_mat(2, vectsize*buffsize)
571 : :
572 : : integer :: shift_fullmat, shift_block_cyclic_mat
573 : : integer :: cur_bsize
574 : :
575 : : ! *************************************************************************
576 : :
577 : : shift_fullmat = iproc*blocksize
578 : : shift_block_cyclic_mat = 0
579 : :
580 : : do while(.true.)
581 : : cur_bsize = MIN(blocksize, nband-shift_fullmat)
582 : : if(cur_bsize <= 0) exit
583 : : block_cyclic_mat(:, shift_block_cyclic_mat*vectsize+1 : shift_block_cyclic_mat*vectsize + cur_bsize*vectsize) = &
584 : : & full_mat(:, shift_fullmat*vectsize+1 : shift_fullmat*vectsize + cur_bsize*vectsize)
585 : :
586 : : shift_block_cyclic_mat = shift_block_cyclic_mat + blocksize
587 : : shift_fullmat = shift_fullmat + blocksize*nprocs
588 : : end do
589 : :
590 : : end subroutine from_mat_to_block_cyclic
591 : : !!***
592 : :
593 : : !!****f* ABINIT/from_block_cyclic_to_mat
594 : : !! NAME
595 : : !! from_block_cyclic_to_mat
596 : : !!
597 : : !! FUNCTION
598 : : !! Fills the columns of full_mat owned by iproc with block_cyclic_mat, using a 1D block-cyclic distribution
599 : : !!
600 : : !! INPUTS
601 : : !! block_cyclic_mat(2,vectsize*buffsize)=local part of full_mat
602 : : !! vectsize=number of rows
603 : : !! nband=number of columns of the full matrix
604 : : !! buffsize=size of the local part of the matrix
605 : : !! blocksize=size of block for the block-cyclic distribution
606 : : !! iproc=local processor number
607 : : !! nprocs=total number of processors
608 : : !!
609 : : !! OUTPUT
610 : : !! full_mat(2,vectsize*nband)=full matrix
611 : : !!
612 : : !! SIDE EFFECTS
613 : : !!
614 : : !! SOURCE
615 : :
616 : : subroutine from_block_cyclic_to_mat(full_mat, vectsize, nband, block_cyclic_mat, buffsize, blocksize, iproc, nprocs)
617 : :
618 : : integer, intent(in) :: vectsize, nband, buffsize, blocksize, iproc, nprocs
619 : : real(dp), intent(inout) :: full_mat(2, vectsize*nband)
620 : : real(dp), intent(in) :: block_cyclic_mat(2, vectsize*buffsize)
621 : :
622 : : integer :: shift_fullmat, shift_block_cyclic_mat
623 : : integer :: cur_bsize
624 : :
625 : : ! *************************************************************************
626 : :
627 : : shift_fullmat = iproc*blocksize
628 : : shift_block_cyclic_mat = 0
629 : :
630 : : do while(.true.)
631 : : cur_bsize = MIN(blocksize, nband-shift_fullmat)
632 : : if(cur_bsize <= 0) exit
633 : : full_mat(:, shift_fullmat*vectsize+1 : shift_fullmat*vectsize + cur_bsize*vectsize) =&
634 : : & full_mat(:, shift_fullmat*vectsize+1 : shift_fullmat*vectsize + cur_bsize*vectsize) +&
635 : : & block_cyclic_mat(:, shift_block_cyclic_mat*vectsize+1 : shift_block_cyclic_mat*vectsize + cur_bsize*vectsize)
636 : :
637 : : shift_block_cyclic_mat = shift_block_cyclic_mat + blocksize
638 : : shift_fullmat = shift_fullmat + blocksize*nprocs
639 : : end do
640 : :
641 : : end subroutine from_block_cyclic_to_mat
642 : : !!***
643 : :
644 : : #endif
645 : : !ScaLapack
646 : :
647 : : end module m_rayleigh_ritz
648 : : !!***
|
