module cgstep_mod
  implicit none
  real, dimension (:), allocatable, private :: s, ss
  contains
  integer function cgstep( forget, x, g, rr, gg,tm1,tm2)
    real, dimension (:) :: x, g, rr, gg
    logical :: forget
    integer :: i,tm1,tm2
    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 ( dprod_omp(gg, gg) .eq. 0 ) then
        call erexit('cgstep: grad vanishes identically')
      end if
      alfa = -  dprod_omp( gg, rr) /  dprod_omp( gg, gg)
    else
      gdg =  dprod_omp( gg, gg)
      ! search plane by solving 2-by-2
      sds = dprod_omp( ss, ss)
      ! G . (R - G*alfa - S*beta) = 0
      gds =  dprod_omp( 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 = - ( dprod_omp( gg, rr))
      sdr = - ( dprod_omp( ss, rr))
      alfa = ( sds * gdr - gds * sdr ) / determ
      beta = (-gds * gdr + gdg * sdr ) / determ
    end if
 !$OMP PARALLEL DO
    do i=1,size(s)
      s(i)=s(i)*beta+alfa*g(i)
      x(i)=x(i)+s(i)
    end do
 !$OMP END PARALLEL DO
 !$OMP PARALLEL DO
    do i=1,size(ss)
      ss(i)=ss(i)*beta+alfa*gg(i)
      rr(i)=rr(i)+ss(i)
    end do
 !$OMP END PARALLEL DO
    forget = .false.
    cgstep = 0
  end function
  subroutine cgstep_close ( )
    if ( allocated( s)) then
      deallocate( s, ss)
    end if
  end subroutine
  double precision function dprod_omp(vec1,vec2)
    real :: vec1(:),vec2(:)
    double precision mys(88),v1,v2
    integer :: i
    integer :: ith
    integer,external :: omp_get_thread_num


  mys=0
!!$OMP PARALLEL  default(shared) private(i,v1,v2,ith)
!  ith=omp_get_thread_num()+1
!!$OMP DO
ith=1
  do i=1,size(vec1)
     v1=vec1(i)
     v2=vec2(i)
     mys(ith)=mys(ith)+v1*v2
  end do
!!$OMP END DO
!!$OMP END PARALLEL 
  dprod_omp=sum(mys)
  end function
  
end module
