module nl_solver
  use chain_mod
     use arnoldi_mod
  use llist_mod
  use ddot_mod
	 use hub
	use  ds

  implicit none
  logical, parameter, private  :: T = .true., F = .false.

contains

  subroutine nlcg_solver(oper,func,func_init,x, dat, niter, reg, nreg, eps &
  ,             x0,err,verb,wt,known)
    optional :: x0,err,verb,reg,nreg,eps,wt,known
   interface
       function oper (adj, add, x, dat) result (stat)
         integer              :: stat
         logical, intent (in) :: adj, add
         real, dimension (:)  :: x, dat
       end function oper
       function reg (adj, add, x, dat) result (stat)
         integer              :: stat
         logical, intent (in) :: adj, add
         real, dimension (:)  :: x, dat
       end function reg
			 subroutine func_init(data,wt,known,eps)
				real :: data(:)
				real, optional :: wt(:),eps
				logical, optional :: known(:)
			end subroutine
 				real  function func(oper,xval,grad,reg,nreg) result(val)
			    optional :: reg,nreg,grad
			    real :: val
			    real :: grad(:),xval(:)
			    interface
			      integer function oper(adj,add,model,data) result(stat)
			        integer :: stat
  			      logical :: adj ,add
      			  real    :: model(:),data(:)
						end function
			      integer function reg(adj,add,model,data) result(stat)
    			    integer :: stat
   			      logical :: adj ,add
      			  real    :: model(:),data(:)
						end function
					end interface
				end function
    end interface
    real, dimension (:),    intent (in)         :: dat, x0,wt
    real, dimension (:),    intent (out)        :: x, err !, res, resm
    real,                   intent (in)         :: eps
    integer,                intent (in)         :: niter, nreg
!    real, dimension (:,:),  intent (out)        :: xmov, rmov   
    logical,                intent (in)         :: verb
    logical, dimension (:), intent (in)         :: known

    real, dimension (size (x))                  :: g,xold,pdir,gnew
    real, dimension (:), pointer            :: rr, gg
    real, dimension (:), pointer                :: rd, rm, gd, gm, wht
    integer                                     :: iter, stat, nd,nr
    logical                                     :: forget, v,stuck
    integer                                     :: dprd, dprm, dpx, dpg,status
		integer                                     :: old_status
		real                                        :: tempr,new_eval,beta,epsilon
	
	

		if(present(nreg) .and. .not. present(reg)) call erexit("must supply&
      &reg and nreg")
		if(.not. present(nreg) .and. present(reg)) call erexit("must supply&
      &reg and nreg")

		if(present(nreg))  then
			nr=nreg
		else 
			nr=0
		end if
		

		if(present(eps)) then
			 epsilon=eps
		else 
			epsilon=1.
		end if

    nd = size (dat)
		allocate(rr(nd+nr))
    rd => rr (1 : nd) ; 
    gd => gg (1 : nd) ;
		if(nr>0) then
			rm => rr (1 + nd :)
 			gm => gg (1 + nd :)
		end if


    rd = dat
		if( nr .ne. 0) rm=0.
		if(present(wt)) then 
			if(present(known)) then 
				call func_init(rr,wt=wt,known=known,eps=epsilon)
			else
				call func_init(rr,wt=wt,eps=epsilon)
			end if
		else
			if(present(known)) then 
				call func_init(rr,known=known,eps=epsilon)
			else
				call func_init(rr,eps=epsilon)
			end if
		end if
		call ds_size(size(x),nd+nr)

    forget = F
    v = F ; if (present (verb)) v = verb
		xold=x;iter=1; stuck=F; old_status=0

		do while(iter<= niter .and.  .not. stuck)
			if(iter==1 .or. old_status .ne. 0) then
				if(nr > 0) then
					tempr= func(oper, x, g, reg, nr)
				else
					tempr= func(oper, x, g)
				end if
				pdir=-g
			end if
			if(nr >0) then
				status=ds_search_simple(oper,func,xold,x,pdir,gnew,new_eval,reg,nr)
			else
				status=ds_search_simple(oper,func,xold,x,pdir,gnew,new_eval)
			end if

			call srite("grad.H",gnew,size(gnew)*4)
			if(status==0) then  !Succesful CGstep
				  !	beta=dot_product(gnew- g,gnew)/sum(g**2)
					beta=dot_product(gnew,gnew)/sum(g**2)  !Fletcher Reeves
					if(v) write(0,*) "nl solve", iter,new_eval
					pdir=-gnew+ beta * pdir
					xold=x; g=gnew
					if(present(err)) err(iter)=new_eval
					iter=iter+1
			else if(status==1)  then !CG not going down, use steepest descent
				if(old_status .ne. 0) then !last time we had the same problem
					write(0,*) "Stuck at a minimum after ",iter," itterations"
					stuck=T
				else if(v) then
           write(0,*) "cgstep did not descend switching to steepest descent"
				end if
			else if(status==2)  then !can't descend fast enough try steepest descent
				if(old_status.ne.0) then  !last time same problem
					write(0,*) "can't descend fast enough try changing minimum length"
					stuck=T
				else if(v) then
					write(0,*) "cgstep did not descend fast enough switching to steepest"
				end if
			else if(status==3) then !the newton length ==0 
				if(old_status.ne.0) then  !last time same problem
					write(0,*) "newton length 0 for steepest, quitting after",&
             iter," itterations"
					stuck=T
				else if(v) then
					write(0,*) "objective function did not decrease switching to steepest descent"
				end if
			end if
			old_status=status
		end do

	deallocate(rr)
  end subroutine nlcg_solver
end module