Line data Source code
1 : !!****m* ABINIT/m_gwls_Projected_BT
2 : !! NAME
3 : !! m_gwls_Projected_BT
4 : !!
5 : !! FUNCTION
6 : !! .
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2009-2026 ABINIT group (JLJ, BR, MC)
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 :
23 : module m_gwls_Projected_BT
24 : !----------------------------------------------------------------------------------------------------
25 : ! This module contains routines to compute the projections of the Sternheimer equations for the
26 : ! BT operator, which appears in the Numerical integral.
27 : !----------------------------------------------------------------------------------------------------
28 : ! local modules
29 : use m_gwls_utility
30 : use m_gwls_wf
31 : use m_gwls_TimingLog
32 : use m_gwls_hamiltonian
33 : use m_gwls_lineqsolver
34 : use m_gwls_GWlanczos
35 : use m_gwls_LanczosBasis
36 : use m_gwls_DielectricArray
37 : use m_gwls_LanczosResolvents
38 : ! abinit modules
39 : use defs_basis
40 : use m_abicore
41 : use m_xmpi
42 :
43 : implicit none
44 : save
45 : private
46 : !!***
47 :
48 : ! Module arrays, to be set and used in this module
49 : integer, public :: BT_lsternheimer
50 : complex(dp), public, allocatable :: projected_BT_A_matrix(:,:)
51 : complex(dp), public, allocatable :: projected_BT_L_matrix(:,:)
52 : complex(dp), public, allocatable :: projected_BT_BETA_matrix(:,:)
53 :
54 : integer, public :: BT_lsternheimer_Lanczos
55 : ! We use pointers so that arrays may be allocated in a subroutine
56 : complex(dp), public, pointer:: projected_BT_A_matrix_Lanczos(:,:)
57 : complex(dp), public, pointer:: projected_BT_L_matrix_Lanczos(:,:)
58 : complex(dp), public, pointer:: projected_BT_BETA_matrix_Lanczos(:,:)
59 :
60 :
61 : integer, public :: BT_lsternheimer_model_Lanczos
62 : ! We use pointers so that arrays may be allocated in a subroutine
63 : complex(dp), public, pointer :: projected_BT_A_matrix_model_Lanczos(:,:)
64 : complex(dp), public, pointer :: projected_BT_L_matrix_model_Lanczos(:,:)
65 : complex(dp), public, pointer :: projected_BT_BETA_matrix_model_Lanczos(:,:)
66 : !!***
67 :
68 :
69 : ! Public methods
70 : public :: compute_projected_BT_shift_Lanczos
71 : public :: compute_projected_BT_shift_Lanczos_DISTRIBUTED
72 : !!***
73 : contains
74 :
75 : !!****f* m_hamiltonian/compute_projected_BT_shift_Lanczos
76 : !! NAME
77 : !! compute_projected_BT_shift_Lanczos
78 : !!
79 : !! FUNCTION
80 : !! .
81 : !!
82 : !! INPUTS
83 : !!
84 : !! OUTPUT
85 : !!
86 : !! SOURCE
87 :
88 7 : subroutine compute_projected_BT_shift_Lanczos(nfreq, list_external_omega, lmax, modified_Lbasis, &
89 7 : kmax_numeric, npt_gauss, dielectric_array, array_integrand )
90 : !----------------------------------------------------------------------------------------------------
91 : !
92 : ! This function returns the integrand
93 : !
94 : ! I(w', w) = sum_{l1,l2} DielectricArray_{l1,l2}(w') * B_{l1,l2}(w',w)
95 : !
96 : ! where BT is obtained by Shift Lanczos for all frequencies. It is assumed that modified_Lbasis
97 : ! already contains the properly modified basis vectors.
98 : !
99 : !----------------------------------------------------------------------------------------------------
100 : integer, intent(in) :: nfreq
101 : real(dp), intent(in) :: list_external_omega(nfreq)
102 : integer, intent(in) :: lmax
103 : integer, intent(in) :: kmax_numeric
104 : integer, intent(in) :: npt_gauss
105 : complex(dp), intent(in) :: modified_Lbasis(npw_k,lmax)
106 : complex(dp), intent(in) :: dielectric_array(lmax,lmax,npt_gauss+1)
107 : complex(dp), intent(out):: array_integrand(npt_gauss+1,nfreq)
108 :
109 :
110 : ! local variables
111 :
112 : logical :: prec
113 :
114 : integer :: l, l1, iw_ext, iw_prime, iw, mb, iblk, nbdblock_lanczos
115 : integer :: ierr
116 :
117 : integer :: mpi_band_rank
118 :
119 : integer :: k
120 : integer :: iz, nz
121 : complex(dp) :: z
122 :
123 7 : complex(dp), allocatable :: list_z(:)
124 :
125 : real(dp) :: external_omega, omega_prime
126 :
127 :
128 7 : complex(dp), allocatable :: matrix_elements_resolvent(:,:)
129 :
130 7 : real(dp), allocatable :: psik_wrk(:,:)
131 7 : real(dp), allocatable :: psikb_wrk(:,:)
132 7 : real(dp), allocatable :: psikg_wrk(:,:)
133 :
134 7 : complex(dp), allocatable :: seed_vector(:)
135 :
136 7 : complex(dp), allocatable :: right_vec_FFT(:)
137 :
138 7 : complex(dp), allocatable :: right_vec_LA(:,:)
139 7 : complex(dp), allocatable :: LR_M_matrix_LA(:,:,:)
140 7 : complex(dp), allocatable :: Hamiltonian_Qk_LA(:,:,:)
141 : complex(dp), allocatable :: left_vecs_LA(:,:)
142 7 : complex(dp), allocatable :: shift_lanczos_matrix(:,:)
143 :
144 7 : complex(dp), allocatable :: work_vec(:)
145 :
146 :
147 : real(dp), allocatable :: real_wrk_vec(:), imag_wrk_vec(:)
148 : real(dp), allocatable :: real_wrk_mat(:,:), imag_wrk_mat(:,:)
149 :
150 : ! *************************************************************************
151 :
152 : ! prepare the complex frequency array
153 7 : nz = 2*nfreq*npt_gauss
154 21 : ABI_MALLOC(list_z,(nz))
155 :
156 28 : ABI_MALLOC(matrix_elements_resolvent, (lmax,nz))
157 :
158 21 : ABI_MALLOC(psik_wrk, (2,npw_k))
159 21 : ABI_MALLOC(psikb_wrk, (2,npw_kb))
160 21 : ABI_MALLOC(psikg_wrk, (2,npw_g))
161 21 : ABI_MALLOC(seed_vector, (npw_g))
162 :
163 :
164 21 : ABI_MALLOC(real_wrk_vec, (kmax_numeric*blocksize))
165 14 : ABI_MALLOC(imag_wrk_vec, (kmax_numeric*blocksize))
166 28 : ABI_MALLOC(real_wrk_mat, (kmax_numeric,kmax_numeric*blocksize))
167 21 : ABI_MALLOC(imag_wrk_mat, (kmax_numeric,kmax_numeric*blocksize))
168 :
169 28 : ABI_MALLOC(shift_lanczos_matrix, (kmax_numeric,kmax_numeric))
170 :
171 21 : ABI_MALLOC(work_vec, (kmax_numeric))
172 :
173 14 : ABI_MALLOC(right_vec_FFT,(kmax_numeric) )
174 28 : ABI_MALLOC(right_vec_LA,(kmax_numeric, blocksize) )
175 35 : ABI_MALLOC(LR_M_matrix_LA,(kmax_numeric,kmax_numeric, blocksize) )
176 35 : ABI_MALLOC(Hamiltonian_Qk_LA,(npw_k, kmax_numeric, blocksize) )
177 :
178 :
179 28 : ABI_MALLOC(left_vecs_LA,(kmax_numeric, lmax) )
180 :
181 7 : iw = 0
182 14 : do iw_ext = 1, nfreq
183 :
184 7 : external_omega = list_external_omega(iw_ext)
185 :
186 84 : do iw_prime = 1, npt_gauss
187 :
188 : ! Remember! the first element of the list_omega array, which contain the imaginary
189 : ! integration frequencies, is zero (which need not be computed explicitly)
190 :
191 70 : omega_prime = list_omega(iw_prime+1)
192 :
193 70 : iw = iw+1
194 70 : list_z(iw) = cmplx_1*external_omega-cmplx_i*omega_prime
195 :
196 70 : iw = iw+1
197 77 : list_z(iw) = cmplx_1*external_omega+cmplx_i*omega_prime
198 :
199 : end do
200 :
201 : end do
202 :
203 :
204 : ! initialize the array to zero
205 91 : array_integrand(:,:) = cmplx_0
206 :
207 :
208 : ! prepare the shift lanczos scheme
209 7 : prec = .false. ! let's not precondition for now
210 7 : call setup_LanczosResolvents(kmax_numeric,prec)
211 :
212 : ! Number of blocks of lanczos vectors
213 7 : nbdblock_lanczos = lmax/blocksize
214 7 : if (modulo(lmax,blocksize) /= 0) nbdblock_lanczos = nbdblock_lanczos + 1
215 :
216 7 : mpi_band_rank = mpi_enreg%me_band
217 :
218 :
219 :
220 : !-------------------------------------------------------------------
221 : !
222 : ! The shift lanczos scheme will be implemented explicitly in the
223 : ! loop below instead of being wrapped in routines in the module
224 : ! gwls_LanczosResolvents. The task to be performed is subtle
225 : ! because of the different data distributions and the need to
226 : ! project on all Lanczos vectors.
227 : !
228 : !-------------------------------------------------------------------
229 :
230 : ! loop on all blocks of Lanczos vectors
231 47 : do iblk = 1, nbdblock_lanczos
232 :
233 : ! Convert a block of Lanczos vectors to the FFT configuration
234 :
235 : ! Change the configuration of the data
236 96 : do mb =1, blocksize
237 56 : l = (iblk-1)*blocksize+mb
238 56 : if (l <= lmax) then
239 8280 : psik_wrk(1,:) = dble ( modified_Lbasis(:,l) )
240 8280 : psik_wrk(2,:) = dimag( modified_Lbasis(:,l) )
241 : else
242 0 : psik_wrk(:,:) = zero
243 : end if
244 :
245 24768 : psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:)
246 :
247 : end do ! mb
248 :
249 : ! change configuration of the data, from LA to FFT
250 40 : call wf_block_distribute(psikb_wrk, psikg_wrk, 1) ! LA -> FFT
251 :
252 : ! compute the seed vector for this block
253 8264 : seed_vector(:) = cmplx_1*psikg_wrk(1,:) + cmplx_i*psikg_wrk(2,:)
254 :
255 :
256 : ! Compute the Lanczos basis for the Hamiltonian, in FFT configuration,
257 : ! given the seed vector. Project the right vector onto the thus produced Lanczos basis.
258 40 : call compute_resolvent_column_shift_lanczos_right_vectors(seed_vector, right_vec_FFT)
259 :
260 : ! Distribute the data from the FFT configuration back to the LA configuration
261 : ! We will use some public data from the LanczosResolvent module.
262 : !
263 : ! PATTERN: xmpi_allgather(xval,nelem,recvbuf,spaceComm,ier)
264 : ! unfortunately, there are only interfaces for real arguments, not complex.
265 :
266 680 : call xmpi_allgather( dble(right_vec_FFT), kmax_numeric, real_wrk_vec, mpi_enreg%comm_band, ierr)
267 680 : call xmpi_allgather(dimag(right_vec_FFT), kmax_numeric, imag_wrk_vec, mpi_enreg%comm_band, ierr)
268 :
269 :
270 :
271 10920 : call xmpi_allgather( dble(LR_M_matrix), kmax_numeric**2, real_wrk_mat, mpi_enreg%comm_band, ierr)
272 10920 : call xmpi_allgather(dimag(LR_M_matrix), kmax_numeric**2, imag_wrk_mat, mpi_enreg%comm_band, ierr)
273 :
274 :
275 96 : do mb =1, blocksize
276 : right_vec_LA (:,mb) = cmplx_1*real_wrk_vec((mb-1)*kmax_numeric+1:mb*kmax_numeric) + &
277 952 : & cmplx_i*imag_wrk_vec((mb-1)*kmax_numeric+1:mb*kmax_numeric)
278 :
279 :
280 : LR_M_matrix_LA(:,:,mb) = cmplx_1*real_wrk_mat(:,(mb-1)*kmax_numeric+1:mb*kmax_numeric) + &
281 15328 : & cmplx_i*imag_wrk_mat(:,(mb-1)*kmax_numeric+1:mb*kmax_numeric)
282 : end do ! mb
283 :
284 680 : do k = 1, kmax_numeric
285 :
286 132224 : psikg_wrk(1,:) = dble ( Hamiltonian_Qk(:,k) )
287 132224 : psikg_wrk(2,:) = dimag( Hamiltonian_Qk(:,k) )
288 :
289 : ! change configuration of the data, from FFT to LA
290 640 : call wf_block_distribute(psikb_wrk, psikg_wrk, 2) ! FFT -> LA
291 :
292 1576 : do mb=1, blocksize
293 :
294 : Hamiltonian_Qk_LA(:, k, mb) = cmplx_1*psikb_wrk(1,(mb-1)*npw_k+1:mb*npw_k) &
295 133120 : & +cmplx_i*psikb_wrk(2,(mb-1)*npw_k+1:mb*npw_k)
296 :
297 : end do ! mb
298 :
299 : end do ! k
300 :
301 :
302 : ! Data is now available in LA configuration, perfect for linear algebra!
303 :
304 :
305 103 : do mb = 1, blocksize
306 : ! loop on all vectors in this block
307 :
308 56 : l1 = (iblk-1)*blocksize+mb
309 :
310 56 : if (l1 > lmax) cycle
311 :
312 : ! Compute Q^dagger | left_vectors >
313 :
314 : ! computes C = alpha * op(A).op(B) + beta * C
315 : call ZGEMM( 'C', & ! First array is hermitian conjugated
316 : 'N', & ! second array is taken as is
317 : kmax_numeric, & ! number of rows of matrix op(A)
318 : lmax, & ! number of columns of matrix op(B)
319 : npw_k, & ! number of columns of op(A) == number of rows of matrix op(B)
320 : cmplx_1, & ! alpha
321 : Hamiltonian_Qk_LA(:, :, mb), & ! A matrix
322 : npw_k, & ! LDA
323 : modified_Lbasis, & ! B matrix
324 : npw_k, & ! LDB
325 : cmplx_0, & ! beta
326 : left_vecs_LA, & ! C matrix
327 56 : kmax_numeric) ! LDC
328 :
329 56 : call xmpi_sum(left_vecs_LA,mpi_enreg%comm_bandfft,ierr) ! sum on all processors working on LA
330 :
331 : ! Use shift Lanczos to compute all matrix elements!
332 :
333 : ! THE FOLLOWING COULD BE DONE IN PARALLEL INSTEAD OF HAVING EVERY PROCESSOR DUMBLY DO THE SAME THING
334 1176 : do iz = 1, nz
335 :
336 1120 : z = list_z(iz)
337 :
338 : ! Generate the matrix to be inverted
339 305760 : shift_lanczos_matrix(:,:) = LR_M_matrix_LA(:,:,mb)
340 :
341 19040 : do k = 1, kmax_numeric
342 19040 : shift_lanczos_matrix(k,k) = shift_lanczos_matrix(k,k)-z
343 : end do
344 :
345 : ! since z could be complex, the matrix is not necessarily hermitian. Invert using general Lapack scheme
346 1120 : call invert_general_matrix(kmax_numeric,shift_lanczos_matrix)
347 : ! the matrix now contains the inverse!
348 :
349 :
350 : ! | work_vec > = M^{-1} . | right_vec >
351 :
352 : ! compute y = alpha op(A).x + beta y
353 :
354 : call ZGEMV( 'N', &! A matrix is as is
355 : kmax_numeric, &! number of rows of matrix A
356 : kmax_numeric, &! number of columns of matrix A
357 : cmplx_1, &! alpha
358 : shift_lanczos_matrix, &! matrix A
359 : kmax_numeric, &! LDA
360 : right_vec_LA(:,mb), &! array X
361 : 1, &! INC X
362 : cmplx_0, &! beta
363 : work_vec, &! Y array
364 1120 : 1) ! INC Y
365 :
366 :
367 : ! matrix_elements = < right_vecs | work_vec >
368 :
369 : call ZGEMV( 'C', &! A matrix is hermitan conjugate
370 : kmax_numeric, &! number of rows of matrix A
371 : lmax, &! number of columns of matrix A
372 : cmplx_1, &! alpha
373 : left_vecs_LA, &! matrix A
374 : kmax_numeric, &! LDA
375 : work_vec, &! array X
376 : 1, &! INC X
377 : cmplx_0, &! beta
378 : matrix_elements_resolvent(:,iz), &! Y array
379 1176 : 1) ! INC Y
380 :
381 :
382 :
383 : end do ! iz
384 :
385 : ! update integrand
386 : iw = 0
387 208 : do iw_ext = 1, nfreq
388 672 : do iw_prime = 1, npt_gauss
389 :
390 560 : iw = iw+1
391 :
392 : ! this expression will be normalized by 1/2pi at the end
393 : array_integrand(iw_prime+1,iw_ext) = array_integrand(iw_prime+1,iw_ext) + &
394 : sum(dielectric_array(:,l1,iw_prime+1)* &
395 5040 : (matrix_elements_resolvent(:,iw)+matrix_elements_resolvent(:,iw+1)))
396 :
397 616 : iw = iw+1
398 :
399 : end do !iw_prime
400 : end do !iw_ext
401 :
402 :
403 :
404 : end do !mb
405 : end do ! iblk
406 :
407 : ! normalize !
408 91 : array_integrand(:,:) = array_integrand(:,:)/(2.0_dp*pi)
409 :
410 :
411 :
412 7 : call cleanup_LanczosResolvents
413 :
414 :
415 7 : ABI_FREE(psik_wrk)
416 7 : ABI_FREE(psikb_wrk)
417 7 : ABI_FREE(psikg_wrk)
418 :
419 7 : ABI_FREE(seed_vector)
420 7 : ABI_FREE(work_vec)
421 :
422 :
423 7 : ABI_FREE(right_vec_FFT)
424 7 : ABI_FREE(right_vec_LA)
425 :
426 7 : ABI_FREE(LR_M_matrix_LA)
427 :
428 7 : ABI_FREE(Hamiltonian_Qk_LA)
429 :
430 7 : ABI_FREE(real_wrk_vec)
431 7 : ABI_FREE(imag_wrk_vec)
432 7 : ABI_FREE(real_wrk_mat)
433 7 : ABI_FREE(imag_wrk_mat)
434 :
435 :
436 7 : ABI_FREE(shift_lanczos_matrix)
437 7 : ABI_FREE(left_vecs_LA)
438 :
439 :
440 7 : ABI_FREE(list_z)
441 7 : ABI_FREE(matrix_elements_resolvent)
442 :
443 7 : end subroutine compute_projected_BT_shift_Lanczos
444 : !!***
445 :
446 :
447 :
448 : !!****f* m_hamiltonian/compute_projected_BT_shift_Lanczos_DISTRIBUTED
449 : !! NAME
450 : !! compute_projected_BT_shift_Lanczos_DISTRIBUTED
451 : !!
452 : !! FUNCTION
453 : !! .
454 : !!
455 : !! INPUTS
456 : !!
457 : !! OUTPUT
458 : !!
459 : !! SOURCE
460 :
461 6 : subroutine compute_projected_BT_shift_Lanczos_DISTRIBUTED(nfreq, list_external_omega, lmax,blocksize_eps, &
462 6 : model_lanczos_vector_belongs_to_this_node, model_lanczos_vector_index, &
463 6 : modified_Lbasis, kmax_numeric, npt_gauss, dielectric_array, array_integrand )
464 : !----------------------------------------------------------------------------------------------------
465 : !
466 : ! This function returns the integrand
467 : !
468 : ! I(w', w) = sum_{l1,l2} DielectricArray_{l1,l2}(w') * B_{l1,l2}(w',w)
469 : !
470 : ! where BT is obtained by shift Lanczos for all frequencies. It is assumed that modified_Lbasis
471 : ! already contains the properly modified basis vectors.
472 : !
473 : ! This routine performs the same function as compute_projected_BT_shift_Lanczos, but
474 : ! takes into account that the MODEL dielectric array is distributed over the processors.
475 : !
476 : !----------------------------------------------------------------------------------------------------
477 : integer, intent(in) :: nfreq
478 : real(dp), intent(in) :: list_external_omega(nfreq)
479 : integer, intent(in) :: lmax, blocksize_eps
480 : integer, intent(in) :: kmax_numeric
481 : integer, intent(in) :: npt_gauss
482 : complex(dp), intent(in) :: modified_Lbasis(npw_k,lmax)
483 : complex(dp), intent(in) :: dielectric_array(lmax, blocksize_eps, npt_gauss+1)
484 :
485 : logical, intent(in) :: model_lanczos_vector_belongs_to_this_node(lmax)
486 : integer, intent(in) :: model_lanczos_vector_index(lmax)
487 :
488 : complex(dp), intent(out):: array_integrand(npt_gauss+1,nfreq)
489 :
490 :
491 :
492 : ! local variables
493 :
494 : logical :: prec
495 :
496 : integer :: l, l1, lb, iw_ext, iw_prime, iw, mb, iblk, nbdblock_lanczos
497 : integer :: ierr
498 :
499 : integer :: mpi_band_rank
500 :
501 : integer :: k
502 : integer :: iz, nz
503 : complex(dp) :: z
504 :
505 6 : complex(dp), allocatable :: list_z(:)
506 :
507 : real(dp) :: external_omega, omega_prime
508 :
509 :
510 6 : complex(dp), allocatable :: matrix_elements_resolvent(:,:)
511 :
512 6 : real(dp), allocatable :: psik_wrk(:,:)
513 6 : real(dp), allocatable :: psikb_wrk(:,:)
514 6 : real(dp), allocatable :: psikg_wrk(:,:)
515 :
516 6 : complex(dp), allocatable :: seed_vector(:)
517 6 : complex(dp), allocatable :: right_vec_FFT(:)
518 6 : complex(dp), allocatable :: right_vec_LA(:,:)
519 6 : complex(dp), allocatable :: LR_M_matrix_LA(:,:,:)
520 6 : complex(dp), allocatable :: Hamiltonian_Qk_LA(:,:,:)
521 : complex(dp), allocatable :: left_vecs_LA(:,:)
522 6 : complex(dp), allocatable :: shift_lanczos_matrix(:,:)
523 6 : complex(dp), allocatable :: work_vec(:)
524 :
525 :
526 : real(dp), allocatable :: real_wrk_vec(:), imag_wrk_vec(:)
527 : real(dp), allocatable :: real_wrk_mat(:,:), imag_wrk_mat(:,:)
528 :
529 : ! *************************************************************************
530 :
531 : ! prepare the complex frequency array
532 6 : nz = 2*nfreq*npt_gauss
533 18 : ABI_MALLOC(list_z,(nz))
534 :
535 24 : ABI_MALLOC(matrix_elements_resolvent, (lmax,nz))
536 :
537 18 : ABI_MALLOC(psik_wrk, (2,npw_k))
538 18 : ABI_MALLOC(psikb_wrk, (2,npw_kb))
539 18 : ABI_MALLOC(psikg_wrk, (2,npw_g))
540 18 : ABI_MALLOC(seed_vector, (npw_g))
541 :
542 :
543 18 : ABI_MALLOC(real_wrk_vec, (kmax_numeric*blocksize))
544 12 : ABI_MALLOC(imag_wrk_vec, (kmax_numeric*blocksize))
545 24 : ABI_MALLOC(real_wrk_mat, (kmax_numeric,kmax_numeric*blocksize))
546 18 : ABI_MALLOC(imag_wrk_mat, (kmax_numeric,kmax_numeric*blocksize))
547 :
548 24 : ABI_MALLOC(shift_lanczos_matrix, (kmax_numeric,kmax_numeric))
549 :
550 18 : ABI_MALLOC(work_vec, (kmax_numeric))
551 :
552 12 : ABI_MALLOC(right_vec_FFT,(kmax_numeric) )
553 24 : ABI_MALLOC(right_vec_LA,(kmax_numeric, blocksize) )
554 30 : ABI_MALLOC(LR_M_matrix_LA,(kmax_numeric,kmax_numeric, blocksize) )
555 30 : ABI_MALLOC(Hamiltonian_Qk_LA,(npw_k, kmax_numeric, blocksize) )
556 :
557 :
558 24 : ABI_MALLOC(left_vecs_LA,(kmax_numeric, lmax) )
559 :
560 6 : iw = 0
561 12 : do iw_ext = 1, nfreq
562 :
563 6 : external_omega = list_external_omega(iw_ext)
564 :
565 72 : do iw_prime = 1, npt_gauss
566 :
567 : ! Remember! the first element of the list_omega array, which contain the imaginary
568 : ! integration frequencies, is zero (which need not be computed explicitly)
569 :
570 60 : omega_prime = list_omega(iw_prime+1)
571 :
572 60 : iw = iw+1
573 60 : list_z(iw) = cmplx_1*external_omega-cmplx_i*omega_prime
574 :
575 60 : iw = iw+1
576 66 : list_z(iw) = cmplx_1*external_omega+cmplx_i*omega_prime
577 :
578 : end do
579 :
580 : end do
581 :
582 :
583 : ! initialize the array to zero
584 78 : array_integrand(:,:) = cmplx_0
585 :
586 :
587 : ! prepare the shift lanczos scheme
588 6 : prec = .false. ! let's not precondition for now
589 6 : call setup_LanczosResolvents(kmax_numeric,prec)
590 :
591 : ! Number of blocks of lanczos vectors
592 6 : nbdblock_lanczos = lmax/blocksize
593 6 : if (modulo(lmax,blocksize) /= 0) nbdblock_lanczos = nbdblock_lanczos + 1
594 :
595 6 : mpi_band_rank = mpi_enreg%me_band
596 :
597 :
598 :
599 : !-------------------------------------------------------------------
600 : !
601 : ! The shift lanczos scheme will be implemented explicitly in the
602 : ! loop below instead of being wrapped in routines in the module
603 : ! gwls_LanczosResolvents. The task to be performed is subtle
604 : ! because of the different data distributions and the need to
605 : ! project on all Lanczos vectors.
606 : !
607 : !-------------------------------------------------------------------
608 :
609 : ! loop on all blocks of Lanczos vectors
610 38 : do iblk = 1, nbdblock_lanczos
611 :
612 : ! Convert a block of Lanczos vectors to the FFT configuration
613 :
614 : ! Change the configuration of the data
615 80 : do mb =1, blocksize
616 48 : l = (iblk-1)*blocksize+mb
617 48 : if (l <= lmax) then
618 6216 : psik_wrk(1,:) = dble ( modified_Lbasis(:,l) )
619 6216 : psik_wrk(2,:) = dimag( modified_Lbasis(:,l) )
620 : else
621 0 : psik_wrk(:,:) = zero
622 : end if
623 :
624 18584 : psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:)
625 :
626 : end do ! mb
627 :
628 : ! change configuration of the data, from LA to FFT
629 32 : call wf_block_distribute(psikb_wrk, psikg_wrk, 1) ! LA -> FFT
630 :
631 : ! compute the seed vector for this block
632 6200 : seed_vector(:) = cmplx_1*psikg_wrk(1,:) + cmplx_i*psikg_wrk(2,:)
633 :
634 :
635 : ! Compute the Lanczos basis for the Hamiltonian, in FFT configuration,
636 : ! given the seed vector. Project the right vector onto the thus produced Lanczos basis.
637 32 : call compute_resolvent_column_shift_lanczos_right_vectors(seed_vector, right_vec_FFT)
638 :
639 : ! Distribute the data from the FFT configuration back to the LA configuration
640 : ! We will use some public data from the LanczosResolvent module.
641 : !
642 : ! PATTERN: xmpi_allgather(xval,nelem,recvbuf,spaceComm,ier)
643 : ! unfortunately, there are only interfaces for real arguments, not complex.
644 :
645 544 : call xmpi_allgather( dble(right_vec_FFT), kmax_numeric, real_wrk_vec, mpi_enreg%comm_band, ierr)
646 544 : call xmpi_allgather(dimag(right_vec_FFT), kmax_numeric, imag_wrk_vec, mpi_enreg%comm_band, ierr)
647 :
648 :
649 :
650 8736 : call xmpi_allgather( dble(LR_M_matrix), kmax_numeric**2, real_wrk_mat, mpi_enreg%comm_band, ierr)
651 8736 : call xmpi_allgather(dimag(LR_M_matrix), kmax_numeric**2, imag_wrk_mat, mpi_enreg%comm_band, ierr)
652 :
653 :
654 80 : do mb =1, blocksize
655 : right_vec_LA (:,mb) = cmplx_1*real_wrk_vec((mb-1)*kmax_numeric+1:mb*kmax_numeric) + &
656 816 : & cmplx_i*imag_wrk_vec((mb-1)*kmax_numeric+1:mb*kmax_numeric)
657 :
658 :
659 : LR_M_matrix_LA(:,:,mb) = cmplx_1*real_wrk_mat(:,(mb-1)*kmax_numeric+1:mb*kmax_numeric) + &
660 13136 : & cmplx_i*imag_wrk_mat(:,(mb-1)*kmax_numeric+1:mb*kmax_numeric)
661 : end do ! mb
662 :
663 544 : do k = 1, kmax_numeric
664 :
665 99200 : psikg_wrk(1,:) = dble ( Hamiltonian_Qk(:,k) )
666 99200 : psikg_wrk(2,:) = dimag( Hamiltonian_Qk(:,k) )
667 :
668 : ! change configuration of the data, from FFT to LA
669 512 : call wf_block_distribute(psikb_wrk, psikg_wrk, 2) ! FFT -> LA
670 :
671 1312 : do mb=1, blocksize
672 :
673 : Hamiltonian_Qk_LA(:, k, mb) = cmplx_1*psikb_wrk(1,(mb-1)*npw_k+1:mb*npw_k) &
674 99968 : & +cmplx_i*psikb_wrk(2,(mb-1)*npw_k+1:mb*npw_k)
675 :
676 : end do ! mb
677 :
678 : end do ! k
679 :
680 :
681 : ! Data is now available in LA configuration, perfect for linear algebra!
682 :
683 :
684 86 : do mb = 1, blocksize
685 : ! loop on all vectors in this block
686 :
687 48 : l1 = (iblk-1)*blocksize+mb
688 :
689 48 : if (l1 > lmax) cycle
690 :
691 : ! Compute Q^dagger | left_vectors >
692 :
693 : ! computes C = alpha * op(A).op(B) + beta * C
694 : call ZGEMM( 'C', & ! First array is hermitian conjugated
695 : 'N', & ! second array is taken as is
696 : kmax_numeric, & ! number of rows of matrix op(A)
697 : lmax, & ! number of columns of matrix op(B)
698 : npw_k, & ! number of columns of op(A) == number of rows of matrix op(B)
699 : cmplx_1, & ! alpha
700 : Hamiltonian_Qk_LA(:, :, mb), & ! A matrix
701 : npw_k, & ! LDA
702 : modified_Lbasis, & ! B matrix
703 : npw_k, & ! LDB
704 : cmplx_0, & ! beta
705 : left_vecs_LA, & ! C matrix
706 48 : kmax_numeric) ! LDC
707 :
708 48 : call xmpi_sum(left_vecs_LA,mpi_enreg%comm_bandfft,ierr) ! sum on all processors working on LA
709 :
710 : ! Use shift Lanczos to compute all matrix elements!
711 :
712 : ! THE FOLLOWING COULD BE DONE IN PARALLEL INSTEAD OF HAVING EVERY PROCESSOR DUMBLY DO THE SAME THING
713 1008 : do iz = 1, nz
714 :
715 960 : z = list_z(iz)
716 :
717 : ! Generate the matrix to be inverted
718 262080 : shift_lanczos_matrix(:,:) = LR_M_matrix_LA(:,:,mb)
719 :
720 16320 : do k = 1, kmax_numeric
721 16320 : shift_lanczos_matrix(k,k) = shift_lanczos_matrix(k,k)-z
722 : end do
723 :
724 : ! since z could be complex, the matrix is not necessarily hermitian. Invert using general Lapack scheme
725 960 : call invert_general_matrix(kmax_numeric,shift_lanczos_matrix)
726 : ! the matrix now contains the inverse!
727 :
728 :
729 : ! | work_vec > = M^{-1} . | right_vec >
730 :
731 : ! compute y = alpha op(A).x + beta y
732 :
733 : call ZGEMV( 'N', &! A matrix is as is
734 : kmax_numeric, &! number of rows of matrix A
735 : kmax_numeric, &! number of columns of matrix A
736 : cmplx_1, &! alpha
737 : shift_lanczos_matrix, &! matrix A
738 : kmax_numeric, &! LDA
739 : right_vec_LA(:,mb), &! array X
740 : 1, &! INC X
741 : cmplx_0, &! beta
742 : work_vec, &! Y array
743 960 : 1) ! INC Y
744 :
745 :
746 : ! matrix_elements = < right_vecs | work_vec >
747 :
748 : call ZGEMV( 'C', &! A matrix is hermitan conjugate
749 : kmax_numeric, &! number of rows of matrix A
750 : lmax, &! number of columns of matrix A
751 : cmplx_1, &! alpha
752 : left_vecs_LA, &! matrix A
753 : kmax_numeric, &! LDA
754 : work_vec, &! array X
755 : 1, &! INC X
756 : cmplx_0, &! beta
757 : matrix_elements_resolvent(:,iz), &! Y array
758 1008 : 1) ! INC Y
759 :
760 :
761 :
762 : end do ! iz
763 :
764 128 : if ( model_lanczos_vector_belongs_to_this_node(l1) ) then
765 :
766 24 : lb = model_lanczos_vector_index(l1)
767 :
768 :
769 :
770 : ! update integrand
771 24 : iw = 0
772 48 : do iw_ext = 1, nfreq
773 288 : do iw_prime = 1, npt_gauss
774 :
775 240 : iw = iw+1
776 :
777 : ! this expression will be normalized by 1/2pi at the end
778 : array_integrand(iw_prime+1,iw_ext) = array_integrand(iw_prime+1,iw_ext) + &
779 : sum(dielectric_array(:,lb,iw_prime+1)* &
780 2160 : (matrix_elements_resolvent(:,iw)+matrix_elements_resolvent(:,iw+1)))
781 :
782 264 : iw = iw+1
783 :
784 : end do !iw_prime
785 : end do !iw_ext
786 :
787 : end if
788 :
789 : end do !mb
790 : end do ! iblk
791 :
792 : ! sum on processors !
793 :
794 6 : call xmpi_sum(array_integrand, mpi_enreg%comm_bandfft, ierr) ! sum on all processors for LA configuration
795 :
796 : ! normalize !
797 78 : array_integrand(:,:) = array_integrand(:,:)/(2.0_dp*pi)
798 :
799 :
800 :
801 6 : call cleanup_LanczosResolvents
802 :
803 :
804 6 : ABI_FREE(psik_wrk)
805 6 : ABI_FREE(psikb_wrk)
806 6 : ABI_FREE(psikg_wrk)
807 :
808 6 : ABI_FREE(seed_vector)
809 6 : ABI_FREE(work_vec)
810 :
811 :
812 6 : ABI_FREE(right_vec_FFT)
813 6 : ABI_FREE(right_vec_LA)
814 :
815 6 : ABI_FREE(LR_M_matrix_LA)
816 :
817 6 : ABI_FREE(Hamiltonian_Qk_LA)
818 :
819 6 : ABI_FREE(real_wrk_vec)
820 6 : ABI_FREE(imag_wrk_vec)
821 6 : ABI_FREE(real_wrk_mat)
822 6 : ABI_FREE(imag_wrk_mat)
823 :
824 :
825 6 : ABI_FREE(shift_lanczos_matrix)
826 6 : ABI_FREE(left_vecs_LA)
827 :
828 :
829 6 : ABI_FREE(list_z)
830 6 : ABI_FREE(matrix_elements_resolvent)
831 :
832 6 : end subroutine compute_projected_BT_shift_Lanczos_DISTRIBUTED
833 : !!***
834 :
835 : end module m_gwls_Projected_BT
836 : !!***
|
