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