TABLE OF CONTENTS
ABINIT/gwls_lineqsolver [ Modules ]
NAME
gwls_lineqsolver
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 !--------------------------------------------------------------------- 30 ! Module to solve A.x = b efficiently, where A will involve 31 ! the Hamiltonian. 32 !--------------------------------------------------------------------- 33 34 35 module gwls_lineqsolver 36 !---------------------------------------------------------------------------------------------------- 37 ! This module contains routines to solve A x = b interatively to solve the Sternheimer equation 38 ! in various contexts... 39 !---------------------------------------------------------------------------------------------------- 40 41 42 ! local modules 43 use m_gwls_utility 44 use gwls_wf 45 use gwls_hamiltonian 46 47 ! abinit modules 48 use defs_basis 49 use m_profiling_abi 50 use m_xmpi 51 use m_bandfft_kpt 52 use m_cgtools 53 54 use m_io_tools, only : get_unit 55 56 implicit none 57 save 58 private
m_hamiltonian/qmr [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
qmr
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_polarisability
CHILDREN
hpsikc,precondition_cplx,unset_precondition,xmpi_sum
SOURCE
639 subroutine qmr(b,x,lambda) !,project_on_what) 640 !-------------------------------------------------------------------------------- 641 ! This subroutine solves the linear algebra problem 642 ! 643 ! A x = b 644 ! 645 ! where A := (H - lambda) can be non-hermitian. 646 ! Thus, complex values of lambda are allowed and 647 ! non-hermitian operators could be handled instead of H. 648 ! 649 ! Arguments : 650 ! INPUT 651 ! ----- 652 ! real(dp) b(2,npw_k) right-hand-side of the equation to be solved 653 ! real(dp) lambda(2) value to be subtracted from the Hamiltonian (complex). 654 ! integer project_on_what flag which determines the projection scheme. 655 ! 656 ! OUTPUT 657 ! ----- 658 ! real(dp) x(2,npw_k) solution 659 ! 660 ! project_on_what action 661 ! ------------------------------------------------------ 662 ! 0 no projection 663 ! 1 projection on conduction states 664 ! 2 projection out of subspace degenerate with lambda 665 ! 3 projection on states beyond all the states explicitly stored 666 !-------------------------------------------------------------------------------- 667 668 !This section has been created automatically by the script Abilint (TD). 669 !Do not modify the following lines by hand. 670 #undef ABI_FUNC 671 #define ABI_FUNC 'qmr' 672 !End of the abilint section 673 674 implicit none 675 676 !External variables 677 real(dp), intent(in) :: b(2,npw_k) 678 real(dp), intent(in) :: lambda(2) 679 real(dp), intent(out) :: x(2,npw_k) 680 !integer, intent(in) :: project_on_what !Unused yet, no projections done. 681 682 !Local variables 683 complex(dpc), allocatable :: xc(:), r(:), v(:), w(:), z(:), p(:), q(:), y(:), t(:), d(:), s(:) 684 complex(dpc), allocatable :: beta(:), eta(:), delta(:), epsilonn(:) 685 complex(dpc) :: lambdac 686 real(dp), allocatable :: rho(:), zeta(:), gama(:), theta(:), resid(:) 687 integer :: i 688 integer :: ierr 689 690 integer :: mpi_communicator 691 !logical :: precondition_on 692 693 ! ************************************************************************* 694 695 if(sum(b**2) < tol12) then 696 x=zero 697 else 698 699 !Allocation 700 ABI_ALLOCATE(xc,(npw_k)) 701 ABI_ALLOCATE(r ,(npw_k)) 702 ABI_ALLOCATE(v ,(npw_k)) 703 ABI_ALLOCATE(w ,(npw_k)) 704 ABI_ALLOCATE(z ,(npw_k)) 705 ABI_ALLOCATE(p ,(npw_k)) 706 ABI_ALLOCATE(q ,(npw_k)) 707 ABI_ALLOCATE(y ,(npw_k)) 708 ABI_ALLOCATE(t ,(npw_k)) 709 ABI_ALLOCATE(d ,(npw_k)) 710 ABI_ALLOCATE(s ,(npw_k)) 711 712 ABI_ALLOCATE(beta ,(nline)) 713 ABI_ALLOCATE(rho ,(nline+1)) 714 ABI_ALLOCATE(zeta ,(nline+1)) 715 ABI_ALLOCATE(gama ,(nline+1)) 716 ABI_ALLOCATE(eta ,(nline+1)) 717 ABI_ALLOCATE(theta ,(nline+1)) 718 ABI_ALLOCATE(delta ,(nline)) 719 ABI_ALLOCATE(epsilonn,(nline)) 720 ABI_ALLOCATE(resid ,(nline+1)) 721 722 !Initialization 723 x = zero 724 xc = zero 725 r = zero 726 v = zero 727 w = zero 728 z = zero 729 p = zero 730 q = zero 731 y = zero 732 t = zero 733 d = zero 734 s = zero 735 736 beta = zero 737 rho = zero 738 zeta = zero 739 gama = zero 740 eta = zero 741 theta = zero 742 delta = zero 743 epsilonn = zero 744 resid = zero 745 746 call unset_precondition() 747 748 749 !mpi_communicator = mpi_enreg%comm_fft 750 mpi_communicator = mpi_enreg%comm_bandfft 751 752 lambdac = dcmplx(lambda(1),lambda(2)) 753 754 i = 1 755 r = dcmplx(b(1,:),b(2,:)) 756 v = r 757 758 rho(i) = norm_kc(v) 759 call xmpi_sum(rho(i),mpi_communicator ,ierr) ! sum on all processors working on FFT! 760 761 w = r 762 call precondition_cplx(z,w) 763 zeta(i) = norm_kc(z) 764 call xmpi_sum(zeta(i),mpi_communicator,ierr) ! sum on all processors working on FFT! 765 766 gama(i) = one 767 eta(i) = -one 768 !theta(i) = zero 769 !p = zero 770 !q = zero 771 772 do i=1,nline 773 v = v/rho(i) 774 w = w/zeta(i) 775 z = z/zeta(i) 776 delta(i) = scprod_kc(z,v) 777 call xmpi_sum(delta(i),mpi_communicator,ierr) ! sum on all processors working on FFT! 778 779 call precondition_cplx(y,v) 780 if(i/=1) then 781 p = y - (zeta(i)*delta(i)/epsilonn(i-1))*p 782 q = z - ( rho(i)*delta(i)/epsilonn(i-1))*q 783 else 784 p = y 785 q = z 786 end if 787 call Hpsikc(t,p,lambdac) 788 epsilonn(i) = scprod_kc(q,t) 789 call xmpi_sum(epsilonn(i),mpi_communicator,ierr) ! sum on all processors working on FFT! 790 791 beta(i) = epsilonn(i)/delta(i) 792 v = t - beta(i)*v 793 rho(i+1) = norm_kc(v) 794 call xmpi_sum(rho(i+1),mpi_communicator,ierr) ! sum on all processors working on FFT! 795 796 call Hpsikc(z,q,conjg(lambdac)) ; w = z - beta(i)*w 797 call precondition_cplx(z,w) 798 zeta(i+1) = norm_kc(z) 799 call xmpi_sum(zeta(i+1), mpi_communicator,ierr) ! sum on all processors working on FFT! 800 801 theta(i+1) = rho(i+1)/(gama(i)*abs(beta(i))) 802 gama(i+1) = 1./sqrt(1+theta(i+1)**2) 803 eta(i+1) = -eta(i)*rho(i)*gama(i+1)**2/(beta(i)*gama(i)**2) 804 d = eta(i+1)*p + ((theta(i)*gama(i+1))**2)*d 805 s = eta(i+1)*t + ((theta(i)*gama(i+1))**2)*s 806 xc = xc + d 807 r = r - s 808 resid(i) = norm_kc(r)**2 809 call xmpi_sum(resid(i), mpi_communicator,ierr) ! sum on all processors working on FFT! 810 811 !write(std_out,*) "QMR residual**2 = ",resid(i),"; i = ",i-1 812 if(resid(i) < tolwfr) exit 813 end do 814 815 if(i>=nline) then 816 write(std_out,*) " **** Iterations were not enough to converge! ****" 817 end if 818 819 call Hpsikc(r,xc,lambdac) 820 v = dcmplx(b(1,:),b(2,:)) 821 resid(nline+1) = norm_kc(r - v)**2 822 call xmpi_sum(resid(nline+1), mpi_communicator,ierr) ! sum on all processors working on FFT! 823 824 write(std_out,*) "QMR residual**2 (at end) = ",resid(nline+1),"; # iterations = ",i-1 825 826 x(1,:) = dble(xc) 827 x(2,:) = dimag(xc) 828 829 !Deallocate 830 ABI_DEALLOCATE(xc) 831 ABI_DEALLOCATE(r) 832 ABI_DEALLOCATE(v) 833 ABI_DEALLOCATE(w) 834 ABI_DEALLOCATE(z) 835 ABI_DEALLOCATE(p) 836 ABI_DEALLOCATE(q) 837 ABI_DEALLOCATE(y) 838 ABI_DEALLOCATE(t) 839 ABI_DEALLOCATE(d) 840 ABI_DEALLOCATE(s) 841 842 ABI_DEALLOCATE(beta ) 843 ABI_DEALLOCATE(rho ) 844 ABI_DEALLOCATE(zeta ) 845 ABI_DEALLOCATE(gama ) 846 ABI_DEALLOCATE(eta ) 847 ABI_DEALLOCATE(theta ) 848 ABI_DEALLOCATE(delta ) 849 ABI_DEALLOCATE(epsilonn) 850 ABI_DEALLOCATE(resid ) 851 852 end if 853 854 end subroutine qmr
m_hamiltonian/sqmr [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
sqmr
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_DielectricArray,gwls_polarisability
CHILDREN
hpsikc,precondition_cplx,unset_precondition,xmpi_sum
SOURCE
89 subroutine sqmr(b,x,lambda,project_on_what,omega,omega_imaginary,kill_Pc_x) 90 !-------------------------------------------------------------------------------- 91 ! This subroutine solves the linear algebra problem 92 ! 93 ! A x = b 94 ! 95 ! Where: 96 ! INPUT 97 ! ----- 98 ! real(dp) b right-hand-side of the equation to be solved 99 ! real(dp) omega *OPTIONAL* frequency used in building A 100 ! logical omega_imaginary *OPTIONAL* is the frequency imaginary? 101 ! real(dp) lambda value to be subtracted from the Hamiltonian 102 ! integer project_on_what flag which determines the projection scheme. 103 ! 104 ! OUTPUT 105 ! ----- 106 ! real(dp) x solution 107 ! 108 ! Note that blocksize corresponds to the number of band processors; it is a global 109 ! variable defined in gwls_hamiltonian. The name is inspired from lobpcgwf.F90. 110 ! 111 ! with: 112 ! omega omega_imaginary Operator 113 ! ------------------------------------------------------ 114 ! absent - A = (H - lambda) 115 ! present - A = (H - lambda)^2 - omega^2 116 ! present present, true A = (H - lambda)^2 + omega^2 117 ! 118 ! project_on_what action 119 ! ------------------------------------------------------ 120 ! 0 no projection 121 ! 1 projection on conduction states 122 ! 2 projection out of subspace degenerate with lambda 123 ! 3 projection on states beyond all the states explicitly stored 124 ! 125 ! NOTE: It is the developper's responsibility to apply (H-ev) on the input 126 ! if the frequency is not zero. 127 !-------------------------------------------------------------------------------- 128 129 !This section has been created automatically by the script Abilint (TD). 130 !Do not modify the following lines by hand. 131 #undef ABI_FUNC 132 #define ABI_FUNC 'sqmr' 133 use interfaces_18_timing 134 !End of the abilint section 135 136 implicit none 137 138 !External variables 139 real(dp), intent(in) :: b(2,npw_g) 140 real(dp), intent(in) :: lambda 141 real(dp), intent(out) :: x(2,npw_g) 142 integer, intent(in) :: project_on_what 143 real(dp), intent(in), optional :: omega 144 logical, optional :: omega_imaginary, kill_Pc_x 145 146 !Local variables 147 real(dp) :: norm, tmp(2), residual 148 real(dp), allocatable :: g(:), theta(:), rho(:), sigma(:), c(:) 149 real(dp), allocatable :: t(:,:), delta(:,:), r(:,:), d(:,:), w(:,:), wmb(:,:) 150 integer :: k, l 151 real(dp):: signe 152 real(dp):: norm_Axb 153 154 155 real(dp):: norm_b, tol14 156 157 integer :: min_index 158 logical :: singular 159 logical :: precondition_on 160 161 162 integer :: pow 163 164 logical :: imaginary 165 166 integer,save :: counter = 0 167 integer :: io_unit 168 character(128) :: filename 169 logical :: file_exists 170 logical :: head_node 171 172 integer :: ierr 173 174 integer :: mpi_communicator, mpi_rank, mpi_group 175 176 177 178 ! timing 179 real(dp) :: tsec(2) 180 integer :: GWLS_TIMAB, OPTION_TIMAB 181 182 ! ************************************************************************* 183 184 ! The processors communicate over FFT! 185 mpi_communicator = mpi_enreg%comm_fft 186 187 ! what is the rank of this processor, within its group? 188 mpi_rank = mpi_enreg%me_fft 189 190 ! Which group does this processor belong to, given the communicator? 191 mpi_group = mpi_enreg%me_band 192 193 194 ! Test if the input has finite norm 195 tol14 = 1.0D-14 196 tmp = cg_zdotc(npw_g,b,b) 197 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 198 norm_b = tmp(1) 199 200 if (norm_b < tol14) then 201 ! Because of band parallelism, it is possible that sqmr gets a zero norm argument. 202 ! A | x> = 0 implies |x > = 0. 203 x(:,:) = zero 204 return 205 end if 206 207 GWLS_TIMAB = 1523 208 OPTION_TIMAB = 1 209 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 210 211 ! only the head node should write to the log file 212 head_node = ( mpi_rank == 0 ) 213 214 !Memory allocation for local variables 215 ABI_ALLOCATE(g, (nline)) 216 ABI_ALLOCATE(theta,(nline)) 217 ABI_ALLOCATE(rho, (nline)) 218 ABI_ALLOCATE(sigma,(nline)) 219 ABI_ALLOCATE(c, (nline)) 220 221 ABI_ALLOCATE(t, (2,npw_g)) 222 ABI_ALLOCATE(delta,(2,npw_g)) 223 ABI_ALLOCATE(r, (2,npw_g)) 224 ABI_ALLOCATE(d, (2,npw_g)) 225 ABI_ALLOCATE(w, (2,npw_g)) 226 ABI_ALLOCATE(wmb, (2,npw_g)) 227 228 229 230 !Some vectors won't be filled (first iteration missing) so it's useful to initialise them. 231 g = zero 232 theta = zero 233 rho = zero 234 sigma = zero 235 c = zero 236 x = zero 237 delta = zero 238 r = zero 239 d = zero 240 w = zero 241 242 243 ! Determine if the frequency is imaginary 244 if (present(omega_imaginary) .and. present(omega)) then 245 imaginary = omega_imaginary 246 else 247 imaginary = .false. 248 end if 249 250 251 ! Define the sign in front of (H-eig(v))**2. 252 ! If omega_imaginary is not given, we assume that omega is real (and sign=-1). 253 if (present(omega)) then 254 if ( imaginary ) then 255 signe = one 256 else 257 signe =-one 258 end if 259 end if 260 261 262 !Check for singularity problems 263 if (present(omega)) then 264 norm = minval(abs((eig(1:nbandv)-lambda)**2 + signe*(omega)**2)) 265 min_index = minloc(abs((eig(1:nbandv)-lambda)**2 + signe*(omega)**2),1) 266 else 267 norm = minval(abs(eig(1:nbandv)-lambda)) 268 min_index = minloc(abs(eig(1:nbandv)-lambda),1) 269 end if 270 singular = norm < 1.0d-12 271 272 !-------------------------------------------------------------------------------- 273 ! If the linear operator has a kernel, then the intermediate vectors obtained in 274 ! SQMR must be projected out of this subspace several time at each iterations, 275 ! otherwise SQMR is unstable. 276 ! 277 ! This is true even if the seed vector has been initially projected out of this 278 ! subspace, since the preconditionning will re-introduce a non-zero component in 279 ! the subspace of the kernel of the linear operator. 280 ! ===> Use project_on_what==2 in such cases. 281 ! 282 ! Here, test if the operator is singular and if we are NOT projecting out of 283 ! the kernel. 284 ! ===> If true, stop the code. 285 ! 286 !-------------------------------------------------------------------------------- 287 288 ! Quit if the operator has an uncontrolled kernel, a sign that the routine is being 289 ! misused by a developper... 290 if (singular .and. ( (project_on_what==1 .and. (min_index > nbandv)) .or. project_on_what==0 )) then 291 write(std_out,*) "ERROR - SQMR: Quasi-singuar problem treated, min. eigenvalue of A is ", norm," < 1d-12." 292 write(std_out,*) " Yet, there is no projection out of the kernel of A. " 293 294 if (project_on_what==1 .and. (min_index > nbandv)) then 295 write(std_out,*) " " 296 write(std_out,*) " There is a projection on the conduction states, but A is singular in this " 297 write(std_out,*) " subspace (the kernel contains state i=",min_index," > ",nbandv,"=# of valence states)." 298 end if 299 300 write(std_out,*) " " 301 write(std_out,*) " In this situation, SQMR will be unstable. Use project_on_what==2 as an " 302 write(std_out,*) " input argument of SQMR." 303 write(std_out,*) " " 304 write(std_out,*) " Decision taken to exit..." 305 stop 306 end if 307 308 !-------------------------------------------------------------------------------- 309 ! Open a log file for the output of SQMR; only write if head of group! 310 !-------------------------------------------------------------------------------- 311 if (head_node) then 312 313 io_unit = get_unit() 314 315 write(filename,'(A,I0.4,A)') "SQMR_GROUP=",mpi_group,".log" 316 317 inquire(file=filename,exist=file_exists) 318 319 if (file_exists) then 320 open( io_unit,file=filename,position='append',status=files_status_old) 321 else 322 open( io_unit,file=filename,status=files_status_new) 323 write(io_unit,10) "#=======================================================================================" 324 write(io_unit,10) "# " 325 write(io_unit,10) "# This file contains information regarding the application of the SQMR scheme, " 326 write(io_unit,10) "# for this MPI group. " 327 write(io_unit,10) "#=======================================================================================" 328 flush(io_unit) 329 end if 330 331 counter = counter + 1 332 write(io_unit,10) "# " 333 write(io_unit,11) "# Call # ", counter 334 write(io_unit,12) "# lambda = ",lambda," Ha " 335 if (present(omega)) then 336 write(io_unit,12) "# omega = ",omega," Ha " 337 if (imaginary) then 338 write(io_unit,10) "# omega is imaginary " 339 else 340 write(io_unit,10) "# omega is real " 341 end if 342 else 343 write(io_unit,10) "# omega is absent " 344 end if 345 346 write(io_unit,13) "# project_on_what = ",project_on_what," " 347 write(io_unit,13) "# " 348 if (present(omega) ) then 349 if (imaginary) then 350 write(io_unit,10) "# SOLVE ((H-lambda)^2 + omega^2) x = b" 351 else 352 write(io_unit,10) "# SOLVE ((H-lambda)^2 - omega^2) x = b" 353 end if 354 else 355 write(io_unit,10) "# SOLVE (H-lambda) x = b" 356 end if 357 358 flush(io_unit) 359 end if ! head_node 360 !-------------------------------------------------------------------------------- 361 ! Precondition to accelerate convergence 362 !-------------------------------------------------------------------------------- 363 precondition_on = .true. 364 if(imaginary) then 365 if(omega > 10.0_dp) then 366 precondition_on = .false. 367 end if 368 end if 369 370 ! DEBUG 371 pow = project_on_what 372 !precondition_on = .false. 373 374 !Prepare to precondition 375 if (precondition_on) then 376 if ( imaginary ) then 377 call set_precondition(lambda,omega) 378 else 379 call set_precondition() 380 end if 381 else 382 call unset_precondition() 383 end if 384 385 !-------------------------------------------------------------------------------- 386 !Initialisation 387 !-------------------------------------------------------------------------------- 388 389 k = 1 390 l = 1 391 r = b 392 393 if (head_node) then 394 write(io_unit,10) "# " 395 write(io_unit,10) "# iteration approximate residual" 396 write(io_unit,10) "#----------------------------------------" 397 flush(io_unit) 398 end if 399 400 do ! outer loop 401 call precondition(d,r) 402 403 ! g(k) = norm_k(r) 404 tmp = cg_zdotc(npw_g,r,r) 405 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 406 g(k) = dsqrt(tmp(1)) 407 408 409 !tmp = scprod_k(r,d) 410 tmp = cg_zdotc(npw_g,r,d) 411 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 412 rho(k) = tmp(1) 413 414 if (head_node) then 415 write(io_unit,16) k, g(k)**2 416 flush(io_unit) 417 end if 418 419 do ! inner loop 420 k=k+1 421 l=l+1 422 423 ! Apply the A operator 424 if (present(omega)) then 425 call Hpsik(w,d,lambda) 426 call Hpsik(w,cte=lambda) 427 w = w + d*signe*omega**2 428 else 429 call Hpsik(w,d,lambda) 430 end if 431 432 ! Apply projections, if requested 433 !if(dtset%gwcalctyp /= 2) then !This is a test to obtain the time taken by the orthos. 434 if(pow == 1) call pc_k_valence_kernel(w) 435 !if(pow == 2) call pc_k(w,eig_e=lambda) 436 !if(pow == 3) call pc_k(w,n=nband,above=.true.) 437 !end if 438 439 ! Apply SQMR scheme 440 !tmp = scprod_k(d,w) 441 tmp = cg_zdotc(npw_g,d,w) 442 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 443 444 sigma(k-1) = tmp(1) 445 r = r-(rho(k-1)/sigma(k-1))*w 446 447 ! The following two lines must have a bug! We cannot distribute the norm this way! 448 ! theta(k) = norm_k(r)/g(k-1) 449 ! call xmpi_sum(theta(k), mpi_communicator,ierr) ! sum on all processors working on FFT! 450 tmp = cg_zdotc(npw_g,r,r) 451 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 452 theta(k) = dsqrt(tmp(1))/g(k-1) 453 454 c(k) = one/dsqrt(one+theta(k)**2) 455 g(k) = g(k-1)*theta(k)*c(k) 456 delta = delta*(c(k)*theta(k-1))**2+d*rho(k-1)/sigma(k-1)*c(k)**2 457 458 x = x+delta 459 460 461 if (head_node) then 462 write(io_unit,16) k, g(k)**2 463 flush(io_unit) 464 end if 465 466 ! Test exit condition 467 if(g(k)**2<tolwfr .or. k>= nline) exit 468 !if(k>=nline) exit 469 470 ! Safety test every 100 iterations, check that estimated residual is of the right order of magnitude. 471 ! If not, restart SQMR. 472 if(mod(l,100)==0) then 473 if(present(omega)) then 474 call Hpsik(w,x,lambda) 475 call Hpsik(w,cte=lambda) 476 w = w + x*signe*omega**2 477 else 478 call Hpsik(w,x,lambda) 479 end if 480 481 if(pow == 1) call pc_k_valence_kernel(w) 482 !if(pow == 2) call pc_k(w,eig_e=lambda) 483 !if(pow == 3) call pc_k(w,n=nband,above=.true.) 484 485 !if(norm_k(w-b)**2 > 10*g(k)**2) exit 486 wmb = w-b 487 tmp = cg_zdotc(npw_g,wmb,wmb) 488 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 489 if(tmp(1) > 10*g(k)**2) exit 490 491 end if 492 493 ! Get ready for next cycle 494 call precondition(w,r) 495 !if(dtset%gwcalctyp /= 2) then 496 if(pow == 1) call pc_k_valence_kernel(w) 497 !if(pow == 2) call pc_k(w,eig_e=lambda) 498 !if(pow == 3) call pc_k(w,n=nband,above=.true.) 499 !end if 500 501 !tmp = scprod_k(r,w) 502 tmp = cg_zdotc(npw_g,r,w) 503 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 504 rho(k) = tmp(1) 505 d = w+d*rho(k)/rho(k-1) 506 507 end do ! end inner loop 508 509 ! Exit condition 510 if(g(k)**2<tolwfr .or. k>=nline) exit 511 !if(k>=nline) exit 512 513 if (head_node) write(io_unit,10) " ---- RESTART of SQMR -----" 514 515 !norm_Axb = norm_k(w-b)**2 516 wmb = w-b 517 tmp = cg_zdotc(npw_g,wmb,wmb) 518 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 519 norm_Axb = tmp(1) 520 521 if (head_node) then 522 write(io_unit,*) "|Ax-b|^2 :",norm_Axb 523 write(io_unit,*) "g(k)^2 :",g(k)**2 524 flush(io_unit) 525 end if 526 527 k = k+1 528 l = 1 529 530 ! Apply the operator 531 if(present(omega)) then 532 call Hpsik(r,x,lambda) 533 call Hpsik(r,cte=lambda) 534 r = r + x*signe*omega**2 535 else 536 call Hpsik(r,x,lambda) 537 end if 538 539 if(pow == 1) call pc_k_valence_kernel(r) 540 !if(pow == 2) call pc_k(r,eig_e=lambda) 541 !if(pow == 3) call pc_k(r,n=nband,above=.true.) 542 543 r=b-r 544 end do 545 546 547 ktot = ktot+k 548 if(k >= nline .and. head_node ) then 549 write(io_unit,10) " **** Iterations were not enough to converge! ****" 550 end if 551 552 if( present(kill_Pc_x) ) then 553 if (.not. kill_Pc_x .and. pow == 1) call pc_k_valence_kernel(x) 554 end if 555 556 if( .not. present(kill_Pc_x) .and. pow == 1 ) call pc_k_valence_kernel(x) 557 558 559 560 561 !if(pow == 2) call pc_k(x,eig_e=lambda) 562 !if(pow == 3) call pc_k(x,n=nband,above=.true.) 563 564 if(present(omega)) then 565 call Hpsik(w,x,lambda) 566 call Hpsik(w,cte=lambda) 567 w = w + x*signe*omega**2 568 else 569 call Hpsik(w,x,lambda) 570 end if 571 if(pow == 1) call pc_k_valence_kernel(w) 572 !if(pow == 2) call pc_k(w,eig_e=lambda) 573 !if(pow == 3) call pc_k(w,n=nband,above=.true.) 574 575 !residual = norm_k(w-b)**2 576 wmb = w-b 577 tmp = cg_zdotc(npw_g,wmb,wmb) 578 call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT! 579 residual = tmp(1) 580 581 if (head_node) then 582 write(io_unit,15) "iterations :", k 583 write(io_unit,14) "tolwfr :", tolwfr 584 write(io_unit,14) "residuals (estimated) :", g(k)**2 585 write(io_unit,14) "residuals : |Ax-b|^2 :", residual 586 close(io_unit) 587 end if 588 589 590 591 ABI_DEALLOCATE(g) 592 ABI_DEALLOCATE(theta) 593 ABI_DEALLOCATE(rho) 594 ABI_DEALLOCATE(sigma) 595 ABI_DEALLOCATE(c) 596 ABI_DEALLOCATE(t) 597 ABI_DEALLOCATE(delta) 598 ABI_DEALLOCATE(r) 599 ABI_DEALLOCATE(d) 600 ABI_DEALLOCATE(w) 601 ABI_DEALLOCATE(wmb) 602 603 604 OPTION_TIMAB = 2 605 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 606 607 608 609 10 format(A) 610 11 format(A,I8) 611 12 format(A,E24.16,A) 612 13 format(A,I2,A) 613 14 format(20X,A,E24.16) 614 15 format(20X,A,I8) 615 16 format(5X,I5,15X,E12.3) 616 617 end subroutine sqmr