module cgstep_mod1
  implicit none
  real, dimension (:), allocatable, private   :: s, ss
  contains
  integer function cgstep( forget, x, g, rr, gg)
    real, dimension (:)  :: x, g, rr, gg
    logical,intent(inout)              :: forget 
    double precision     :: sds, gdg, gds, determ, gdr, sdr, alfa, beta

    if ( .not. allocated (s)) then
      allocate ( s (size ( x))) 
      allocate (ss (size (rr)))
    end if
    write(0,*) "passei 00",size(x),size(rr),allocated (s)
    if (allocated (s)) then
       write(0,*) "passei 01"
       forget=.true.
       write(0,*) "passei 02"
    end if
    if ( forget) then
      write(0,*) "passei 10"
      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))
      write(0,*) "passei 11"
    else
      write(0,*) "passei 20"
      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
      write(0,*) "passei 21"
    end if 
    write(0,*) "passei 30"
    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
    write(0,*) "passei 31"
    if ( allocated( s)) then
      deallocate( s, ss)
    end if
  end function 
end module 
