! Copyright 2007 Stanford University
!All rights reserved
!
!Author: Claudio Guerra (modified from cgstep by Sergey Fomel)
!
module cgstep1_mod   
  implicit none
  real, dimension (:), allocatable, private	   :: s, ss, dr, ra
  logical, parameter, private			   :: T = .true., F = .false.
  integer					   :: stat
  contains

  integer function cgstep1( forget, x, g, rr, gg, Fop)
  interface
       integer function Fop(adj,add,m,d) result(stat)
         real::m(:),d(:)
         logical,intent(in)::adj,add
       end function 
  end interface
  real, dimension (:)  :: x, g, rr, gg
  logical              :: forget 
  real  		 :: alfa, beta

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

    if ( forget) then

      s = 0.
      ss = 0.
      beta = 0e0    ! steepest descent

      if ( dot_product(gg, gg) .eq. 0 ) then
        call erexit('cgstep: grad vanishes identically')
      end if

      alfa = - dot_product( rr, rr) / dot_product( g, gg)

    else
      beta = dot_product(rr,rr)/dot_product(ra,ra)
      stat = Fop(F,F,dr,s)
      alfa = dot_product(rr,rr)/dot_product(s,dr)
    end if 

    ra = rr
    s = g + beta*s                        ! update solution step
    ss = alfa * dr                        ! update residual step
    x  =  x +  s                          ! update solution
    rr = rr + ss                          ! update residual
    forget = .false.
    cgstep1 = 0

  end function 

  subroutine cgstep1_close ( )
    if ( allocated( s)) then
      deallocate( s, ss, dr, ra)
    end if
  end subroutine 
end module 
