module bgmalgt

    use modminres, only: minres, lss_type
    use iso_c_binding, only: c_ptr
    
    implicit none
  
    interface
       subroutine fsub(n,x,f,ierr,pdata)
         use iso_c_binding, only: c_ptr
         implicit none
         integer, intent(in) :: n
         real(kind=8), intent(in) :: x(n)
         real(kind=8), intent(out) :: f 
         integer, intent(inout) :: ierr
         type(c_ptr), optional, intent(in) :: pdata
       end subroutine fsub
  
       subroutine gsub(n,x,g,ierr,pdata)
         use iso_c_binding, only: c_ptr
         implicit none
         integer, intent(in) :: n
         real(kind=8), intent(in) :: x(n)
         real(kind=8), intent(out) :: g(n)
         integer, intent(inout) :: ierr
         type(c_ptr), optional, intent(in) :: pdata
       end subroutine gsub
  
       subroutine hsub(n,x,lim,hnnz,hrow,hcol,hval,ierr,pdata)
         use iso_c_binding, only: c_ptr
         implicit none
         integer, intent(in) :: n,lim
         real(kind=8), intent(in) :: x(n)
         integer, intent(out) :: hnnz
         integer, intent(out) :: hrow(lim),hcol(lim)
         real(kind=8), intent(out) :: hval(lim)
         integer, intent(inout) :: ierr
         type(c_ptr), optional, intent(in) :: pdata
       end subroutine hsub
  
       subroutine stpsub(n,x,gsupn,inhdefstp,stp,iter,ierr,pdata)
         use iso_c_binding, only: c_ptr
         implicit none
         integer, intent(in) :: n,iter
         real(kind=8), intent(in) :: gsupn
         real(kind=8), intent(in) :: x(n)
         logical, intent(out) :: inhdefstp,stp
         integer, intent(inout) :: ierr
         type(c_ptr), optional, intent(in) :: pdata
       end subroutine stpsub
  
       subroutine tabsub(n,f,gsupn,iter,pdata)
          use iso_c_binding, only: c_ptr
          implicit none
          integer, intent(in) :: iter,n
          real(kind=8), intent(in) :: f,gsupn
          type(c_ptr), optional, intent(in) :: pdata
        end subroutine tabsub
    end interface
  
    private
    
    ! SCALAR PARAMETERS
    integer, parameter :: dimmax = 1000000, texmax = 20
    real(kind=8), parameter :: lrgstreal = huge( 1.0d0 ), macheps = epsilon( 1.0d0 ), &
         macheps12 = sqrt( macheps ), macheps14 = sqrt( macheps12 ), beta = 0.5d0, &
         gamma = 1.0d-04, eta = 1.0d+04, lspgmin = 1.0d-16, lspgmax = 1.0d+16, r = 0.1d0, &
         tau1 = 0.1d0, tau2 = 0.9d0, sgmsmall = 1.0d-08, sgmmin = 1.0d-08, sgmmax = 1.0d+16, &
         hmin = 1.0d-08, hmax = 1.0d+08
  
    public algt, algtb, project
  
  contains
  
    subroutine algt(evalf,evalg,evalh,hnnzmax,n,x,lind, &
      lbnd,uind,ubnd,f,g,gpsupn,ftarget,eps,maxit,extallowed,iter, &
      ierr,istop,evalstp,tabline,pdata)
  
   implicit none
  
   ! PROCEDURE ARGUMENTS
   procedure(fsub) :: evalf
   procedure(gsub) :: evalg
   procedure(hsub) :: evalh
   procedure(stpsub), optional :: evalstp
   procedure(tabsub), optional :: tabline
  
   ! SCALAR ARGUMENTS
   logical, intent(in) :: extallowed
   integer, intent(in) :: hnnzmax,maxit,n
   integer, intent(out) :: ierr,istop,iter
   real(kind=8), intent(in) :: ftarget,eps
   real(kind=8), intent(out) :: f,gpsupn
   type(c_ptr), optional, intent(in) :: pdata
  
   ! ARRAY ARGUMENTS
   logical, intent(in) :: lind(n),uind(n)
   real(kind=8), intent(in) :: lbnd(n),ubnd(n)
   real(kind=8), intent(inout) :: x(n)
   real(kind=8), intent(out) :: g(n)
  
   ! LOCAL SCALARS
   integer :: allocerr
   type(lss_type) :: lss

    ierr = 0
    lss%nnzmax = hnnzmax
    allocate(lss%matrix%row(hnnzmax),lss%matrix%col(hnnzmax),lss%matrix%val(hnnzmax),stat=allocerr)
    if ( allocerr .ne. 0 ) then
      write(*,*) 'ERROR: In algt, allocation error.'
      stop
    end if
    
    call algtb(evalf,evalg,evalh,n,x,lind,lbnd,uind,ubnd, &
            f,g,gpsupn,lss,ftarget,eps,maxit,extallowed,iter,ierr,istop, &
            evalstp,tabline,pdata)
    
    deallocate(lss%matrix%row,lss%matrix%col,lss%matrix%val,stat=allocerr)
      if ( allocerr .ne. 0 ) then
        write(*,*) 'ERROR: In algt, allocation error.'
        stop
    end if 
  end subroutine algt
  
  ! *****************************************************************
  ! *****************************************************************
  
  subroutine algtb(evalf,evalg,evalh,n,x,lind,lbnd,uind,ubnd, &
      f,g,gpsupn,lss,ftarget,eps,maxit,extallowed,iter,ierr,istop,evalstp, &
      tabline,pdata)
  
   implicit none
  
   ! PROCEDURE ARGUMENTS
   procedure(fsub) :: evalf
   procedure(gsub) :: evalg
   procedure(hsub) :: evalh
   procedure(stpsub), optional :: evalstp
   procedure(tabsub), optional :: tabline
  
   ! SCALAR ARGUMENTS
   logical, intent(in) :: extallowed
   integer, intent(in) :: maxit,n
   integer, intent(inout) :: ierr
   integer, intent(out) :: istop,iter
   real(kind=8), intent(in) :: ftarget,eps
   real(kind=8), intent(out) :: f,gpsupn
   type(c_ptr), optional, intent(in) :: pdata
  
   ! ARRAY ARGUMENTS
   logical, intent(in) :: lind(n),uind(n)
   real(kind=8), intent(in) :: lbnd(n),ubnd(n)
   real(kind=8), intent(inout) :: x(n)
   real(kind=8), intent(out) :: g(n)
   type(lss_type), intent(inout) :: lss
  
   ! Solves the box-constrained minimization problem
   !
   !        Minimize f(x) subject to l <= x <= u
   ! istop =  0 means gsupn <= eps
   ! istop =  6 means f <= ftarget
   ! istop =  8 means lack of progress (f did not improve in 3 consecuitove iterations)
   ! istop = 10 means iter >= maxit
   ! istop = 20 means user-provided stopping criterion satisfied
   
   ! LOCAL SCALARS
   logical :: extinhibited,extrecom,infrontier,inhdefstp,newface,sgminiundef,ustp,useregcub,calcgtrial
   integer :: extfails,i,k,navoidcy,nfnoprogr,nfr,tex,iterregcub,istopregcub,sigma,iternmr,iterrc,dim,&
        nkeepel
   real(kind=8) :: eps12,eps14,eps18,fbest,fext,fref,ftrial,gisupn,gtd,etaini,etamult, &
        gpsupnb,t,ttmp,eta,rho,rhomin,alpha,tau,amax,tref,a1,a2
  
   ! LOCAL ARRAYS
   logical :: frvar(n),fvpre(n)
   integer :: countg(4),ifr(n),renumbered(dimmax),idiag(n)
   real(kind=8) :: d(n),gp(n),gprev(n),gptrial(n),gtrial(n),&
        xext(n),xprev(n),xref(n),xtrial(n),resmin(n),d2(n),origv(n),diagv(n)

   ! ==========================================================================
   ! Initialization
   ! ==========================================================================
  
   iter = 0
   iternmr = 0
   iterrc = 0
  
   extfails = 0
   extinhibited = .false.
   
   navoidcy = 0
   nfnoprogr = 0
   fbest = lrgstreal
   gpsupnb = lrgstreal
  
   countg(1:4) = 0
   
   eps12 = sqrt( eps )
   eps14 = sqrt( eps12 )
   eps18 = sqrt( eps14 )
   a1 = 1.0d+8 
   a2 = 1.0d-16
   etaini  = 1.0d-8
   etamult = 0.9d0

   rhomin  = 1.0d-6
   alpha   = 1.0d-8
   tau     = 10

   eta = etaini
   rho = 0.0d0 

   sgminiundef = .true.
   useregcub   = .false.
  
   call project(n,lind,lbnd,uind,ubnd,x)
  
   call evalf(n,x,f,ierr,pdata)
   if ( ierr .ne. 0 ) return
  
   call evalg(n,x,g,ierr,pdata)
   if ( ierr .ne. 0 ) return


   ! ==========================================================================
   ! Main loop
   ! ==========================================================================
  
  100 continue  
   ! Check stopping criterion

    gp(1:n) = x(1:n) - g(1:n)
    
    call project(n,lind,lbnd,uind,ubnd,gp)

    gp(1:n) = gp(1:n) - x(1:n)

    gpsupn = maxval( abs( gp(1:n) ) )

    write(*,*) 'iter = ',iter,' f = ',f,' gpsupn = ',gpsupn
  
     if ( present( tabline ) ) then
         call tabline(n,f,gpsupn,iter,pdata)
      if ( ierr .ne. 0 ) return
    end if
  
    inhdefstp = .false.
   
    if ( present( evalstp ) ) then
      call evalstp(n,x,gpsupn,inhdefstp,ustp,iter,ierr,pdata)
      if ( ierr .ne. 0 ) return
  
      if ( ustp ) then
         write(*,*) 'User-defined stopping criterion satisfied'
         istop = 20
         return
      end if
    end if
  
    if ( .not. inhdefstp ) then
      if ( gpsupn .le. eps ) then
         write(*,*) 'Stopping criterion gpsupn <= eps satisfied'
         istop = 0
         return
      end if
      
      if ( f .le. ftarget ) then
         write(*,*) 'Stopping criterion f <= ftarget satisfied'
         istop = 6
         return
      end if
  
      if ( f .lt. fbest ) then
         nfnoprogr = 0
      
      else
         nfnoprogr = nfnoprogr + 1
         if ( nfnoprogr .gt. 3 ) then
            write(*,*) 'Lack of progress.'
            istop = 8
            return
         end if
      end if
      
      if ( f .lt. fbest .or. ( f .eq. fbest .and. gpsupn .lt. gpsupnb ) ) then
         fbest = f
         gpsupnb = gpsupn
      end if
  
      if ( iter .ge. maxit ) then
         write(*,*) 'Maximum of iterations reached'
         istop = 10
         return
      end if
    end if
   
    ! opt by an inner-to-the-face or a leaving-face iteration
    
    frvar(1:n) = ( .not. lind(1:n) .or. ( x(1:n) - lbnd(1:n) ) / max( 1.0d0, lbnd(1:n) ) .gt. 1.0d-12 ) .and. &
                    ( .not. uind(1:n) .or. ( ubnd(1:n) - x(1:n) ) / max( 1.0d0, ubnd(1:n) ) .gt. 1.0d-12 )
    
    newface = .false.
    if ( iter .eq. 0 .or. any( frvar(1:n) .neqv. fvpre(1:n) .or. &
            ( .not. frvar(1:n) .and. x(1:n) .ne. xprev(1:n) ) ) ) then
        newface = .true.
    end if
    
    write(*,*) 'Value of newface = ',newface
    
    ! iterate
    
    iter = iter + 1
    
    extrecom = .false.
    
    xprev(1:n) = x(1:n)
    gprev(1:n) = g(1:n)
    fvpre(1:n) = frvar(1:n)
    
    gisupn = max( 0.0d0, maxval( abs( gp(1:n) ), frvar(1:n) ) )
    
    write(*,*) 'gisupn = ',gisupn
    
    nfr = count( frvar(1:n) )
    
    write(*,*) 'Number of free variables: ',nfr
   
    call evalh(n,x,lss%nnzmax,lss%matrix%ne,lss%matrix%row,lss%matrix%col,lss%matrix%val,ierr,pdata)
    if ( ierr .ne. 0 ) return

    dim = max( n, maxval(lss%matrix%row(1:lss%matrix%ne)), maxval(lss%matrix%col(1:lss%matrix%ne)) )
    if ( dim .gt. dimmax )  then
        write(*,*) 'Increase dimmax in bgmalgt to at least ',dim,' and re-run.'
        stop
    end if

    lss%matrix%n = dim - n + nfr

    renumbered(1:n) = unpack( (/ (i,i=1,nfr) /), frvar(1:n), 0 )
    renumbered(n+1:dim) = (/ (nfr+i,i=1,dim-n) /)

    where ( renumbered(lss%matrix%row(1:lss%matrix%ne)) .ne. 0 .and. &
            renumbered(lss%matrix%col(1:lss%matrix%ne)) .ne. 0 )
        lss%matrix%row(1:lss%matrix%ne) = renumbered(lss%matrix%row(1:lss%matrix%ne))
        lss%matrix%col(1:lss%matrix%ne) = renumbered(lss%matrix%col(1:lss%matrix%ne))
    elsewhere
        lss%matrix%row(1:lss%matrix%ne) = 0
        lss%matrix%col(1:lss%matrix%ne) = 0
    end where

    nkeepel = count( lss%matrix%row(1:lss%matrix%ne) .ne. 0 )

    lss%matrix%val(1:nkeepel) = pack( lss%matrix%val(1:lss%matrix%ne), lss%matrix%row(1:lss%matrix%ne) .ne. 0 )
    
    lss%matrix%row(1:nkeepel) = pack( lss%matrix%row(1:lss%matrix%ne), lss%matrix%row(1:lss%matrix%ne) .ne. 0 )
    lss%matrix%col(1:nkeepel) = pack( lss%matrix%col(1:lss%matrix%ne), lss%matrix%col(1:lss%matrix%ne) .ne. 0 )

    lss%matrix%ne = nkeepel
    
    ! Set pointers to the diagonal elements

    idiag(1:nfr) = 0
    
    do k = 1,lss%matrix%ne
        if ( lss%matrix%row(k) .eq. lss%matrix%col(k) .and. lss%matrix%row(k) .le. nfr ) then
            if ( idiag(lss%matrix%row(k)) .eq. 0 ) then
                idiag(lss%matrix%row(k)) = k
                origv(lss%matrix%row(k)) = lss%matrix%val(k)
                diagv(lss%matrix%row(k)) = lss%matrix%val(k)
            else
                diagv(lss%matrix%row(k)) = diagv(lss%matrix%row(k)) + lss%matrix%val(k)
            end if
        end if
    end do

    do i = 1,nfr
        if ( idiag(i) .eq. 0 ) then
            if ( lss%matrix%ne + 1 .gt. lss%nnzmax ) then
                write(*,*) 'In bmalgt: lack of memory!'
                ierr = -98
                return
            end if
        
        lss%matrix%ne = lss%matrix%ne + 1
        
        lss%matrix%row(lss%matrix%ne) = i
        lss%matrix%col(lss%matrix%ne) = i
        lss%matrix%val(lss%matrix%ne) = 0.0d0
        
        idiag(i) = lss%matrix%ne
        origv(i) = 0.0d0
        diagv(i) = 0.0d0
        end if
    end do

    if ( gisupn .ge. r * gpsupn .and. .not. useregcub ) then
      iternmr = iternmr + 1
  
      ifr(1:nfr) = pack( (/ (i, i=1,n) /), frvar(1:n) )
  
      write(*,*) 'An inner-to-the-face iteration will be done'
  
      call newtonls()
      if ( ierr .ne. 0 ) return   
    else  
      write(*,*) 'A leaving-face iteration will be done'
      iterrc = iterrc + 1
      sigma = -1
      useregcub = .false.

