TABLE OF CONTENTS
- ABINIT/m_gwls_GWlanczos
- m_gwls_GWlanczos/block_lanczos_algorithm
- m_gwls_GWlanczos/diagonalize_lanczos_banded
- m_gwls_GWlanczos/get_seeds
ABINIT/m_gwls_GWlanczos [ Modules ]
NAME
m_gwls_GWlanczos
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_GWlanczos 24 !---------------------------------------------------------------------------------------------------- 25 ! This module implements the Lanczos scheme to band diagonalize an implicit operator. 26 !---------------------------------------------------------------------------------------------------- 27 !local modules 28 use m_gwls_utility 29 use m_gwls_TimingLog 30 use m_gwls_wf 31 use m_gwls_hamiltonian 32 use m_gwls_lineqsolver 33 use m_gwls_polarisability 34 use m_gwls_QR_factorization 35 36 !abinit modules 37 use defs_basis 38 use defs_wvltypes 39 use m_abicore 40 use m_xmpi 41 use m_pawang 42 use m_errors 43 44 use m_io_tools, only : get_unit 45 use m_paw_dmft, only : paw_dmft_type 46 use m_ebands, only : ebands_init, ebands_free 47 48 implicit none 49 save 50 private
m_gwls_GWlanczos/block_lanczos_algorithm [ Functions ]
[ Top ] [ m_gwls_GWlanczos ] [ Functions ]
NAME
block_lanczos_algorithm
FUNCTION
.
INPUTS
OUTPUT
SOURCE
171 subroutine block_lanczos_algorithm(mpi_communicator,matrix_function,kmax,nseeds,Hsize,seeds,alpha,beta,Lbasis,X0,beta0,Qk) 172 !---------------------------------------------------------------------------------------------------- 173 ! 174 ! This subroutine implements the Block Lanczos algorithm for an arbitrary, user-supplied function which 175 ! returns the action of the implicit matrix, which is assumed to be Hermitian. 176 ! 177 ! 178 !---------------------------------------------------------------------------------------------------- 179 implicit none 180 181 !----------------------------------------- 182 ! interface with implicit matrix function 183 !----------------------------------------- 184 !******************************************************** 185 !*** NOTE: *** 186 !*** *** 187 !*** There *appears* to be a bug in makemake; *** 188 !*** when the name of the function in the interface *** 189 !*** below contains the word "implicit", makemake *** 190 !*** yields an error. I suppose makemake simply *** 191 !*** parses for the word "implicit" without regards *** 192 !*** for the fact that it may simply be part of a *** 193 !*** naive developper's chosen name for the routine. *** 194 !*** Correspondingly, I change the name to *** 195 !*** "matrix_function" to avoid problems. *** 196 !*** *** 197 !*** Bruno Rousseau *** 198 !*** 08/13/2012 *** 199 !******************************************************** 200 interface 201 subroutine matrix_function(v_out,v_in,l) 202 203 use defs_basis 204 205 integer, intent(in) :: l 206 complex(dpc), intent(out) :: v_out(l) 207 complex(dpc), intent(in) :: v_in(l) 208 209 end subroutine matrix_function 210 end interface 211 212 !------------------------------ 213 ! input/output variables 214 !------------------------------ 215 216 integer, intent(in) :: mpi_communicator 217 integer, intent(in) :: kmax ! number of Lanczos blocks 218 integer, intent(in) :: nseeds ! size of each blocks 219 integer, intent(in) :: Hsize ! size of the Hilbert space in which the matrix lives 220 221 complex(dpc), intent(inout):: seeds(Hsize,nseeds) ! seed vectors for the algorithm 222 ! overwritten by X_{k+1} on output 223 224 !logical, intent(in) :: ortho ! should the Lanczos vector be orthogonalized? 225 226 complex(dpc), intent(out) :: alpha(nseeds,nseeds,kmax) ! the alpha array from the Lanczos algorithm 227 complex(dpc), intent(out) :: beta(nseeds,nseeds,kmax) ! the beta array from the Lanczos algorithm 228 complex(dpc), intent(out) :: Lbasis(Hsize,nseeds*kmax) ! array containing the Lanczos basis 229 230 231 complex(dpc), intent(in),optional :: X0(Hsize,nseeds) 232 complex(dpc), intent(in),optional :: beta0(nseeds,nseeds) 233 complex(dpc), intent(in),optional :: Qk(:,:) ! array containing vectors to which 234 235 ! the basis must be orthonormalized 236 237 238 239 !------------------------------ 240 ! local variables 241 !------------------------------ 242 integer :: k, seed1 243 integer :: dum(2), lk 244 245 complex(dpc), allocatable :: xk(:,:), xkm1(:,:), rk(:,:) 246 247 integer :: ntime, itime 248 real(dp) :: total_time1, total_time2 249 real(dp) :: time1, time2 250 integer :: ierr 251 real(dp),allocatable :: list_time(:) 252 253 ! ************************************************************************* 254 255 call cpu_time(total_time1) 256 257 258 ntime = 7 259 ABI_MALLOC(list_time,(ntime)) 260 list_time(:) = zero 261 262 if(present(Qk)) then 263 dum = shape(Qk) 264 lk = dum(2) 265 end if 266 267 ABI_MALLOC( xk, (Hsize,nseeds)) 268 ABI_MALLOC( xkm1,(Hsize,nseeds)) 269 ABI_MALLOC( rk ,(Hsize,nseeds)) 270 271 272 alpha = cmplx_0 273 beta = cmplx_0 274 275 !------------------------------------------------ 276 ! Orthonormalize the seeds 277 !------------------------------------------------ 278 ! initialize the xk array with the seeds 279 xk(:,:) = seeds(:,:) 280 281 282 ! orthonormalize the block using the QR algorithm 283 ! xk is overwritten by Q, the array of orthonormal vectors 284 285 call extract_QR(mpi_communicator, Hsize,nseeds,xk) 286 287 !------------------------------------------------ 288 ! Loop on all blocks 289 !------------------------------------------------ 290 291 292 do k = 1, kmax 293 294 itime = 0 295 296 ! tabulate basis, computed at previous step 297 Lbasis(:,nseeds*(k-1)+1:nseeds*k) = xk(:,:) 298 299 300 itime = itime+1 301 call cpu_time(time1) 302 303 ! Initialize the residual array 304 do seed1 = 1, nseeds 305 ! If we are constructing the $\hat \epsilon(i\omega = 0)$ matrix (and the Lanczos basis at the same time), 306 ! note the index in which the Sternheimer solutions will be stored (for use in the projected Sternheimer section). 307 if(write_solution) index_solution = (k-1)*nseeds + seed1 308 309 call matrix_function(rk(:,seed1),xk(:,seed1),Hsize) 310 end do 311 312 call cpu_time(time2) 313 list_time(itime) = list_time(itime) + time2-time1 314 315 itime = itime+1 316 call cpu_time(time1) 317 ! compute the alpha array, alpha = X^d.A.X 318 319 call ZGEMM( 'C', & ! take Hermitian conjugate of first array 320 'N', & ! leave second array as is 321 nseeds, & ! the number of rows of the matrix op( A ) 322 nseeds, & ! the number of columns of the matrix op( B ) 323 Hsize, & ! the number of columns of the matrix op( A ) == rows of matrix op( B ) 324 cmplx_1, & ! alpha constant 325 xk, & ! matrix A 326 Hsize, & ! LDA 327 rk, & ! matrix B 328 Hsize, & ! LDB 329 cmplx_0, & ! beta constant 330 alpha(:,:,k), & ! matrix C 331 nseeds) ! LDC 332 call xmpi_sum(alpha(:,:,k),mpi_communicator,ierr) ! sum on all processors 333 334 335 336 call cpu_time(time2) 337 list_time(itime) = list_time(itime) + time2-time1 338 339 itime = itime+1 340 call cpu_time(time1) 341 ! update the residual array, rk = rk-X.alpha 342 call ZGEMM( 'N', & ! leave first array as is 343 'N', & ! leave second array as is 344 Hsize, & ! the number of rows of the matrix op( A ) 345 nseeds, & ! the number of columns of the matrix op( B ) 346 nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B ) 347 -cmplx_1, & ! alpha constant 348 xk, & ! matrix A 349 Hsize, & ! LDA 350 alpha(:,:,k), & ! matrix B 351 nseeds, & ! LDB 352 cmplx_1, & ! beta constant 353 rk, & ! matrix C 354 Hsize) ! LDC 355 356 call cpu_time(time2) 357 list_time(itime) = list_time(itime) + time2-time1 358 359 if (k .eq. 1 .and. present(X0) .and. present(beta0)) then 360 ! if k == 1, and X0,beta0 are present, 361 ! update the residual array, r1 = r1-X_{0}.beta^d_{0} 362 call ZGEMM( 'N', & ! leave first array as is 363 'C', & ! Hermitian conjugate the second array 364 Hsize, & ! the number of rows of the matrix op( A ) 365 nseeds, & ! the number of columns of the matrix op( B ) 366 nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B ) 367 -cmplx_1, & ! alpha constant 368 X0, & ! matrix A 369 Hsize, & ! LDA 370 beta0(:,:), & ! matrix B 371 nseeds, & ! LDB 372 cmplx_1, & ! beta constant 373 rk, & ! matrix C 374 Hsize) ! LDC 375 end if 376 377 378 itime = itime+1 379 call cpu_time(time1) 380 if (k .gt. 1) then 381 382 ! if k > 1, update the residual array, rk = rk-X_{k-1}.beta^d_{k-1} 383 call ZGEMM( 'N', & ! leave first array as is 384 'C', & ! Hermitian conjugate the second array 385 Hsize, & ! the number of rows of the matrix op( A ) 386 nseeds, & ! the number of columns of the matrix op( B ) 387 nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B ) 388 -cmplx_1, & ! alpha constant 389 xkm1, & ! matrix A 390 Hsize, & ! LDA 391 beta(:,:,k-1), & ! matrix B 392 nseeds, & ! LDB 393 cmplx_1, & ! beta constant 394 rk, & ! matrix C 395 Hsize) ! LDC 396 397 end if 398 call cpu_time(time2) 399 list_time(itime) = list_time(itime) + time2-time1 400 401 ! store xk for next iteration 402 xkm1(:,:) = xk(:,:) 403 404 405 ! Orthonormalize THE RESIDUAL to all previously calculated directions 406 407 itime = itime+1 408 call cpu_time(time1) 409 !if ( ortho .and. (dtset%gwcalctyp/=1) ) then !This is a test to obtain the CPU time taken by the orthogonalizations. 410 411 if(present(Qk)) then 412 ! Orthonormalize to all previously calculated directions, if 413 ! this is a restarted Lanczos step 414 call orthogonalize(mpi_communicator, Hsize,lk,nseeds,Qk,rk) 415 end if 416 417 call orthogonalize(mpi_communicator, Hsize,k*nseeds,nseeds,Lbasis(:,1:k*nseeds),rk) 418 419 !end if 420 call cpu_time(time2) 421 list_time(itime) = list_time(itime) + time2-time1 422 423 424 itime = itime+1 425 call cpu_time(time1) 426 427 ! perform QR decomposition to extract X_{k+1} and beta_{k} 428 call extract_QR(mpi_communicator, Hsize,nseeds,rk,beta(:,:,k)) 429 430 call cpu_time(time2) 431 list_time(itime) = list_time(itime) + time2-time1 432 ! copy the Q matrix (written on rk) in xk, which becomes X_{k+1} 433 xk(:,:) = rk(:,:) 434 435 436 end do !end loop on k 437 438 ! overwrite the seeds with the last vector block. 439 seeds(:,:) = xk(:,:) 440 441 ABI_FREE( xk ) 442 ABI_FREE( xkm1) 443 ABI_FREE( rk ) 444 call cpu_time(total_time2) 445 446 list_time(7) = total_time2-total_time1 447 448 call write_block_lanczos_timing_log(list_time,ntime) 449 450 ABI_FREE(list_time) 451 452 end subroutine block_lanczos_algorithm
m_gwls_GWlanczos/diagonalize_lanczos_banded [ Functions ]
[ Top ] [ m_gwls_GWlanczos ] [ Functions ]
NAME
diagonalize_lanczos_banded
FUNCTION
.
INPUTS
OUTPUT
SOURCE
468 subroutine diagonalize_lanczos_banded(kmax,nseeds,Hsize,alpha,beta,Lbasis,eigenvalues,debug) 469 !----------------------------------------------------------------------------------- 470 ! Given the result of the Lanczos algorithm, this subroutine diagonalize the banded 471 ! matrix as well as updates the basis. 472 !----------------------------------------------------------------------------------- 473 implicit none 474 475 integer, intent(in) :: kmax ! number of Lanczos blocks 476 integer, intent(in) :: nseeds ! size of each blocks 477 integer, intent(in) :: Hsize ! size of the Hilbert space in which the matrix lives 478 logical, intent(in) :: debug 479 480 complex(dpc), intent(in) :: alpha(nseeds,nseeds,kmax) ! the alpha array from the Lanczos algorithm 481 complex(dpc), intent(in) :: beta (nseeds,nseeds,kmax) ! the beta array from the Lanczos algorithm 482 483 complex(dpc), intent(inout) :: Lbasis(Hsize,nseeds*kmax) ! array containing the Lanczos basis 484 485 486 real(dp), intent(out) :: eigenvalues(nseeds*kmax) 487 488 489 ! local variables 490 491 integer :: kd ! number of superdiagonal above the diagonal in banded storage 492 integer :: ldab ! dimension of banded storage matrix 493 494 complex(dpc), allocatable :: band_storage_matrix(:,:) 495 complex(dpc), allocatable :: saved_band_storage_matrix(:,:) 496 497 complex(dpc), allocatable :: eigenvectors(:,:) 498 499 complex(dpc), allocatable :: Lbasis_tmp(:,:) 500 501 integer :: i, j 502 integer :: k 503 integer :: s1, s2 504 integer :: info 505 506 507 complex(dpc), allocatable :: work(:) 508 real(dp), allocatable :: rwork(:) 509 510 integer :: io_unit 511 character(128) :: filename 512 logical :: file_exists 513 514 integer :: debug_unit 515 character(50) :: debug_filename 516 517 ! ************************************************************************* 518 519 520 521 ! number of superdiagonals 522 kd = nseeds 523 ldab = kd + 1 524 525 ABI_MALLOC( band_storage_matrix, (ldab,nseeds*kmax)) 526 ABI_MALLOC(saved_band_storage_matrix, (ldab,nseeds*kmax)) 527 !--------------------------------------------------------- 528 ! Store banded matrix in banded format 529 !--------------------------------------------------------- 530 ! for UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 531 532 band_storage_matrix(:,:) = cmplx_0 533 534 !----------------------------------------- 535 ! Store the alpha and beta matrices 536 !----------------------------------------- 537 538 ! loop on all blocks 539 do k=1,kmax 540 541 ! alpha blocks 542 do s2 = 1, nseeds 543 j = (k-1)*nseeds+s2 544 545 do s1 = s2, nseeds 546 i = j+s1-s2 547 band_storage_matrix(1+i-j,j) = alpha(s1,s2,k) 548 end do 549 end do 550 551 ! exit when k = kmax, as this beta block does not contribute. 552 if (k .eq. kmax) exit 553 554 ! beta blocks 555 do s2 = 1, nseeds 556 j = (k-1)*nseeds+s2 557 558 do s1 = 1, s2 559 i = j+s1-s2+nseeds 560 band_storage_matrix(1+i-j,j) = beta(s1,s2,k) 561 end do 562 563 end do 564 end do 565 566 saved_band_storage_matrix(:,:) = band_storage_matrix(:,:) 567 568 !----------------------------------------- 569 ! Diagonalize the banded matrix 570 !----------------------------------------- 571 572 ABI_MALLOC(eigenvectors, (nseeds*kmax,nseeds*kmax)) 573 574 ABI_MALLOC(work,(nseeds*kmax)) 575 ABI_MALLOC(rwork,(3*nseeds*kmax-2)) 576 577 call ZHBEV( 'V', & ! compute eigenvalues and eigenvectors 578 'L', & ! lower triangular part of matrix is stored in banded_matrix 579 nseeds*kmax, & ! dimension of matrix 580 kd, & ! number of superdiagonals in banded matrix 581 band_storage_matrix, & ! matrix in banded storage 582 ldab, & ! leading dimension of banded_matrix 583 eigenvalues, & ! eigenvalues of matrix 584 eigenvectors, & ! eigenvectors of matrix 585 nseeds*kmax, & ! dimension of eigenvector matrix 586 work, rwork, info ) ! work arrays and info 587 588 589 if ( info /= 0) then 590 debug_unit = get_unit() 591 write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log' 592 593 open(debug_unit,file=trim(debug_filename),status='unknown') 594 595 write(debug_unit,'(A)') '*********************************************************************************************' 596 write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZHBEV (1), gwls_GWlanczos' 597 write(debug_unit,'(A)') '*********************************************************************************************' 598 599 close(debug_unit) 600 601 end if 602 603 604 605 !---------------------------------------------------------------------------------- 606 ! update the Lanczos basis to reflect diagonalization of T matrix 607 ! 608 ! Note that by definition 609 ! 610 ! Q^H . A . Q = T ==> A = Q . T . Q^H 611 ! 612 ! where Q (Lbasis) contains the Lanczos basis. 613 ! 614 ! Diagonalizing T, ie T = U . LAMBDA . U^H leads to 615 ! 616 ! A = [ Q.U] . LAMBDA . [Q.U]^H 617 ! 618 ! The updated basis is thus Q.U == Lbasis . eigenvectors 619 !---------------------------------------------------------------------------------- 620 621 ! NEVER use matmul!!! It sends temporary arrays to the stack, which can be much smaller 622 ! than needed; this leads to mysterious segfaults! 623 624 ! Lbasis = matmul(Lbasis,eigenvectors) 625 626 ! use temporary array, which is PROPERLY ALLOCATED, to perform matrix multiplication 627 ABI_MALLOC(Lbasis_tmp, (Hsize,nseeds*kmax)) 628 629 ! Compute C = A * B, where A = Lbasis, B = eigenvectors, and C = Lbasis_tmp 630 call ZGEMM( 'N', & ! leave array A as is 631 'N', & ! leave array B as is 632 Hsize, & ! number of rows of A 633 nseeds*kmax, & ! number of columns of B 634 nseeds*kmax, & ! number of columns of A /rows of B 635 cmplx_1, & ! constant alpha 636 Lbasis, & ! matrix A 637 Hsize, & ! LDA 638 eigenvectors, & ! matrix B 639 nseeds*kmax, & ! LDB 640 cmplx_0, & ! constant beta 641 Lbasis_tmp, & ! matrix C 642 Hsize) ! LDC 643 644 ! overwrite initial array 645 Lbasis(:,:) = Lbasis_tmp(:,:) 646 647 648 ABI_FREE(Lbasis_tmp) 649 650 if ( debug .and. mpi_enreg%me == 0) then 651 !---------------------------------------------------------------------------------- 652 ! For the purpose of debugging, print relevant results to file to check 653 ! data consistency. This may not be necessary once the code has been shown to 654 ! work properly. 655 !---------------------------------------------------------------------------------- 656 657 io_unit = get_unit() 658 i = 0 659 file_exists = .true. 660 do while (file_exists) 661 i = i+1 662 write(filename,'(A,I0.4,A)') "diagonalize_banded_matrix_",i,".log" 663 inquire(file=filename,exist=file_exists) 664 end do 665 666 667 open(io_unit,file=filename,status=files_status_new) 668 write(io_unit,10) "#=======================================================================================" 669 write(io_unit,10) "# " 670 write(io_unit,10) "# This file contains information pertaining to the diagonalization of a banded " 671 write(io_unit,10) "# matrix, expressed in terms of the Lanczos alpha and beta block arrays. " 672 write(io_unit,10) "# " 673 write(io_unit,10) "#=======================================================================================" 674 write(io_unit,10) "# " 675 write(io_unit,12) "# diagonalization info : ",info," " 676 write(io_unit,10) "# " 677 write(io_unit,10) "# l lambda_l " 678 write(io_unit,10) "#=======================================================================================" 679 680 do i = 1, nseeds*kmax 681 write(io_unit,13) i, eigenvalues(i) 682 end do 683 684 write(io_unit,10) " " 685 write(io_unit,10) "# " 686 write(io_unit,10) "# alpha and beta blocks " 687 write(io_unit,10) "# " 688 write(io_unit,10) "#=======================================================================================" 689 690 ! loop on all blocks 691 do k=1,kmax 692 write(io_unit,10) "# " 693 write(io_unit,12) "# block k = ",k," " 694 write(io_unit,10) "# " 695 write(io_unit,15) "# alpha: ||alpha^H-alpha|| = ", & 696 sqrt(sum(abs(alpha(:,:,k)-transpose(conjg(alpha(:,:,k))))**2)) 697 do s1 = 1, nseeds 698 write(io_unit,14) alpha(s1,:,k) 699 end do 700 701 write(io_unit,10) "# " 702 write(io_unit,10) "# beta " 703 do s1 = 1, nseeds 704 write(io_unit,14) beta(s1,:,k) 705 end do 706 707 708 end do 709 710 711 712 713 write(io_unit,10) "# " 714 write(io_unit,10) "# band storage matrix: " 715 write(io_unit,10) "#=======================================================================================" 716 717 do i = 1, ldab 718 write(io_unit,14) saved_band_storage_matrix(i,:) 719 end do 720 721 close(io_unit) 722 end if 723 724 725 ! clean up memory 726 ABI_FREE(eigenvectors) 727 ABI_FREE( work) 728 ABI_FREE(rwork) 729 ABI_FREE(band_storage_matrix) 730 ABI_FREE(saved_band_storage_matrix) 731 732 733 10 format(A) 734 12 format(A,I5,A) 735 13 format(I5,5X,F24.12) 736 14 format(4X,1000(F12.8,SP,F12.8,1X,'i',2X)) 737 15 format(A,ES24.8) 738 739 end subroutine diagonalize_lanczos_banded
m_gwls_GWlanczos/get_seeds [ Functions ]
[ Top ] [ m_gwls_GWlanczos ] [ Functions ]
NAME
get_seeds
FUNCTION
.
INPUTS
OUTPUT
SOURCE
75 subroutine get_seeds(first_seed, nseeds, seeds) 76 !---------------------------------------------------------------------------------------------------- 77 ! 78 ! This subroutine compute the seeds using the eigenstates of the Hamiltonian 79 ! 80 !---------------------------------------------------------------------------------------------------- 81 implicit none 82 83 integer, intent(in) :: first_seed, nseeds 84 complex(dpc), intent(out) :: seeds(npw_k,nseeds) 85 86 real(dp) , allocatable :: psik_out(:,:) 87 real(dp) , allocatable :: psikb_e(:,:) 88 real(dp) , allocatable :: psig_e(:,:) 89 real(dp) , allocatable :: psikb_s(:,:) 90 real(dp) , allocatable :: psig_s(:,:) 91 92 93 94 ! local variables 95 integer :: n 96 integer :: i, j, nsblk 97 98 ! ************************************************************************* 99 100 ! Generate the seeds for the Lanczos algorithm 101 ABI_MALLOC(psik_out,(2,npw_k)) 102 ABI_MALLOC(psikb_e,(2,npw_kb)) 103 ABI_MALLOC(psig_e,(2,npw_g)) 104 ABI_MALLOC(psikb_s,(2,npw_kb)) 105 ABI_MALLOC(psig_s,(2,npw_g)) 106 107 nsblk = ceiling(1.0*nseeds/blocksize) 108 109 do i=1,nsblk 110 do j=1,blocksize 111 psikb_e(:,(j-1)*npw_k+1:j*npw_k) = cg(:,(e-1)*npw_k+1:e*npw_k) 112 end do 113 114 psig_e = zero 115 call wf_block_distribute(psikb_e, psig_e,1) ! LA -> FFT 116 117 do j=1,blocksize 118 n = (i-1)*blocksize + j-1 + first_seed 119 if ((i-1)*blocksize + j <= nseeds) then 120 psikb_s(:,(j-1)*npw_k+1:j*npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k) 121 else 122 psikb_s(:,(j-1)*npw_k+1:j*npw_k) = zero 123 end if 124 end do 125 126 psig_s = zero 127 call wf_block_distribute(psikb_s, psig_s,1) ! LA -> FFT 128 129 ! Fourier transform valence wavefunction, to real space 130 call g_to_r(psir1,psig_s) 131 132 psir1(2,:,:,:) = -psir1(2,:,:,:) 133 134 call gr_to_g(psig_s,psir1,psig_e) 135 136 ! return to LA configuration, in order to apply Coulomb potential 137 call wf_block_distribute(psikb_s, psig_s,2) ! FFT -> LA 138 139 do j=1, blocksize 140 n = (i-1)*blocksize + j 141 if(n<=nseeds) then 142 psik_out = psikb_s(:,(j-1)*npw_k+1:j*npw_k) 143 call sqrt_vc_k(psik_out) 144 seeds(:,n) = cmplx_1*psik_out(1,:) + cmplx_i*psik_out(2,:) 145 end if 146 end do 147 end do 148 149 ABI_FREE(psik_out) 150 ABI_FREE(psikb_e) 151 ABI_FREE(psig_e) 152 ABI_FREE(psikb_s) 153 ABI_FREE(psig_s) 154 155 end subroutine get_seeds