module cgstep_mod001 
  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,soma,dgr,dgg
    integer              :: i

write(0,*) "cgstep000: passei 01"
    if ( .not. allocated (s)) then
write(0,*) "cgstep000: passei 0a"
      forget = .true.
write(0,*) "cgstep000: passei 0b"
      allocate ( s (size ( x))) 
write(0,*) "cgstep000: passei 0c"
      allocate (ss (size (rr)))
    end if
write(0,*) "cgstep001: passei 01"
    if ( forget) then
      s = 0.
      ss = 0.
      beta = 0.d0    ! steepest descent
write(0,*) "cgstep002: passei 01"

!!!!!!!!!!!!!!!!!!!!!!!
!      soma=dble(0.)
!      DO i = 1, size(gg)
!         if (not(isnan(gg(i)))) soma = soma + dble(gg(i)) * dble(gg(i))
!      END DO
!      dgg=soma

!      soma=dble(0.)
!      DO i = 1, size(gg)
!         if (not(isnan(gg(i))).and.not(isnan(rr(i)))) then
!            soma = soma + dble(rr(i)) * dble(gg(i))
!         end if
!      END DO
!      dgr=soma
!!!!!!!!!!!!!!!!!!!!!!!!
      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,*) "cgstep: passei 01", forget, alfa
    else
write(0,*) "2) cgstep: passei 000", forget
      gdg = sum( dprod( gg, gg))       
      ! search plane by solving 2-by-2
write(0,*) "2) cgstep: passei 001", forget,gdg
      sds = sum( dprod( ss, ss))       
write(0,*) "2) cgstep: passei 002", forget,sds
      !  G . (R - G*alfa - S*beta) = 0
      gds = sum( dprod( gg, ss))       
write(0,*) "2) cgstep: passei 003", forget,gdg
      !  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,*) "cgstep: passei 02", forget, alfa,beta
    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
write(0,*) "cgstep: passei 03", forget
  end function 
  subroutine cgstep_close ( )
    if ( allocated( s)) then
      deallocate( s, ss)
    end if
  end subroutine 
end module 