400  continue

      d(1:n) = 0.0d0
      call gradprojtomodquad(n,g(1:n),lind(1:n),lbnd(1:n) - x(1:n),uind(1:n),ubnd(1:n) - x(1:n),rho,d(1:n), &
      lss,ftarget,eps,maxit,extallowed,iterregcub,istopregcub)

      write(*,*) 'dsupn = ',maxval( abs( d(1:n) ) ),' deucn = ',norm2( d(1:n) )

      xtrial(1:n) = x(1:n) + d(1:n)
      call evalf(n,xtrial,ftrial,ierr,pdata)
      if ( ierr .ne. 0 ) return

      write(*,*) 'ftrial =', ftrial, ' f = ', f, 'gtd = ', dot_product(g,d)

      if ( ftrial - f .gt. - alpha * norm2(d(1:n))**3 ) then 
        rho = max(rhomin,tau*rho)
        go to 400
      else
        rho = max(rhomin, 0.5d0*rho) 
      end if
      gtd = dot_product( g(1:n), d(1:n) )

      call evalg(n,xtrial,gtrial,ierr,pdata)
      if ( ierr .ne. 0 ) return

      f = ftrial
      x(1:n) = xtrial(1:n)
      g(1:n) = gtrial(1:n)
      go to 100
    end if
  
    if ( .not. calcgtrial) call evalg(n,xtrial,gtrial,ierr,pdata)
    if ( ierr .ne. 0 ) return
    
    if ( extallowed .and. .not. extinhibited .and. extrecom ) then
      gptrial(1:n) = xtrial(1:n) - gtrial(1:n)
      call project(n,lind,lbnd,uind,ubnd,gptrial)
      gptrial(1:n) = gptrial(1:n) - xtrial(1:n)

      if ( .not. maxval( abs( gptrial(1:n) ) ) .le. eps ) then
        if ( infrontier .or. dot_product( gtrial(ifr(1:nfr)), d(1:nfr) ) .le. beta * gtd ) then
          fref = ftrial
          xref(1:n) = xtrial(1:n)
          xext(1:n) = x(1:n)
 
          tex = 1
