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