module steppr_mod
  implicit none
  real   , dimension(:), allocatable, private :: s, ss
  complex, dimension(:), allocatable, private :: p, pp
  interface cgstep
    module procedure cgstep_real
    module procedure cgstep_cplx
  end interface
contains

  integer function cgstep_real(forget, x, g, rr, gg) result(stat)
    real, dimension(:) :: x, g, rr, gg
    logical            :: forget
    double precision   :: sds, gdg, gds, determ, gdr, sdr
    double precision   :: 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.d0 
      ! steepest descent
      if ( dot_product(gg, gg) .eq. 0 ) then
        call erexit('cgstep: grad vanishes identically')
      end if
      alfa = - sum(dprod(gg, rr)) / sum(dprod(gg, gg))
    else
      gdg = sum(dprod(gg, gg))
      ! search plane by solving 2-by-2
      sds = sum(dprod(ss, ss))
      ! G . (R - G*alfa - S*beta) = 0
      gds = sum(dprod(gg, ss))
      ! S . (R - G*alfa - S*beta) = 0
      if (gdg.eq.0. .or. sds.eq.0.) then
        stat = 1
        return
      end if
      determ = gdg * sds * max(1.d0-(gds/gdg)*(gds/sds), 1.d-12)
      gdr = -sum( dprod( gg, rr))
      sdr = -sum( dprod( ss, rr))
      alfa = ( sds * gdr - gds * sdr ) / determ
      beta = (-gds * gdr + gdg * sdr ) / determ
    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.
    stat = 0
  end function

  integer function cgstep_cplx(forget, x, g, rr, gg) result(stat)
    complex, dimension(:) :: x, g, rr, gg
    logical               :: forget
    double precision      :: sds, gdg, gds, determ, gdr, sdr
    double precision      :: alfa, beta
    if (.not.allocated (p)) then
      forget = .true.
      allocate (p (size(x )))
      allocate (pp(size(rr)))
    end if
    if (forget) then
      p  = 0.
      pp = 0.
      beta = 0.d0 
      ! steepest descent
      if ( dble(dot_product(gg, gg)) .eq. 0 ) then
        call erexit('cgstep: grad vanishes identically')
      end if
      alfa = - dble(dot_product(gg, rr))/dble(dot_product(gg,gg))
    else
      gdg = dble(dot_product(gg, gg))
      ! search plane by solving 2-by-2
      sds = dble(dot_product(pp, pp))
      ! Re{ G'. (R + G*alfa + S*beta)} = 0
      gds = dble(dot_product(gg, pp))
      ! Re{ S'. (R + G*alfa + S*beta)} = 0
      if ( gdg.eq.0. .or. sds.eq.0.) then
        stat = 1
        return
      end if
      determ = gdg*sds*max((1.d0-(gds/gdg)*(gds/sds)),1.d-12)
      gdr = -dble(dot_product(gg, rr))
      sdr = -dble(dot_product(pp, rr))
      alfa = ( sds * gdr - gds * sdr ) / determ
      beta = (-gds * gdr + gdg * sdr ) / determ
    end if
    p  = alfa * g  + beta * p  ! update solution step
    pp = alfa * gg + beta * pp ! update residual step
    x  = x  + p                ! update solution
    rr = rr + pp               ! update residual
    forget = .false.
    stat = 0
  end function

  subroutine cgstep_close ( )
    if (allocated(s)) deallocate(s, ss)
    if (allocated(p)) deallocate(p, pp)
  end subroutine

end module