320 continue
          if ( tex .le. texmax ) then
            if ( .not. fref .le. ftarget ) then
              xext(ifr(1:nfr)) = x(ifr(1:nfr)) + ( 2.0d0 ** tex )* tref * d(1:nfr)
              call project(n,lind,lbnd,uind,ubnd,xext)

              call evalf(n,xext,fext,ierr,pdata)
              if ( ierr .ne. 0 ) return

              if ( fext .lt. fref ) then
                if ( .not. infrontier .and. ( ( 2.0d0 ** tex )* tref ) .gt. amax ) then 
                  infrontier = .true.    
                  sigma = 2
                end if
                fref = fext
                xref(1:n) = xext(1:n)
                write(*,*) 'Extrapolation with 2^tex = ',2.0d0 ** tex,' gave descent fext = ',fext
                tex = tex + 1
                go to 320
              else
                write(*,*) 'Extrapolation with 2^tex = ',2.0d0 ** tex,' failed.'
                if ( tex .gt. 1 )  then
                  extfails = 0
                else
                  extfails = extfails + 1
                  if ( extfails .eq. 3 ) then
                    extinhibited = .true.
                  end if
                end if
              end if
            else
              write(*,*) 'fref <= ftarget'
            end if
          else
            write(*,*) 'tex > texmax (tex =',tex,', texmax = ',texmax,')'
          end if

          if ( fref .lt. ftrial ) then
            ftrial = fref
            xtrial(1:n) = xref(1:n)
             
            call evalg(n,xtrial,gtrial,ierr,pdata)
            if ( ierr .ne. 0 ) return
          end if
        end if
      end if
    end if

   f = ftrial
   x(1:n) = xtrial(1:n)
   g(1:n) = gtrial(1:n)
     
   go to 100

  contains
   
   ! *****************************************************************
   ! *****************************************************************
   
  subroutine newtonls()
  
    implicit none
    
    ! LOCAL SCALARS
    logical :: armijo,flagmr
    integer :: i,j,z,rbdnnz,mriter,maxitmr,Dtype
    real(kind=8) :: beta1,beta2,normd1,famax,sHs,rHr,gtr,gts,rtr,sts
    
    ! LOCAL ARRAYS
    integer :: rbdind(nfr)
    real(kind=8) :: amaxv(nfr),rhs(dimmax),xamax(n),s(n),Hr(nfr),Hs(nfr)
    character :: amaxt(nfr),rbdtype(nfr)
 
    maxitmr = 2*nfr
    calcgtrial = .false.
    ! Compute search direction

    rhs(1:nfr) = - g(ifr(1:nfr))
    rhs(nfr+1:lss%matrix%n) = 0.0d0

    write (*,*)
    write (*,*) "minreseps   = ", eta
    write (*,*) "maxitminres = ", maxitmr

    call minres (dim,nfr,rhs(1:dim),eta,s(1:dim),resmin(1:dim),mriter,maxitmr,1,Dtype,lss)
    if ( Dtype .eq. 2 ) then ! MINRES reach maximum iterations.
        d(1:nfr) = s(1:nfr)
        flagmr = .true. 
    else
        flagmr = .false.
        ! calc Hs and Hr
        Hs(1:nfr) = 0.0d0 
        Hr(1:nfr) = 0.0d0
        do z = 1,lss%matrix%ne
            i = lss%matrix%row(z)
            j = lss%matrix%col(z)
            if (j .ne. i) then
                Hs(i) = Hs(i) + lss%matrix%val(z) * s(j)
                Hs(j) = Hs(j) + lss%matrix%val(z) * s(i)
                Hr(i) = Hr(i) + lss%matrix%val(z) * resmin(j)
                Hr(j) = Hr(j) + lss%matrix%val(z) * resmin(i)
            else
                Hs(i) = Hs(i) + lss%matrix%val(z) * s(j)
                Hr(i) = Hr(i) + lss%matrix%val(z) * resmin(j)
            end if 
        end do 
        sHs = dot_product( s(1:nfr)     , Hs(1:nfr)     ) 
        rHr = dot_product( resmin(1:nfr), Hr(1:nfr)     ) 
        gts = dot_product( g(ifr(1:nfr)), s(1:nfr)      )
        gtr = dot_product( g(ifr(1:nfr)), resmin(1:nfr) )
        sts = dot_product( s(1:nfr)     , s(1:nfr)      )
        rtr = dot_product( resmin(1:nfr), resmin(1:nfr) )

        if ( Dtype .eq. 0 ) then ! MINRES stop with SOL
            d(1:nfr) = s(1:nfr)
            if ( .not. ( gts + sHs .le. 0.0d0 .or. ( gts + sHs .le. 10.0d0*macheps .and. gts .lt. 0.0d0 ) ) &
              .or. norm2( Hr(1:nfr) ) .gt. eta*norm2( Hs(1:nfr) ) .or. min( sHs/sts, rHr/rtr )  .lt. eta ) then
                flagmr = .true.
            end if
        else
            d(1:nfr) = resmin(1:nfr)
            if ( .not. ( abs( gtr + rtr ) .le. macheps .and. gtr .lt. 0.0d0 ) .or. rHr .gt. 0.0d0 .or. &
            norm2( Hr(1:nfr) ) .le. eta*norm2( Hs(1:nfr) )) then 
                flagmr = .true. 
            end if
        end if
    end if
    write(*,*) 'flag = ', flagmr

    if ( flagmr ) then 
        eta = max(eta * etamult, eps)
        normd1 = norm2( d(1:nfr) )
        if ( normd1 .gt. a1*gpsupn ) then
            beta1 = ( a1 * gpsupn ) / normd1
            d2(1:nfr) = beta1 * d(1:nfr)
        else
            d2(1:nfr) = d(1:nfr)
        end if
        gtd = dot_product( g(ifr(1:nfr)), d2(1:nfr) )
        if ( gtd .gt. -a2 * ( gpsupn**2 ) ) then 
            beta2 = ( 1.0d0 - a2 ) / ( 1.0d0 + ( gtd / ( gpsupn**2 ) ) )
            d(1:nfr) = beta2 * d2(1:nfr) + (1.0d0 - beta2) * ( -g(ifr(1:nfr)) )
        else
            d(1:nfr) = d2(1:nfr)
        end if
    end if
    gtd = dot_product( g(ifr(1:nfr)), d(1:nfr) )
    write(*,*) 'dsupn = ',maxval( abs( d(1:nfr) ) ),' deucn = ',norm2( d(1:nfr) )

    ! Compute maximum feasible step
 
    amaxv(1:nfr) = lrgstreal
    
    where ( d(1:nfr) .lt. 0.0d0 .and. lind(ifr(1:nfr)) )
       amaxt(1:nfr) = 'L'
       amaxv(1:nfr) = ( lbnd(ifr(1:nfr)) - x(ifr(1:nfr)) ) / d(1:nfr)
    end where
 
    where ( d(1:nfr) .gt. 0.0d0 .and. uind(ifr(1:nfr)) )
       amaxt(1:nfr) = 'U'
       amaxv(1:nfr) = ( ubnd(ifr(1:nfr)) - x(ifr(1:nfr)) ) / d(1:nfr)
    end where
 
    amax = minval( amaxv(1:nfr) )
 
    write(*,*) 'amax = ',amax
 
    if ( amax .lt. lrgstreal ) then
       rbdnnz = count( amaxv(1:nfr) .eq. amax )
       rbdind(1:rbdnnz) = pack( (/ (i,i=1,nfr) /), amaxv(1:nfr) .eq. amax )
       rbdtype(1:rbdnnz) = pack( amaxt(1:nfr), amaxv(1:nfr) .eq. amax )
       write(*,*) 'rbdnnz = ',rbdnnz,' rbdind = ',rbdind(1:rbdnnz),' rbdtype = ',rbdtype(1:rbdnnz)
    else
       rbdnnz = 0
       write(*,*) 'rbdnnz = ',rbdnnz
    end if
    
    ! Test projected unitary step
    if ( amax .lt. 1.0d0 ) then
       write(*,*) 'Unitary step is infeasible. Trying projected unitary step.'
       
       xtrial(1:n) = x(1:n)
       xtrial(ifr(1:nfr)) = xtrial(ifr(1:nfr)) + d(1:nfr)
 
       call projinface(n,lind,lbnd,uind,ubnd,xtrial,nfr,ifr)
      
       call evalf(n,xtrial,ftrial,ierr,pdata)
       if ( ierr .ne. 0 ) return
  
       write(*,*) 'ftrial = ',ftrial
 
       if ( ftrial .le. ftarget ) then
          write(*,*) 'ftrial <= ftarget'
          return
 
       else if ( ftrial .le. f ) then
          write(*,*) 'Simple decrease was obtained and since we are in the frontier, this is enough.'
          extrecom   = .true.
          infrontier = .true.
          sigma = 2
          tref = 1.0d0
          return

       else if ( ftrial .le. f + macheps14 * max( 1.0d0, abs( f ) ) .and. &
           norm2( xtrial(ifr(1:nfr)) - x(ifr(1:nfr)) ) .le. macheps14 * norm2( x(ifr(1:nfr)) ) ) then
         if ( navoidcy .lt. 2 ) then
            navoidcy = navoidcy + 1
            useregcub = .false.
            write(*,*) 'The unitary step is infeasible and its projection did NOT give simple decrease.'
            write(*,*) 'However, the Euclidian norm of the step is small and the change in f is negligible.'
            write(*,*) 'Therefore, we accept it anyway.'
            return
         end if
      end if
    end if

    ! Backtracking
  
    write(*,*) 'gtd = ',gtd
    
    t = min( 1.0d0, amax )
    armijo = .false.
 
 310   continue
    
    xtrial(1:n) = x(1:n)
    xtrial(ifr(1:nfr)) = xtrial(ifr(1:nfr)) + t * d(1:nfr)
 
    if ( t .eq. amax ) then
       where ( rbdtype(1:rbdnnz) .eq. 'L' )
          xtrial(ifr(rbdind(1:rbdnnz))) = lbnd(ifr(rbdind(1:rbdnnz)))
       elsewhere
          xtrial(ifr(rbdind(1:rbdnnz))) = ubnd(ifr(rbdind(1:rbdnnz)))
       end where
    end if
    
    call evalf(n,xtrial,ftrial,ierr,pdata)
    if ( ierr .ne. 0 ) return
  
    write(*,*) 'in backtracking, t = ',t,' ftrial = ',ftrial
    
     if ( t .eq. 0.0d0 ) then
        write(*,*) 'null step in backtracking'
        xtrial(1:n) = x(1:n)
        ftrial = f
        return
        
    else if ( ftrial .le. ftarget ) then
       write(*,*) 'ftrial <= ftarget'
       return
 
    else if ( t .eq. amax .and. ftrial .le. f ) then
       write(*,*) 'Simple decrease was obtained and since we are in the frontier, this is enough.'
       extrecom   = .true.
       infrontier = .true.
       sigma  = 2
       tref = amax
       return
 
    else if ( ftrial .le. f + t * gamma * gtd ) then
       write(*,*) 'Armijo holds'
       armijo = .true.
       if ( t .eq. min( 1.0d0, amax ) ) then
          tref = t ! = min(1.0d0, amax)
          if ( flagmr ) then ! Step 5
            extrecom   = .true.
            if ( t .eq. amax ) then 
               infrontier = .true.
               sigma = 2
            else ! t eq. 1.0d0
               infrontier = .false.
               sigma = 1
            end if
            return ! step 5
          end if
       else 
         extrecom = .false.
         infrontier = .false.
         sigma = merge(1,0,flagmr) ! Set sigma to 1 if flagmr is true, otherwise set it to 0
         return ! step 4
       end if

    else if ( t .eq. min( 1.0d0, amax ) .and. &
       ftrial .le. f + macheps14 * max( 1.0d0, abs( f ) ) .and. &
       norm2( xtrial(ifr(1:nfr)) - x(ifr(1:nfr)) ) .le. macheps14 * norm2( x(ifr(1:nfr)) ) ) then
     if ( navoidcy .lt. 2 ) then
           navoidcy = navoidcy + 1
           write(*,*) 'The unitary or the maximum step did NOT give sufficient decrease.'
           write(*,*) 'However, the Euclidian norm of the step is small and the change in f is negligible.'
           write(*,*) 'Therefore, we accept it anyway.'
           return
       end if
    end if

    if ( .not. armijo ) then
        ttmp = ( - gtd * t ** 2 ) / ( 2.0d0 * ( ftrial - f - t * gtd ) )
        
        if ( .not. ( tau1 * t .le. ttmp .and. ttmp .le. tau2 * t ) ) then
        t = 0.5d0 * t
        else
        t = ttmp
        end if
    
        go to 310
    end if

    ! if you got here armijo .eq. true .and. flagmr .eq. 0 .and. t .eq. min(1.0d0,amax)
    ! Step 6, 6.1. 6.2 e 6.3
    sigma = merge(1,0,flagmr) ! Set sigma to 1 if flagmr is true, otherwise set it to 0

    call evalg(n,xtrial,gtrial,ierr,pdata)
    if ( ierr .ne. 0 ) return
    calcgtrial = .true.

    gptrial(1:n) = xtrial(1:n) - gtrial(1:n)
    call project(n,lind,lbnd,uind,ubnd,gptrial)
    gptrial(1:n) = gptrial(1:n) - xtrial(1:n)

    if ( .not. maxval( abs( gptrial(1:n) ) ) .le. eps ) then 
          fref = ftrial
          xref(1:n) = xtrial(1:n)
          xext(1:n) = x(1:n)
 
          tex = 1
