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