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