# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"

module cgstep_mod
  implicit none
  real, dimension (:), allocatable, private :: s, ss
  contains
  integer function cgstep( forget, x, g, rr, gg)
    real, dimension (:) :: x, g, rr, gg
    logical :: forget
    double precision :: 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.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
        cgstep = 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.
    cgstep = 0
  end function
  subroutine cgstep_close ( )
    if ( allocated( s)) then
      deallocate( s, ss)
    end if
  end subroutine
end module