810     if ( .not. fref .le. ftarget ) then
            t = ( 2.0d0 ** tex )* tref                        
            if ( t .gt. amax ) then 
                write(*,*) 'In Extrapolation withou Proj: Not in Omega'
                xamax(1:n) = x(1:n)
                xamax(ifr(1:nfr)) = x(ifr(1:nfr)) + amax * d(1:nfr)
                call evalf(n,xamax,famax,ierr,pdata)
                if ( ierr .ne. 0 ) return 
                write(*,*) 'famax = ', famax, f, fref
                if ( famax .le. f ) then 
                    ftrial = famax
                    xtrial(1:n) = xamax(1:n)
                    extrecom = .true.
                    infrontier = .true.
                    sigma = 2
                    calcgtrial = .false.
                    return
                else
                    ftrial = fref
                    xtrial(1:n) = xref(1:n)
                    sigma = 0
                    return 
                end if
            else 
                xext(ifr(1:nfr)) = x(ifr(1:nfr)) + t * d(1:nfr)
                call evalf(n,xext,fext,ierr,pdata)
                if ( ierr .ne. 0 ) return
                
                if ( fext .ge. f + t * gamma * gtd ) then
                    write(*,*) 'In Extrapolation withou Proj: Not Holds Armijo', fext,  f + t * gamma * gtd 
                    ftrial = fref
                    xtrial(1:n) = xref(1:n)
                    sigma = 0
                    return 
                else
                    write(*,*) 'Extrapolation without proj with 2^tex = ',2.0d0 ** tex,' gave descent fext = ',fext, tref
                    fref = fext
                    xref(1:n) = xext(1:n)
                    tex = tex + 1
                    calcgtrial = .false.
                    go to 810
                end if
            end if
        else
            write(*,*) 'fref <= ftarget'
            ftrial = fref
            xtrial(1:n) = xref(1:n)
            calcgtrial = .false.
            sigma = 0
            return 
        end if        
    end if
 
  end subroutine newtonls
  
  end subroutine algtb
  
  ! *****************************************************************
  ! *****************************************************************
  
  subroutine projinface(n,lind,lbnd,uind,ubnd,x,nind,ind)
  
   implicit none
  
   ! SCALAR ARGUMENTS
   integer, intent(in) :: n,nind
  
   ! ARRAY ARGUMENTS
   logical, intent(in) :: lind(n),uind(n)
   integer, intent(in) :: ind(nind)
   real(kind=8), intent(in) :: lbnd(n),ubnd(n)
   real(kind=8), intent(inout) :: x(n)
  
   where ( lind(ind(1:nind)) )
      x(ind(1:nind)) = max( lbnd(ind(1:nind)), x(ind(1:nind)) )
   end where
  
   where ( uind(ind(1:nind)) )
      x(ind(1:nind)) = min( ubnd(ind(1:nind)), x(ind(1:nind)) )
   end where
  
  end subroutine projinface
  
  ! *****************************************************************
  ! *****************************************************************
  
  subroutine project(n,lind,lbnd,uind,ubnd,x)
  
   implicit none
  
   ! SCALAR ARGUMENTS
   integer, intent(in) :: n
  
   ! ARRAY ARGUMENTS
   logical, intent(in) :: lind(n),uind(n)
   real(kind=8), intent(in) :: lbnd(n),ubnd(n)
   real(kind=8), intent(inout) :: x(n)
  
   where ( lind(1:n) )
      x(1:n) = max( lbnd(1:n), x(1:n) )
   end where
  
   where ( uind(1:n) )
      x(1:n) = min( ubnd(1:n), x(1:n) )
   end where
  
  end subroutine project

  ! *****************************************************************
  ! *****************************************************************

  subroutine evalfquadmod (lss,n,rho,gxk,s,fquadmod)
    implicit none 

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n 
    real(kind=8), intent(in) :: rho
    real(kind=8), intent(out) :: fquadmod

    ! ARRAY ARGUMENTS 
    real(kind=8), intent(in) :: gxk(n),s(n)
    type(lss_type), intent(in) :: lss

    ! LOCAL SCALARS
    integer :: i,j,k

    fquadmod = 0.0d0 
    do k = 1,lss%matrix%ne
        i = lss%matrix%row(k)
        j = lss%matrix%col(k)
        if (j .ne. i) then
            fquadmod = fquadmod + 2.0d0*s(i)*s(j)*lss%matrix%val(k)
        else
            fquadmod = fquadmod + s(i)*s(i)*lss%matrix%val(k)
        end if 
    end do   
    fquadmod = fquadmod + dot_product(gxk,s) + rho * ( norm2(s(1:n))**3 )
    

    end subroutine evalfquadmod


  ! *****************************************************************
  ! *****************************************************************

  subroutine evalgquadmod (lss,n,rho,gxk,s,gquadmod)
    implicit none 
    
    ! SCALAR ARGUMENTS
    integer, intent(in) :: n 
    real(kind=8), intent(in) :: rho
    
    ! ARRAY ARGUMENTS 
    real(kind=8), intent(in) ::gxk(n),s(n)
    real(kind=8), intent(out) :: gquadmod(n)
    type(lss_type), intent(in) :: lss
    
    ! LOCAL SCALARS
    integer :: i,j,k

    gquadmod(1:n) = 0.0d0 
    do k = 1,lss%matrix%ne
       i = lss%matrix%row(k)
       j = lss%matrix%col(k)
       if (j .ne. i) then
          gquadmod(i) = gquadmod(i) + lss%matrix%val(k) * s(j)
          gquadmod(j) = gquadmod(j) + lss%matrix%val(k) * s(i)
       else
          gquadmod(i) = gquadmod(i) + lss%matrix%val(k) * s(j)
       end if 
    end do   
    gquadmod(1:n) = gquadmod(1:n) + gxk(1:n) + ( 3.0d0 * rho * norm2(s(1:n)) ) * s(1:n)
    
    end subroutine evalgquadmod
    
  ! *****************************************************************
  ! *****************************************************************
  
  subroutine gradprojtomodquad(n,gxk,slind,slbnd,suind,subnd,rho,s, &
    lss,ftarget,eps,maxit,extallowed,iter,istop)

    implicit none

    ! SCALAR ARGUMENTS
    logical, intent(in) :: extallowed
    integer, intent(in) :: maxit,n
    integer, intent(out) :: istop,iter
    real(kind=8), intent(in) :: ftarget,eps,rho

    ! ARRAY ARGUMENTS
    logical, intent(in) :: slind(n),suind(n)
    real(kind=8), intent(in) :: gxk(n),slbnd(n),subnd(n)
    real(kind=8), intent(out) :: s(n)
    type(lss_type), intent(inout) :: lss

    ! OPTIMATION OF T2(xk,s) = g(xk)^t s + 1/2 s^t H(xk) s

    ! istop =  0 means gsupn <= eps or gsupn <= s^t s
    ! istop =  1 means gsupn <= eps^1/2 during 100 consecutive iterations
    ! istop =  2 means gsupn <= eps^1/4 during 5,000 consecutive iterations
    ! istop =  3 means gsupn <= eps^1/8 during 10,000 consecutive iterations
    ! istop =  6 means f <= ftarget
    ! istop =  8 means lack of progress (f did not improve in 3 consecuitove iterations)
    ! istop = 10 means iter >= maxit
    ! istop = 20 means user-provided stopping criterion satisfied
    
    ! LOCAL SCALARS
    logical :: extinhibited,extrecom,infrontier,inhdefstp,sgminiundef
    integer :: extfails,navoidcy,nfnoprogr,nfr,tex,inform
    real(kind=8) :: eps12,eps14,eps18,fbest,fext,fref,ftrial,gtd, &
        gpsupnb,t,ttmp,f,gpsupn

    ! LOCAL ARRAYS
    logical :: frvar(n),fvpre(n)
    integer :: countg(4),ifr(n)
    real(kind=8) :: g(n),d(n),gp(n),gprev(n),gptrial(n),gtrial(n), &
        sext(n),sprev(n),sref(n),strial(n)

    ! ==========================================================================
    ! Initialization
    ! ==========================================================================

    write(*,*) '--------------------------------------------------------------------'
    write(*,*) 'Try Cubic Regularization with rho = ', rho 
    write(*,*) '--------------------------------------------------------------------'

    iter = 0

    extfails = 0
    inform = 0
    extinhibited = .false.
    
    navoidcy = 0
    nfnoprogr = 0
    fbest = lrgstreal
    gpsupnb = lrgstreal

    countg(1:4) = 0
    
    eps12 = sqrt( eps )
    eps14 = sqrt( eps12 )
    eps18 = sqrt( eps14 )

    sgminiundef = .true.

    call project(n,slind,slbnd,suind,subnd,s)

    call evalfquadmod(lss,n,rho,gxk,s,f)
    
    call evalgquadmod (lss,n,rho,gxk,s,g)

 
    ! ==========================================================================
    ! Main loop
    ! ==========================================================================

    100 continue

    ! Check stopping criterion

    gp(1:n) = s(1:n) - g(1:n)

    call project(n,slind,slbnd,suind,subnd,gp)

    gp(1:n) = gp(1:n) - s(1:n)

    gpsupn = maxval( abs( gp(1:n) ) )  

    write(*,*) 'iter = ',iter,' fquadmod(xk,s) = ',f,' gquadmod(xk,s)psupn = ',gpsupn
    inhdefstp = .false.
    
    if ( gpsupn .le. eps ) then
        countg(1:4) = countg(1:4) + 1
    else
        countg(1:1) = 0
        if ( gpsupn .le. eps12 ) then
            countg(2:4) = countg(2:4) + 1
        else
            countg(2:2) = 0
            if ( gpsupn .le. eps14 ) then
                countg(3:4) = countg(3:4) + 1
            else
                countg(3:3) = 0
                if ( gpsupn .le. eps18 ) then
                    countg(4:4) = countg(4:4) + 1
                else
                    countg(4:4) = 0
                end if
            end if
        end if
    end if

    if ( .not. inhdefstp ) then

        ! stop criterion equation (32) with theta = 1
        if ( gpsupn .le. norm2( s(1:n ) )**2 ) then 
            write(*,*) 'Stopping criterion gpsupn <= theta*||s||^2'
            istop = 0 
            return 
        end if

        if ( countg(1) .ge. 1 ) then
            write(*,*) 'Stopping criterion gpsupn <= eps satisfied'
            istop = 0
            return
        else if ( countg(2) .ge. 100 ) then
            write(*,*) 'Stopping because gpsupn <= eps^1/2 during 100 consecutive iterations'
            istop = 1
            return
        else if ( countg(3) .ge. 5000 ) then
            write(*,*) 'Stopping because gpsupn <= eps^1/4 during 5,000 consecutive iterations'
            istop = 2
            return
        else if ( countg(4) .ge. 10000 ) then
            write(*,*) 'Stopping because gpsupn <= eps^1/8 during 10,000 consecutive iterations'
            istop = 3
            return
        end if
        
        if ( f .le. ftarget ) then
            write(*,*) 'Stopping criterion f <= ftarget satisfied'
            istop = 6
            return
        end if

        if ( f .lt. fbest ) then
             nfnoprogr = 0
        
        else
            nfnoprogr = nfnoprogr + 1
            if ( nfnoprogr .gt. 3 ) then
                write(*,*) 'Lack of progress.'
                istop = 8
                return
            end if
        end if
        
        if ( f .lt. fbest .or. ( f .eq. fbest .and. gpsupn .lt. gpsupnb ) ) then
            fbest = f
            gpsupnb = gpsupn
        end if

        if ( iter .ge. maxit ) then
            write(*,*) 'Maximum of iterations reached'
            istop = 10
            return
        end if
    end if
    
    ! iterate

    iter = iter + 1

    extrecom = .false.
    
    sprev(1:n) = s(1:n)
    gprev(1:n) = g(1:n)
    fvpre(1:n) = frvar(1:n)

    t = 1.0d0

    strial(1:n) = s(1:n) - g(1:n)
    call project(n,slind,slbnd,suind,subnd,strial)
    d(1:n) = strial(1:n) - s(1:n)

    gtd = dot_product( g(1:n), d(1:n) )

