TABLE OF CONTENTS


ABINIT/m_gwls_lineqsolver [ Modules ]

[ Top ] [ Modules ]

NAME

 m_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 m_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 m_gwls_wf
45 use m_gwls_hamiltonian
46 
47 ! abinit modules
48 use defs_basis
49 use m_abicore
50 use m_xmpi
51 use m_bandfft_kpt
52 use m_cgtools
53 
54 use m_time,      only : timab
55 use m_io_tools,  only : get_unit
56 
57 implicit none
58 save
59 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

 90 subroutine sqmr(b,x,lambda,project_on_what,omega,omega_imaginary,kill_Pc_x)
 91 !--------------------------------------------------------------------------------
 92 ! This subroutine solves the linear algebra problem
 93 !
 94 !                                 A x = b
 95 ! 
 96 ! Where:
 97 !                INPUT
 98 !                -----
 99 !        real(dp) b                          right-hand-side of the equation to be solved
100 !     real(dp) omega                         *OPTIONAL* frequency used in building A
101 !     logical  omega_imaginary               *OPTIONAL* is the frequency imaginary?
102 !     real(dp) lambda                        value to be subtracted from the Hamiltonian
103 !     integer  project_on_what               flag which determines the projection scheme.
104 !
105 !                OUTPUT
106 !                -----
107 !        real(dp) x                          solution
108 !
109 !    Note that blocksize corresponds to the number of band processors; it is a global
110 !    variable defined in gwls_hamiltonian. The name is inspired from lobpcgwf.F90.
111 !
112 ! with:
113 !        omega     omega_imaginary      Operator
114 !        ------------------------------------------------------
115 !     absent          -            A =   (H - lambda)  
116 !     present         -            A =   (H - lambda)^2 - omega^2
117 !     present    present, true     A =   (H - lambda)^2 + omega^2
118 !
119 !        project_on_what                        action
120 !        ------------------------------------------------------
121 !                0                no projection
122 !                1                projection on conduction states
123 !                2                projection out of subspace degenerate with lambda
124 !                3                projection on states beyond all the states explicitly stored
125 ! 
126 ! NOTE: It is the developper's responsibility to apply (H-ev) on the input 
127 !       if the frequency is not zero.
128 !--------------------------------------------------------------------------------
129 
130 !This section has been created automatically by the script Abilint (TD).
131 !Do not modify the following lines by hand.
132 #undef ABI_FUNC
133 #define ABI_FUNC 'sqmr'
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