!
  !-----------------------------------------------------------------------
  subroutine bcgsolve_all (Ax, e, b, x, h_diag, &
     ndmx, ndim, ethr, ik, kter, conv_root, anorm, nbnd, g2kin, vr, evq)
  !-----------------------------------------------------------------------
  !
  !  rewritten by FG following the book
  !  "Templates for the solution of linear systems: building blocks
  !  for iterative methods", by R. Barret et al, SIAM (http://www.siam.org/books)
  !  pag. 13. The only changes are:
  !  1. the product r^T*z on line 4 has been rewritten as z^T*r
  !  t obe consistent with the BiConjugate gradient method on pag. 20
  !  (look for [*] below).
  !  2. Instead of introducing p in the if/then/else, we overwrite z.
  !  Therefore our z at some point corresponds to their p.
  !
  !-----------------------------------------------------------------------
  !
  !     iterative solution of the linear system:
  !
  !                 ( h - e + Q + w + i * eta ) * x = b         
  !
  !                 where h is a complex hermitean matrix, e is a real sca
  !                 x and b are complex vectors
  !
  !     on input:
  !                 Ax    EXTERNAL  name of a subroutine:
  !                          Ax(ndim,psi,psip)
  !                          Calculates  H*psi products.
  !                          Vectors psi and psip should be dimensined
  !                          (ndmx,nvec). nvec=1 is used!
  !
  !                 cg_psi   EXTERNAL  name of a subroutine:
  !                          g_psi(ndmx,ndim,notcnv,psi,e)
  !                          which calculates (h-e)^-1 * psi, with
  !                          some approximation, e.g. (diag(h)-e)
  !
  !                 e        real     unperturbed eigenvalue.
  !
  !                 x     contains an estimate of the solution
  !                          vector.
  !
  !                 b    contains the right hand side vector
  !                          of the system.
  !
  !                 ndmx     integer row dimension of x, ecc.
  !
  !                 ndim     integer actual row dimension of x
  !
  !                 ethr     real     convergence threshold. solution
  !                          improvement is stopped when the error in
  !                          eq (1), defined as l.h.s. - r.h.s., becomes
  !                          less than ethr in norm.
  !
  !     on output:  x     contains the refined estimate of the
  !                          solution vector.
  !
  !                 b    is corrupted on exit
  !
  !   revised (extensively)       6 Apr 1997 by A. Dal Corso & F. Mauri
  !   revised (to reduce memory) 29 May 2004 by S. de Gironcoli
  !
  !-----------------------------------------------------------------------
  !
  use parameters, only : DP, nbnd_occ, eta
  use gspace, only : nr, ngm
  use constants
  implicit none
  integer :: i
  complex(DP) :: vr(nr)
  real(DP) :: g2kin (ngm)
  complex(kind=DP) :: evq (ngm, nbnd_occ)
  !
  !   first the I/O variables
  !
  integer :: ndmx, & ! input: the maximum dimension of the vectors
             ndim, & ! input: the actual dimension of the vectors
             kter, & ! output: counter on iterations
             nbnd, & ! input: the number of bands
             ik      ! input: the k point

  real(kind=DP) :: &
             e(nbnd), &        ! input: the actual eigenvalue
             anorm,   &        ! output: the norm of the error in the solution
             h_diag(ndmx,nbnd), & ! input: preconditioner
             ethr               ! input: the required precision
  complex(kind=DP) :: &
             x (ndmx, nbnd), & ! output: the solution of the linear syst
             y (ndmx, nbnd), & ! output: the solution of the linear syst
             b (ndmx, nbnd), & ! input: the known term
             bt (ndmx, nbnd)    ! input: the known term
  logical :: conv_root         ! output: if true the root is converged
  external Ax                  ! input: the routine computing Ax
  !
  ! local variables
  !
  integer, parameter :: maxter = 200
  integer :: iter, ibnd, lbnd
  complex(kind=DP), allocatable :: r  (:,:), q  (:,:), p  (:,:), pold  (:,:)
  complex(kind=DP), allocatable :: rt (:,:), qt (:,:), pt (:,:), ptold (:,:)
  complex(kind=DP) :: a, c, beta, alpha, ZDOTC
  complex(kind=DP) :: rho (nbnd), rhoold (nbnd)
  real(kind=DP) :: eu (nbnd)
  integer :: conv (nbnd)
  real(kind=DP) :: kter_eff
  !
  allocate ( r (ndmx,nbnd), q (ndmx,nbnd), p (ndmx,nbnd), pold (ndmx ,nbnd) )    
  allocate ( rt(ndmx,nbnd), qt(ndmx,nbnd), pt(ndmx,nbnd), ptold(ndmx ,nbnd) )    
  !
  kter_eff = zero
  conv = 0
  !
  do ibnd = 1, nbnd
    !
    do iter = 1, maxter
       !
       if (iter .eq. 1) then
          call Ax ( x(1,ibnd), r (1,ibnd) , e(ibnd), ik, 1, g2kin, vr, evq, + eta)  
          r  (:,ibnd) = b (:,ibnd) - r  (:,ibnd)  
          rt (:,ibnd) = conjg ( r  (:,ibnd) )
       endif
       !
       anorm = sqrt ( abs ( ZDOTC (ndim, r(1,ibnd), 1, r(1,ibnd), 1)  ) )
       if (anorm.lt.ethr) goto 100
!      write(6,*) anorm
       !
       if (iter.eq.1) then
         p   (:, ibnd) = r   (:, ibnd) 
         pt  (:, ibnd) = rt  (:, ibnd) 
       else
         a = ZDOTC (ndim, qt(1,ibnd), 1, r(1,ibnd), 1)
         c = ZDOTC (ndim, pt(1,ibnd), 1, q(1,ibnd), 1)
         beta = - a / c
         p   (:, ibnd) = r   (:, ibnd) + beta        * pold   (:, ibnd)
         pt  (:, ibnd) = rt  (:, ibnd) + conjg(beta) * ptold  (:, ibnd)
       endif
       !
       call Ax ( p (1,ibnd), q (1,ibnd), e(ibnd), ik, 1, g2kin, vr, evq, + eta )
       call Ax ( pt(1,ibnd), qt(1,ibnd), e(ibnd), ik, 1, g2kin, vr, evq, - eta )
       !
       a = ZDOTC (ndim, rt(1,ibnd), 1, r(1,ibnd), 1)
       c = ZDOTC (ndim, pt(1,ibnd), 1, q(1,ibnd), 1)
       alpha = a / c
       !
       x  (:,ibnd) = x  (:,ibnd) + alpha        * p  (:,ibnd)
       !
       r  (:,ibnd) = r  (:,ibnd) - alpha        * q  (:,ibnd)
       rt (:,ibnd) = rt (:,ibnd) - conjg(alpha) * qt (:,ibnd)
       !
       pold = p
       ptold = pt
       !
    enddo
100 continue
    !
  enddo
  !
  return
  end subroutine bcgsolve_all
  !