200 continue

    call evalfquadmod(lss,n,rho,gxk,strial,ftrial)

    if ( t .eq. 0.0d0 ) then
        write(*,*) 'null step in backtracking'
        strial(1:n) = s(1:n)
        ftrial = f
    
    else if ( ftrial .le. ftarget ) then
        write(*,*) 'ftrial <= ftarget'

    else if ( .not. ( ftrial .le. f + t * gamma * gtd ) ) then
        ttmp = ( - gtd * t ** 2 ) / ( 2.0d0 * ( ftrial - f - t * gtd ) )
    
        if ( .not. ( tau1 * t .le. ttmp .and. ttmp .le. tau2 * t ) ) then
            t = 0.5d0 * t
        else
            t = ttmp
        end if

        strial(1:n) = s(1:n) + t * d(1:n)

        go to 200
    end if

    call evalgquadmod(lss,n,rho,gxk,strial,gtrial)

    if ( extallowed .and. .not. extinhibited .and. extrecom ) then
        gptrial(1:n) = strial(1:n) - gtrial(1:n)
        call project(n,slind,slbnd,suind,subnd,gptrial)
        gptrial(1:n) = gptrial(1:n) - strial(1:n)

        if ( .not. maxval( abs( gptrial(1:n) ) ) .le. eps ) then
            if ( infrontier .or. dot_product( gtrial(ifr(1:nfr)), d(1:nfr) ) .le. beta * gtd ) then
                fref = ftrial
                sref(1:n) = strial(1:n)
        
                tex = 1
