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