module cdfomel_mod
! 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
  n=n_in
  m=m_in
  niter=niter_in
  allocate(ssn(niter));ssn(1)=eps11
  allocate(s(m,niter),ss(n,niter))
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(:),rn		! solution,  residual
real       :: g(:), gg(:)		! direction, conjugate direction
real,optional      :: l(:), u(:)		! bounds
double precision :: ddot, alpha

liter = min0(iter,niter)
s(:,liter) = - g 				! initial direction
ss(:,liter) = - gg
write(0,*) "======================="
write(0,*) liter,sum(abs(gg)),sum(abs(ss(:,liter))),sum(abs(s(:,liter)))
write(0,*) "======================="
do j= 1, liter-1				! update direction
   alpha = ddot(n,gg,ss(:,j))/ssn(j)
   s(:,liter)  =  s(:,liter) + alpha *  s(:,j)
   ss(:,liter) = ss(:,liter) + alpha * ss(:,j)
end do
ssn(liter) = dmax1(ddot(n,ss(:,liter),ss(:,liter)),1.d-35)
if (liter == niter) then
   do j= 1, niter-1
      ssn(j) = ssn(j+1)
      s(:,j) =  s(:,j+1)
      ss(:,j) = ss(:,j+1)
   end do
end if

write(0,*) "======================="
write(0,*) sum(abs(ss(:,liter))),sum(abs(s(:,liter))),sum(abs(rr)),ssn(liter)
write(0,*) "======================="
		
alpha = ddot(n,ss(:,liter),rr)/ssn(liter)

if (present(l)) then
   call random_number(rn);rn =2*rn-1;write(0,*) "RN ",rn
   if (rn > 0) then
      rn = rn*1000
!      if (all(mask=(x - 4*alpha *  s(:,liter) >= l.or.x - 4*alpha *  s(:,liter)  <= u))) then
      do while (any(mask=(x - rn*alpha * s(:,liter) <= l.or.x - rn*alpha * s(:,liter) >= u)))
         rn=rn/2
      write(0,*) "1) Recomputing rn ", rn
      end do
      alpha=rn*alpha
      write(0,*) "1) Recomputing alpha ", alpha
   end if
!   if (any(mask=(x - alpha *  s(:,liter) <= l.or.x - alpha *  s(:,liter)  >= u))) then
!      alpha=0.5*alpha
!      do while (all(mask=(x - alpha * s(:,liter) >= l.or.x - alpha * s(:,liter) <= u)))
!         alpha=0.5*alpha
!         write(0,*) "2) Recomputing alpha ", alpha
!      end do
!   end if
end if
x  =  x - alpha *  s(:,liter)
!rr = rr - alpha * ss(:,liter)
end subroutine
!-----------------------------------------------------
double precision function ddot( nn, x, y)
integer          :: i, nn
real             :: x(nn), y(nn)
double precision :: val
val = 0.;
do i=1,nn
   val = val + x(i) * y(i)
end do
ddot = val;
end function
end module
