TABLE OF CONTENTS


ABINIT/m_gwls_Projected_BT [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_Projected_BT

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2018 ABINIT group (JLJ, BR, MC)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 
28 module m_gwls_Projected_BT
29 !----------------------------------------------------------------------------------------------------
30 ! This module contains routines to compute the projections of the Sternheimer equations for the 
31 ! BT operator, which appears in the Numerical integral.
32 !----------------------------------------------------------------------------------------------------
33 ! local modules
34 use m_gwls_utility
35 use m_gwls_wf
36 use m_gwls_TimingLog
37 use m_gwls_hamiltonian
38 use m_gwls_lineqsolver
39 use m_gwls_GWlanczos
40 use m_gwls_LanczosBasis 
41 use m_gwls_DielectricArray
42 use m_gwls_LanczosResolvents
43 ! abinit modules
44 use defs_basis
45 use m_abicore
46 use m_xmpi
47 
48 implicit none
49 save
50 private

m_hamiltonian/compute_projected_BT_shift_Lanczos [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  compute_projected_BT_shift_Lanczos

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

      cleanup_lanczosresolvents
      compute_resolvent_column_shift_lanczos_right_vectors
      invert_general_matrix,setup_lanczosresolvents,wf_block_distribute
      xmpi_allgather,xmpi_sum,zgemm,zgemv

SOURCE

102 subroutine compute_projected_BT_shift_Lanczos(nfreq, list_external_omega, lmax, modified_Lbasis,         &
103 kmax_numeric, npt_gauss, dielectric_array, array_integrand )
104 !----------------------------------------------------------------------------------------------------
105 !
106 ! This function returns the integrand
107 !
108 !         I(w', w) =  sum_{l1,l2} DielectricArray_{l1,l2}(w') * B_{l1,l2}(w',w)
109 !
110 ! where BT is obtained by Shift Lanczos for all frequencies.  It is assumed that modified_Lbasis
111 ! already contains the properly modified basis vectors.
112 !
113 !----------------------------------------------------------------------------------------------------
114 
115 !This section has been created automatically by the script Abilint (TD).
116 !Do not modify the following lines by hand.
117 #undef ABI_FUNC
118 #define ABI_FUNC 'compute_projected_BT_shift_Lanczos'
119 !End of the abilint section
120 
121 implicit none
122 
123 integer,      intent(in) :: nfreq
124 real(dp),     intent(in) :: list_external_omega(nfreq)
125 integer,      intent(in) :: lmax
126 integer,      intent(in) :: kmax_numeric
127 integer,      intent(in) :: npt_gauss
128 complex(dpc), intent(in) :: modified_Lbasis(npw_k,lmax)
129 complex(dpc), intent(in) :: dielectric_array(lmax,lmax,npt_gauss+1)
130 complex(dpc), intent(out):: array_integrand(npt_gauss+1,nfreq)
131 
132 
133 ! local variables
134 
135 logical   :: prec
136 
137 integer   :: l, l1, iw_ext, iw_prime, iw, mb, iblk, nbdblock_lanczos
138 integer   :: ierr
139 
140 integer      :: mpi_band_rank 
141 
142 integer      :: k
143 integer      :: iz, nz
144 complex(dpc) :: z
145 
146 complex(dpc), allocatable :: list_z(:)
147 
148 real(dp)     :: external_omega, omega_prime
149 
150 
151 complex(dpc), allocatable :: matrix_elements_resolvent(:,:)
152 
153 real(dp),     allocatable :: psik_wrk(:,:)
154 real(dp),     allocatable :: psikb_wrk(:,:)
155 real(dp),     allocatable :: psikg_wrk(:,:)
156 
157 complex(dpc), allocatable :: seed_vector(:)
158 
159 complex(dpc), allocatable :: right_vec_FFT(:)
160 
161 complex(dpc), allocatable :: right_vec_LA(:,:)
162 complex(dpc), allocatable :: LR_M_matrix_LA(:,:,:)
163 complex(dpc), allocatable :: Hamiltonian_Qk_LA(:,:,:) 
164 complex(dpc), allocatable :: left_vecs_LA(:,:) 
165 complex(dpc), allocatable :: shift_lanczos_matrix(:,:) 
166 
167 complex(dpc), allocatable :: work_vec(:) 
168 
169 
170 real(dp),     allocatable :: real_wrk_vec(:),   imag_wrk_vec(:)
171 real(dp),     allocatable :: real_wrk_mat(:,:), imag_wrk_mat(:,:)
172 
173 ! *************************************************************************
174 
175 ! prepare the complex frequency array
176 nz = 2*nfreq*npt_gauss
177 ABI_ALLOCATE(list_z,(nz))
178 
179 ABI_ALLOCATE(matrix_elements_resolvent, (lmax,nz))
180 
181 ABI_ALLOCATE(psik_wrk,  (2,npw_k))
182 ABI_ALLOCATE(psikb_wrk, (2,npw_kb))
183 ABI_ALLOCATE(psikg_wrk, (2,npw_g))
184 ABI_ALLOCATE(seed_vector, (npw_g))
185 
186 
187 ABI_ALLOCATE(real_wrk_vec, (kmax_numeric*blocksize))
188 ABI_ALLOCATE(imag_wrk_vec, (kmax_numeric*blocksize))
189 ABI_ALLOCATE(real_wrk_mat, (kmax_numeric,kmax_numeric*blocksize))
190 ABI_ALLOCATE(imag_wrk_mat, (kmax_numeric,kmax_numeric*blocksize))
191 
192 ABI_ALLOCATE(shift_lanczos_matrix, (kmax_numeric,kmax_numeric))
193 
194 ABI_ALLOCATE(work_vec, (kmax_numeric))
195 
196 ABI_ALLOCATE(right_vec_FFT,(kmax_numeric) )
197 ABI_ALLOCATE(right_vec_LA,(kmax_numeric, blocksize) )
198 ABI_ALLOCATE(LR_M_matrix_LA,(kmax_numeric,kmax_numeric, blocksize) )
199 ABI_ALLOCATE(Hamiltonian_Qk_LA,(npw_k, kmax_numeric, blocksize) )
200 
201 
202 ABI_ALLOCATE(left_vecs_LA,(kmax_numeric, lmax) )
203 
204 iw = 0
205 do iw_ext = 1, nfreq
206 
207 external_omega = list_external_omega(iw_ext)
208 
209 do iw_prime = 1, npt_gauss
210 
211 ! Remember! the first element of the list_omega array, which contain the imaginary 
212 ! integration frequencies, is zero (which need not be computed explicitly)
213 
214 omega_prime = list_omega(iw_prime+1)
215 
216 iw = iw+1
217 list_z(iw) = cmplx_1*external_omega-cmplx_i*omega_prime
218 
219 iw = iw+1
220 list_z(iw) = cmplx_1*external_omega+cmplx_i*omega_prime
221 
222 end do
223 
224 end do 
225 
226 
227 ! initialize the array to zero
228 array_integrand(:,:) = cmplx_0
229 
230 
231 ! prepare the shift lanczos scheme
232 prec = .false.  ! let's not precondition for now
233 call setup_LanczosResolvents(kmax_numeric,prec)
234 
235 ! Number of blocks of lanczos vectors
236 nbdblock_lanczos = lmax/blocksize
237 if (modulo(lmax,blocksize) /= 0) nbdblock_lanczos = nbdblock_lanczos + 1
238 
239 mpi_band_rank = mpi_enreg%me_band
240 
241 
242 
243 !-------------------------------------------------------------------
244 !
245 ! The shift lanczos scheme will be implemented explicitly in the
246 ! loop below instead of being wrapped in routines in the module
247 ! gwls_LanczosResolvents. The task to be performed is subtle 
248 ! because of the different data distributions and the need to 
249 ! project on all Lanczos vectors.
250 !
251 !-------------------------------------------------------------------
252 
253 ! loop on all blocks of Lanczos vectors
254 do iblk = 1, nbdblock_lanczos
255 
256 ! Convert a block of Lanczos vectors to the FFT configuration
257 
258 ! Change the configuration of the data
259 do mb =1, blocksize
260 l = (iblk-1)*blocksize+mb
261 if (l <= lmax) then 
262   psik_wrk(1,:)  = dble ( modified_Lbasis(:,l) )
263   psik_wrk(2,:)  = dimag( modified_Lbasis(:,l) )
264 else
265   psik_wrk(:,:)  = zero
266 end if
267 
268 psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:) 
269 
270 end do ! mb 
271 
272 ! change configuration of the data, from LA to FFT
273 call wf_block_distribute(psikb_wrk,  psikg_wrk, 1) ! LA -> FFT
274 
275 ! compute the seed vector for this block
276 seed_vector(:) = cmplx_1*psikg_wrk(1,:) + cmplx_i*psikg_wrk(2,:)
277 
278 
279 ! Compute the Lanczos basis for the Hamiltonian, in FFT configuration,
280 ! given the seed vector. Project the right vector onto the thus produced Lanczos basis.
281 call compute_resolvent_column_shift_lanczos_right_vectors(seed_vector, right_vec_FFT)
282 
283 ! Distribute the data from the FFT configuration back to the LA configuration
284 ! We will use some public data from the LanczosResolvent module.
285 !
286 ! PATTERN: xmpi_allgather(xval,nelem,recvbuf,spaceComm,ier) 
287 ! unfortunately, there are only interfaces for real arguments, not complex.
288 
289 call xmpi_allgather( dble(right_vec_FFT), kmax_numeric, real_wrk_vec, mpi_enreg%comm_band, ierr) 
290 call xmpi_allgather(dimag(right_vec_FFT), kmax_numeric, imag_wrk_vec, mpi_enreg%comm_band, ierr) 
291 
292 
293 
294 call xmpi_allgather( dble(LR_M_matrix), kmax_numeric**2, real_wrk_mat, mpi_enreg%comm_band, ierr) 
295 call xmpi_allgather(dimag(LR_M_matrix), kmax_numeric**2, imag_wrk_mat, mpi_enreg%comm_band, ierr) 
296 
297 
298 do mb =1, blocksize
299 right_vec_LA    (:,mb) = cmplx_1*real_wrk_vec((mb-1)*kmax_numeric+1:mb*kmax_numeric) + &
300 &                                    cmplx_i*imag_wrk_vec((mb-1)*kmax_numeric+1:mb*kmax_numeric)
301 
302 
303 LR_M_matrix_LA(:,:,mb) = cmplx_1*real_wrk_mat(:,(mb-1)*kmax_numeric+1:mb*kmax_numeric) + &
304 &                                    cmplx_i*imag_wrk_mat(:,(mb-1)*kmax_numeric+1:mb*kmax_numeric)
305 end do ! mb
306 
307 do k = 1, kmax_numeric
308 
309 psikg_wrk(1,:) = dble ( Hamiltonian_Qk(:,k) )
310 psikg_wrk(2,:) = dimag( Hamiltonian_Qk(:,k) )
311 
312 ! change configuration of the data, from FFT to LA
313 call wf_block_distribute(psikb_wrk,  psikg_wrk, 2) ! FFT -> LA 
314 
315 do mb=1, blocksize
316 
317 Hamiltonian_Qk_LA(:, k, mb) =  cmplx_1*psikb_wrk(1,(mb-1)*npw_k+1:mb*npw_k) &
318 &                                              +cmplx_i*psikb_wrk(2,(mb-1)*npw_k+1:mb*npw_k) 
319 
320 end do ! mb
321 
322 end do ! k
323 
324 
325 ! Data is now available in LA configuration, perfect for linear algebra!
326 
327 
328 do mb = 1, blocksize
329 ! loop on all vectors in this block
330 
331 l1 = (iblk-1)*blocksize+mb
332 
333 if (l1 > lmax) cycle
334 
335 !  Compute Q^dagger | left_vectors > 
336 
337 ! computes C = alpha *  op(A).op(B) + beta * C
338 call ZGEMM(                         'C', &  ! First array is hermitian conjugated 
339 'N', &  ! second array is taken as is
340 kmax_numeric, &  ! number of rows of matrix op(A)
341 lmax, &  ! number of columns of matrix op(B)
342 npw_k, &  ! number of columns of op(A) == number of rows of matrix op(B)
343 cmplx_1, &  ! alpha
344 Hamiltonian_Qk_LA(:, :, mb), &  ! A matrix
345 npw_k, &  ! LDA
346 modified_Lbasis, &  ! B matrix
347 npw_k, &  ! LDB
348 cmplx_0, &  ! beta
349 left_vecs_LA, &  ! C matrix
350 kmax_numeric)  ! LDC
351 
352 call xmpi_sum(left_vecs_LA,mpi_enreg%comm_bandfft,ierr) ! sum on all processors working on LA
353 
354 ! Use shift Lanczos to compute all matrix elements! 
355 
356 ! THE FOLLOWING COULD BE DONE IN PARALLEL INSTEAD OF HAVING EVERY PROCESSOR DUMBLY DO THE SAME THING
357 do iz = 1, nz
358 
359 z = list_z(iz)
360 
361 ! Generate the matrix to be inverted
362 shift_lanczos_matrix(:,:) = LR_M_matrix_LA(:,:,mb)
363 
364 do k = 1, kmax_numeric
365 shift_lanczos_matrix(k,k) = shift_lanczos_matrix(k,k)-z
366 end do 
367 
368 ! since z could be complex, the matrix is not necessarily hermitian. Invert using general Lapack scheme
369 call invert_general_matrix(kmax_numeric,shift_lanczos_matrix)
370 ! the matrix now contains the inverse!
371 
372 
373 ! | work_vec > = M^{-1} . | right_vec >
374 
375 ! compute y = alpha op(A).x + beta y
376 
377 call ZGEMV(                       'N', &! A matrix is as is
378 kmax_numeric, &! number of rows of matrix A
379 kmax_numeric, &! number of columns of matrix A
380 cmplx_1, &! alpha
381 shift_lanczos_matrix, &! matrix A
382 kmax_numeric, &! LDA
383 right_vec_LA(:,mb), &! array X
384 1, &! INC X
385 cmplx_0, &! beta
386 work_vec, &! Y array
387 1)  ! INC Y
388 
389 
390 ! matrix_elements = < right_vecs | work_vec > 
391 
392 call ZGEMV(                       'C', &! A matrix is hermitan conjugate
393 kmax_numeric, &! number of rows of matrix A
394 lmax, &! number of columns of matrix A
395 cmplx_1, &! alpha
396 left_vecs_LA, &! matrix A
397 kmax_numeric, &! LDA
398 work_vec, &! array X
399 1, &! INC X
400 cmplx_0, &! beta
401 matrix_elements_resolvent(:,iz), &! Y array
402 1)  ! INC Y
403 
404 
405 
406 end do ! iz
407 
408 ! update integrand
409 iw = 0
410 do iw_ext = 1, nfreq
411 do iw_prime = 1, npt_gauss
412 
413 iw = iw+1
414 
415 ! this expression will be normalized by 1/2pi at the end
416 array_integrand(iw_prime+1,iw_ext) = array_integrand(iw_prime+1,iw_ext)  +    &
417 sum(dielectric_array(:,l1,iw_prime+1)*  &
418 (matrix_elements_resolvent(:,iw)+matrix_elements_resolvent(:,iw+1)))
419 
420 iw = iw+1
421 
422 end do !iw_prime 
423 end do !iw_ext 
424 
425 
426 
427 end do !mb
428 end do ! iblk
429 
430 ! normalize !
431 array_integrand(:,:)  = array_integrand(:,:)/(2.0_dp*pi)
432 
433 
434 
435 call cleanup_LanczosResolvents
436 
437 
438 ABI_DEALLOCATE(psik_wrk)
439 ABI_DEALLOCATE(psikb_wrk)
440 ABI_DEALLOCATE(psikg_wrk)
441 
442 ABI_DEALLOCATE(seed_vector)
443 ABI_DEALLOCATE(work_vec)
444 
445 
446 ABI_DEALLOCATE(right_vec_FFT)
447 ABI_DEALLOCATE(right_vec_LA)
448 
449 ABI_DEALLOCATE(LR_M_matrix_LA)
450 
451 ABI_DEALLOCATE(Hamiltonian_Qk_LA)
452 
453 ABI_DEALLOCATE(real_wrk_vec)
454 ABI_DEALLOCATE(imag_wrk_vec)
455 ABI_DEALLOCATE(real_wrk_mat)
456 ABI_DEALLOCATE(imag_wrk_mat)
457 
458 
459 ABI_DEALLOCATE(shift_lanczos_matrix)
460 ABI_DEALLOCATE(left_vecs_LA)
461 
462 
463 ABI_DEALLOCATE(list_z)
464 ABI_DEALLOCATE(matrix_elements_resolvent)
465 
466 end subroutine compute_projected_BT_shift_Lanczos

m_hamiltonian/compute_projected_BT_shift_Lanczos_DISTRIBUTED [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  compute_projected_BT_shift_Lanczos_DISTRIBUTED

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

      cleanup_lanczosresolvents
      compute_resolvent_column_shift_lanczos_right_vectors
      invert_general_matrix,setup_lanczosresolvents,wf_block_distribute
      xmpi_allgather,xmpi_sum,zgemm,zgemv

SOURCE

493 subroutine compute_projected_BT_shift_Lanczos_DISTRIBUTED(nfreq, list_external_omega, lmax,blocksize_eps,       &
494 model_lanczos_vector_belongs_to_this_node, model_lanczos_vector_index,  &
495 modified_Lbasis, kmax_numeric, npt_gauss, dielectric_array, array_integrand )
496 !----------------------------------------------------------------------------------------------------
497 !
498 ! This function returns the integrand
499 !
500 !         I(w', w) =  sum_{l1,l2} DielectricArray_{l1,l2}(w') * B_{l1,l2}(w',w)
501 !
502 ! where BT is obtained by shift Lanczos for all frequencies.  It is assumed that modified_Lbasis
503 ! already contains the properly modified basis vectors.
504 !
505 ! This routine performs the same function as compute_projected_BT_shift_Lanczos, but 
506 ! takes into account that the MODEL dielectric array is distributed over the processors.
507 !
508 !----------------------------------------------------------------------------------------------------
509 
510 !This section has been created automatically by the script Abilint (TD).
511 !Do not modify the following lines by hand.
512 #undef ABI_FUNC
513 #define ABI_FUNC 'compute_projected_BT_shift_Lanczos_DISTRIBUTED'
514 !End of the abilint section
515 
516 implicit none
517 
518 integer,      intent(in) :: nfreq
519 real(dp),     intent(in) :: list_external_omega(nfreq)
520 integer,      intent(in) :: lmax, blocksize_eps
521 integer,      intent(in) :: kmax_numeric
522 integer,      intent(in) :: npt_gauss
523 complex(dpc), intent(in) :: modified_Lbasis(npw_k,lmax)
524 complex(dpc), intent(in) :: dielectric_array(lmax, blocksize_eps, npt_gauss+1)
525 
526 logical, intent(in) :: model_lanczos_vector_belongs_to_this_node(lmax)
527 integer, intent(in) :: model_lanczos_vector_index(lmax)
528 
529 complex(dpc), intent(out):: array_integrand(npt_gauss+1,nfreq)
530 
531 
532 
533 ! local variables
534 
535 logical   :: prec
536 
537 integer   :: l, l1, lb, iw_ext, iw_prime, iw, mb, iblk, nbdblock_lanczos
538 integer   :: ierr
539 
540 integer      :: mpi_band_rank 
541 
542 integer      :: k
543 integer      :: iz, nz
544 complex(dpc) :: z
545 
546 complex(dpc), allocatable :: list_z(:)
547 
548 real(dp)     :: external_omega, omega_prime
549 
550 
551 complex(dpc), allocatable :: matrix_elements_resolvent(:,:)
552 
553 real(dp),     allocatable :: psik_wrk(:,:)
554 real(dp),     allocatable :: psikb_wrk(:,:)
555 real(dp),     allocatable :: psikg_wrk(:,:)
556 
557 complex(dpc), allocatable :: seed_vector(:)
558 
559 complex(dpc), allocatable :: right_vec_FFT(:)
560 
561 complex(dpc), allocatable :: right_vec_LA(:,:)
562 complex(dpc), allocatable :: LR_M_matrix_LA(:,:,:)
563 complex(dpc), allocatable :: Hamiltonian_Qk_LA(:,:,:) 
564 complex(dpc), allocatable :: left_vecs_LA(:,:) 
565 complex(dpc), allocatable :: shift_lanczos_matrix(:,:) 
566 
567 complex(dpc), allocatable :: work_vec(:) 
568 
569 
570 real(dp),     allocatable :: real_wrk_vec(:),   imag_wrk_vec(:)
571 real(dp),     allocatable :: real_wrk_mat(:,:), imag_wrk_mat(:,:)
572 
573 ! *************************************************************************
574 
575 ! prepare the complex frequency array
576 nz = 2*nfreq*npt_gauss
577 ABI_ALLOCATE(list_z,(nz))
578 
579 ABI_ALLOCATE(matrix_elements_resolvent, (lmax,nz))
580 
581 ABI_ALLOCATE(psik_wrk,  (2,npw_k))
582 ABI_ALLOCATE(psikb_wrk, (2,npw_kb))
583 ABI_ALLOCATE(psikg_wrk, (2,npw_g))
584 ABI_ALLOCATE(seed_vector, (npw_g))
585 
586 
587 ABI_ALLOCATE(real_wrk_vec, (kmax_numeric*blocksize))
588 ABI_ALLOCATE(imag_wrk_vec, (kmax_numeric*blocksize))
589 ABI_ALLOCATE(real_wrk_mat, (kmax_numeric,kmax_numeric*blocksize))
590 ABI_ALLOCATE(imag_wrk_mat, (kmax_numeric,kmax_numeric*blocksize))
591 
592 ABI_ALLOCATE(shift_lanczos_matrix, (kmax_numeric,kmax_numeric))
593 
594 ABI_ALLOCATE(work_vec, (kmax_numeric))
595 
596 ABI_ALLOCATE(right_vec_FFT,(kmax_numeric) )
597 ABI_ALLOCATE(right_vec_LA,(kmax_numeric, blocksize) )
598 ABI_ALLOCATE(LR_M_matrix_LA,(kmax_numeric,kmax_numeric, blocksize) )
599 ABI_ALLOCATE(Hamiltonian_Qk_LA,(npw_k, kmax_numeric, blocksize) )
600 
601 
602 ABI_ALLOCATE(left_vecs_LA,(kmax_numeric, lmax) )
603 
604 iw = 0
605 do iw_ext = 1, nfreq
606 
607 external_omega = list_external_omega(iw_ext)
608 
609 do iw_prime = 1, npt_gauss
610 
611 ! Remember! the first element of the list_omega array, which contain the imaginary 
612 ! integration frequencies, is zero (which need not be computed explicitly)
613 
614 omega_prime = list_omega(iw_prime+1)
615 
616 iw = iw+1
617 list_z(iw) = cmplx_1*external_omega-cmplx_i*omega_prime
618 
619 iw = iw+1
620 list_z(iw) = cmplx_1*external_omega+cmplx_i*omega_prime
621 
622 end do
623 
624 end do 
625 
626 
627 ! initialize the array to zero
628 array_integrand(:,:) = cmplx_0
629 
630 
631 ! prepare the shift lanczos scheme
632 prec = .false.  ! let's not precondition for now
633 call setup_LanczosResolvents(kmax_numeric,prec)
634 
635 ! Number of blocks of lanczos vectors
636 nbdblock_lanczos = lmax/blocksize
637 if (modulo(lmax,blocksize) /= 0) nbdblock_lanczos = nbdblock_lanczos + 1
638 
639 mpi_band_rank = mpi_enreg%me_band
640 
641 
642 
643 !-------------------------------------------------------------------
644 !
645 ! The shift lanczos scheme will be implemented explicitly in the
646 ! loop below instead of being wrapped in routines in the module
647 ! gwls_LanczosResolvents. The task to be performed is subtle 
648 ! because of the different data distributions and the need to 
649 ! project on all Lanczos vectors.
650 !
651 !-------------------------------------------------------------------
652 
653 ! loop on all blocks of Lanczos vectors
654 do iblk = 1, nbdblock_lanczos
655 
656 ! Convert a block of Lanczos vectors to the FFT configuration
657 
658 ! Change the configuration of the data
659 do mb =1, blocksize
660 l = (iblk-1)*blocksize+mb
661 if (l <= lmax) then 
662   psik_wrk(1,:)  = dble ( modified_Lbasis(:,l) )
663   psik_wrk(2,:)  = dimag( modified_Lbasis(:,l) )
664 else
665   psik_wrk(:,:)  = zero
666 end if
667 
668 psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:) 
669 
670 end do ! mb 
671 
672 ! change configuration of the data, from LA to FFT
673 call wf_block_distribute(psikb_wrk,  psikg_wrk, 1) ! LA -> FFT
674 
675 ! compute the seed vector for this block
676 seed_vector(:) = cmplx_1*psikg_wrk(1,:) + cmplx_i*psikg_wrk(2,:)
677 
678 
679 ! Compute the Lanczos basis for the Hamiltonian, in FFT configuration,
680 ! given the seed vector. Project the right vector onto the thus produced Lanczos basis.
681 call compute_resolvent_column_shift_lanczos_right_vectors(seed_vector, right_vec_FFT)
682 
683 ! Distribute the data from the FFT configuration back to the LA configuration
684 ! We will use some public data from the LanczosResolvent module.
685 !
686 ! PATTERN: xmpi_allgather(xval,nelem,recvbuf,spaceComm,ier) 
687 ! unfortunately, there are only interfaces for real arguments, not complex.
688 
689 call xmpi_allgather( dble(right_vec_FFT), kmax_numeric, real_wrk_vec, mpi_enreg%comm_band, ierr) 
690 call xmpi_allgather(dimag(right_vec_FFT), kmax_numeric, imag_wrk_vec, mpi_enreg%comm_band, ierr) 
691 
692 
693 
694 call xmpi_allgather( dble(LR_M_matrix), kmax_numeric**2, real_wrk_mat, mpi_enreg%comm_band, ierr) 
695 call xmpi_allgather(dimag(LR_M_matrix), kmax_numeric**2, imag_wrk_mat, mpi_enreg%comm_band, ierr) 
696 
697 
698 do mb =1, blocksize
699 right_vec_LA    (:,mb) = cmplx_1*real_wrk_vec((mb-1)*kmax_numeric+1:mb*kmax_numeric) + &
700 &                                    cmplx_i*imag_wrk_vec((mb-1)*kmax_numeric+1:mb*kmax_numeric)
701 
702 
703 LR_M_matrix_LA(:,:,mb) = cmplx_1*real_wrk_mat(:,(mb-1)*kmax_numeric+1:mb*kmax_numeric) + &
704 &                                    cmplx_i*imag_wrk_mat(:,(mb-1)*kmax_numeric+1:mb*kmax_numeric)
705 end do ! mb
706 
707 do k = 1, kmax_numeric
708 
709 psikg_wrk(1,:) = dble ( Hamiltonian_Qk(:,k) )
710 psikg_wrk(2,:) = dimag( Hamiltonian_Qk(:,k) )
711 
712 ! change configuration of the data, from FFT to LA
713 call wf_block_distribute(psikb_wrk,  psikg_wrk, 2) ! FFT -> LA 
714 
715 do mb=1, blocksize
716 
717 Hamiltonian_Qk_LA(:, k, mb) =  cmplx_1*psikb_wrk(1,(mb-1)*npw_k+1:mb*npw_k) &
718 &                                             +cmplx_i*psikb_wrk(2,(mb-1)*npw_k+1:mb*npw_k) 
719 
720 end do ! mb
721 
722 end do ! k
723 
724 
725 ! Data is now available in LA configuration, perfect for linear algebra!
726 
727 
728 do mb = 1, blocksize
729 ! loop on all vectors in this block
730 
731 l1 = (iblk-1)*blocksize+mb
732 
733 if (l1 > lmax) cycle
734 
735 !  Compute Q^dagger | left_vectors > 
736 
737 ! computes C = alpha *  op(A).op(B) + beta * C
738 call ZGEMM(                         'C', &  ! First array is hermitian conjugated 
739 'N', &  ! second array is taken as is
740 kmax_numeric, &  ! number of rows of matrix op(A)
741 lmax, &  ! number of columns of matrix op(B)
742 npw_k, &  ! number of columns of op(A) == number of rows of matrix op(B)
743 cmplx_1, &  ! alpha
744 Hamiltonian_Qk_LA(:, :, mb), &  ! A matrix
745 npw_k, &  ! LDA
746 modified_Lbasis, &  ! B matrix
747 npw_k, &  ! LDB
748 cmplx_0, &  ! beta
749 left_vecs_LA, &  ! C matrix
750 kmax_numeric)  ! LDC
751 
752 call xmpi_sum(left_vecs_LA,mpi_enreg%comm_bandfft,ierr) ! sum on all processors working on LA
753 
754 ! Use shift Lanczos to compute all matrix elements! 
755 
756 ! THE FOLLOWING COULD BE DONE IN PARALLEL INSTEAD OF HAVING EVERY PROCESSOR DUMBLY DO THE SAME THING
757 do iz = 1, nz
758 
759 z = list_z(iz)
760 
761 ! Generate the matrix to be inverted
762 shift_lanczos_matrix(:,:) = LR_M_matrix_LA(:,:,mb)
763 
764 do k = 1, kmax_numeric
765 shift_lanczos_matrix(k,k) = shift_lanczos_matrix(k,k)-z
766 end do 
767 
768 ! since z could be complex, the matrix is not necessarily hermitian. Invert using general Lapack scheme
769 call invert_general_matrix(kmax_numeric,shift_lanczos_matrix)
770 ! the matrix now contains the inverse!
771 
772 
773 ! | work_vec > = M^{-1} . | right_vec >
774 
775 ! compute y = alpha op(A).x + beta y
776 
777 call ZGEMV(                       'N', &! A matrix is as is
778 kmax_numeric, &! number of rows of matrix A
779 kmax_numeric, &! number of columns of matrix A
780 cmplx_1, &! alpha
781 shift_lanczos_matrix, &! matrix A
782 kmax_numeric, &! LDA
783 right_vec_LA(:,mb), &! array X
784 1, &! INC X
785 cmplx_0, &! beta
786 work_vec, &! Y array
787 1)  ! INC Y
788 
789 
790 ! matrix_elements = < right_vecs | work_vec > 
791 
792 call ZGEMV(                       'C', &! A matrix is hermitan conjugate
793 kmax_numeric, &! number of rows of matrix A
794 lmax, &! number of columns of matrix A
795 cmplx_1, &! alpha
796 left_vecs_LA, &! matrix A
797 kmax_numeric, &! LDA
798 work_vec, &! array X
799 1, &! INC X
800 cmplx_0, &! beta
801 matrix_elements_resolvent(:,iz), &! Y array
802 1)  ! INC Y
803 
804 
805 
806 end do ! iz
807 
808 if ( model_lanczos_vector_belongs_to_this_node(l1) ) then
809 
810   lb = model_lanczos_vector_index(l1)
811 
812 
813 
814   ! update integrand
815   iw = 0
816   do iw_ext = 1, nfreq
817   do iw_prime = 1, npt_gauss
818 
819   iw = iw+1
820 
821   ! this expression will be normalized by 1/2pi at the end
822   array_integrand(iw_prime+1,iw_ext) = array_integrand(iw_prime+1,iw_ext)  +    &
823   sum(dielectric_array(:,lb,iw_prime+1)*  &
824   (matrix_elements_resolvent(:,iw)+matrix_elements_resolvent(:,iw+1)))
825 
826   iw = iw+1
827 
828   end do !iw_prime 
829   end do !iw_ext 
830 
831 end if
832 
833 end do !mb
834 end do ! iblk
835 
836 ! sum on processors !
837 
838 call xmpi_sum(array_integrand, mpi_enreg%comm_bandfft, ierr) ! sum on all processors for LA configuration
839 
840 ! normalize !
841 array_integrand(:,:)  = array_integrand(:,:)/(2.0_dp*pi)
842 
843 
844 
845 call cleanup_LanczosResolvents
846 
847 
848 ABI_DEALLOCATE(psik_wrk)
849 ABI_DEALLOCATE(psikb_wrk)
850 ABI_DEALLOCATE(psikg_wrk)
851 
852 ABI_DEALLOCATE(seed_vector)
853 ABI_DEALLOCATE(work_vec)
854 
855 
856 ABI_DEALLOCATE(right_vec_FFT)
857 ABI_DEALLOCATE(right_vec_LA)
858 
859 ABI_DEALLOCATE(LR_M_matrix_LA)
860 
861 ABI_DEALLOCATE(Hamiltonian_Qk_LA)
862 
863 ABI_DEALLOCATE(real_wrk_vec)
864 ABI_DEALLOCATE(imag_wrk_vec)
865 ABI_DEALLOCATE(real_wrk_mat)
866 ABI_DEALLOCATE(imag_wrk_mat)
867 
868 
869 ABI_DEALLOCATE(shift_lanczos_matrix)
870 ABI_DEALLOCATE(left_vecs_LA)
871 
872 
873 ABI_DEALLOCATE(list_z)
874 ABI_DEALLOCATE(matrix_elements_resolvent)
875 
876 end subroutine compute_projected_BT_shift_Lanczos_DISTRIBUTED