module ccgstep_mod
  implicit none
  complex, dimension (:), allocatable, private :: s, ss

contains

  integer function ccgstep( forget, x, g, rr, gg)
    complex, dimension (:) :: x, g, rr, gg
    logical                :: forget
    complex                :: sds, gdg, gds, determ, gdr, sdr, alfa, beta

    if ( .not. allocated (s)) then
      forget = .true.
      allocate ( s (size ( x)))
      allocate (ss (size (rr)))
    end if

    if ( forget) then
      s = 0.
      ss = 0.
      beta = 0. ! steepest descent
      if ( dot_product(gg, gg) .eq. 0 ) then
        call erexit('ccgstep: grad vanishes identically')
      end if
      alfa = - dot_product( gg, rr) / dot_product( gg, gg)
    else
      gdg = dot_product( gg, gg)
      ! search plane by solving 2-by-2
      sds = dot_product( ss, ss)
      ! G . (R - G*alfa - S*beta) = 0
      gds = dot_product( gg, ss)
      ! S . (R - G*alfa - S*beta) = 0
      if ( gdg.eq.0. .or. sds.eq.0.) then
        ccgstep = 1
        return
      end if
      determ = gdg * sds * (1. - (gds/gdg)*(gds/sds))
      gdr = - dot_product( gg, rr)
      sdr = - dot_product( ss, rr)
      alfa = ( sds * gdr - gds * sdr )*conjg(determ) / (cabs(determ)+0.000001)
      beta = (-gds * gdr + gdg * sdr )*conjg(determ) / (cabs(determ)+0.000001)
    end if
    s = alfa * g + beta * s      ! update solution step
    ss = alfa * gg + beta * ss   ! update residual step
    x = x + s                    ! update solution
    rr = rr + ss                 ! update residual
    forget = .false.
    ccgstep = 0
  end function

  subroutine ccgstep_close ( )
    if ( allocated( s)) then
      deallocate( s, ss)
    end if
  end subroutine
end module
