module cdfomel_mod1
! A step of conjugate-direction descent.
!
implicit none
integer ,private                        :: n, m, niter
real, private                           :: eps11
double precision,private, allocatable   :: ssn(:)
real,private,allocatable                :: s(:,:), ss(:,:)       ! step,conjugate step

contains
!--------------------------------------------------
subroutine cdplus_init(niter_in,n_in,m_in,eps_in)
integer :: n_in, m_in,niter_in
real    :: eps_in
  eps11=eps_in
    write(0,*) "cdstep ",eps11
  n=n_in
  m=m_in
  niter=niter_in
!  allocate(ssn(niter))
!  allocate(s(m,niter),ss(n,niter))
  allocate(ssn(2))
  allocate(s(n,2),ss(m,2))
end subroutine
!--------------------------------------------------
subroutine cdplus_close()
  deallocate(ssn,s,ss)
end subroutine
!--------------------------------------------------
subroutine cdplus( iter, x, g, rr, gg,l,u)
integer    :: i,j,liter,iter
real       :: x(:), rr(:)		! solution,  residual
real       :: g(:), gg(:)		! direction, conjugate direction
real,optional      :: l(:), u(:)		! bounds
double precision :: ddot, alpha

ssn(1)=eps11
liter = min0(iter,2)
s(:,liter) = - g 				! initial direction
ss(:,liter) = - gg
write(0,*) "01) alpha = ",liter,alpha,minval(gg),maxval(gg),minval(ss(:,liter)),maxval(ss(:,liter))
write(0,*) "01) alpha = ",liter,alpha,minval(g),maxval(g),minval(s(:,liter)),maxval(s(:,liter))
if (liter.eq.2) then
   do j= 1, liter-1				! update direction
      alpha = ddot(m,gg,ss(:,j))/eps11
      s(:,liter)  =  s(:,liter) + alpha *  s(:,j)
      ss(:,liter) = ss(:,liter) + alpha * ss(:,j)
   end do
else
   alpha = ddot(m,gg,ss(:,liter))/eps11
   s(:,liter)  =  s(:,liter) + alpha *  s(:,liter)
   ss(:,liter) = ss(:,liter) + alpha * ss(:,liter)
end if
write(0,*) "01) alpha = ",alpha,minval(gg),maxval(gg),minval(ss(:,liter)),maxval(ss(:,liter))
write(0,*) "01) alpha = ",alpha,minval(g),maxval(g),minval(s(:,liter)),maxval(s(:,liter))
ssn(liter) = dmax1(ddot(m,ss(:,liter),ss(:,liter)),1.d-35)
			
alpha = ddot(m,ss(:,liter),rr)/ssn(liter)
write(0,*) "02) alpha = ",alpha,maxval(ss(:,liter)),minval(ss(:,liter)),maxval(rr),minval(rr),ssn(liter)
if (present(l)) then
   if (all(mask=(x - 2*alpha *  s(:,liter) >= l.or.x - 2*alpha *  s(:,liter)  <= u))) then
      alpha=2*alpha
      write(0,*) "1) Recomputing alpha ", alpha
   end if
end if
x = x - alpha *  s(:,liter)
!rr = rr - alpha * ss(:,liter)
end subroutine
!-----------------------------------------------------
double precision function ddot( n, x, y)
integer          :: i, n
real             :: x(n), y(n)
double precision :: val
val = 0.;
do i=1,n
   val = val + x(i) * y(i)
end do
ddot = val;
end function
end module