330 continue
                if ( tex .le. texmax ) then
                    if ( .not. fref .le. ftarget ) then
                        sext(1:n) = s(1:n) + 2.0d0 ** tex * ( strial(1:n) - s(1:n) )

                        call project(n,slind,slbnd,suind,subnd,sext)

                        call evalfquadmod(lss,n,rho,gxk,sext,fext)

                        if ( fext .lt. fref ) then
                            fref = fext
                            sref(1:n) = sext(1:n)
                            
                            write(*,*) 'Extrapolation with 2^tex = ',2.0d0 ** tex,' gave descent fext = ',fext
                            
                            tex = tex + 1
                            go to 330
                        else
                            write(*,*) 'Extrapolation with 2^tex = ',2.0d0 ** tex,' failed.'
                            if ( tex .gt. 1 )  then
                                extfails = 0
                            else
                                extfails = extfails + 1
                                if ( extfails .eq. 3 ) then
                                    extinhibited = .true.
                                end if
                            end if
                        end if
                    else
                        write(*,*) 'fref <= ftarget'
                    end if
                else
                    write(*,*) 'tex > texmax (tex =',tex,', texmax = ',texmax,')'
                end if

                if ( fref .lt. ftrial ) then
                    ftrial = fref
                    strial(1:n) = sref(1:n)
                    
                    call evalgquadmod (lss,n,rho,gxk,strial,gtrial)
                end if
            end if
        end if
    end if

    f = ftrial
    s(1:n) = strial(1:n)
    g(1:n) = gtrial(1:n)

    go to 100

 end subroutine gradprojtomodquad

  end module bgmalgt
  