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-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