module bgmalgp

    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 algpunc, algpuncb, algp, algpb, project
  
  contains
    
    ! *****************************************************************
    ! *****************************************************************
  
    subroutine algpunc(evalf,evalg,evalh,hnnzmax,n,x,f,g,gsupn, &
         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,gsupn
      type(c_ptr), optional, intent(in) :: pdata
  
      ! ARRAY ARGUMENTS
      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 algp, allocation error.'
        stop
      end if

      call algpuncb(evalf,evalg,evalh,n,x,f,g,gsupn,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 algp, allocation error.'
        stop
      end if      
    end subroutine algpunc
  
    ! *****************************************************************
    ! *****************************************************************
  
    subroutine algpuncb(evalf,evalg,evalh,n,x,f,g,gsupn,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,gsupn
      type(c_ptr), optional, intent(in) :: pdata
  
      ! ARRAY ARGUMENTS
      real(kind=8), intent(inout) :: x(n)
      real(kind=8), intent(out) :: g(n)
      type(lss_type), intent(inout) :: lss
  
      ! 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,inhdefstp,sgminiundef,ustp
      integer :: dim,extfails,i,k,navoidcy,nfnoprogr,tex, &
      mriter,Dtype,maxitmr
      real(kind=8) :: eps12,eps14,eps18,fbest,fbtrial,fext,fref,ftrial, &
           gtd,gsupnb,sgmini,sgmtrial,t,ttmp,amreps,bmreps, &
           mreps,a1,a2,beta1,beta2,normd1
    
      ! LOCAL ARRAYS
      integer :: countg(4),idiag(n),ifr(n)
      real(kind=8) :: d(n),diagv(n),gtrial(n),origv(n),rhs(dimmax),xbtrial(n), &
           xext(n),xref(n),xtrial(n),resmin(n),d2(n)
  
      !PARAMETERS
      real(kind=8), parameter :: mrepsi     =   1.0d-01
      real(kind=8), parameter :: mrepsf     =   1.0d-08
    
      ! ==========================================================================
      ! Initialization
      ! ==========================================================================
  
      iter = 0
  
      extfails = 0
      extinhibited = .false.
  
      navoidcy = 0
      nfnoprogr = 0
      fbest = lrgstreal
      gsupnb = lrgstreal
  
      countg(1:4) = 0
      ifr(1:n) = (/ (i, i=1,n) /)
      
      eps12 = sqrt( eps )
      eps14 = sqrt( eps12 )
      eps18 = sqrt( eps14 )
  
      maxitmr = 2*n
      a1 = 1.0d+8 
      a2 = 1.0d-16
  
      sgminiundef = .true.
      
      call evalf(n,x,f,ierr,pdata)
      if ( ierr .ne. 0 ) return
  
      call evalg(n,x,g,ierr,pdata)
      if ( ierr .ne. 0 ) return
  
      gsupn = maxval( abs( g(1:n) ) )
      if ( gsupn .ne. 0.0d0 ) then
          amreps = log10( mrepsf / mrepsi ) / log10( sqrt(mrepsf) / gsupn )
          bmreps = log10( mrepsi ) - amreps * log10( gsupn )
       else
          amreps = 0.0d0
          bmreps = mrepsf
       end if
     
      ! ==========================================================================
      ! Main loop
      ! ==========================================================================
  100 continue
    
      ! Check stopping criterion
      if (iter .GT. 0) then 
        gsupn = maxval( abs( g(1:n) ) )
      end if
  
      write(*,*) 'iter = ',iter,' f = ',f,' gsupn = ',gsupn
 
      if ( present( tabline ) ) then
          call tabline(n,f,gsupn,iter,pdata)
          if ( ierr .ne. 0 ) return
       end if
  
      inhdefstp = .false.
      
      if ( present( evalstp ) ) then
         call evalstp(n,x,gsupn,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 ( gsupn .le. eps  ) then
            write(*,*) 'Stopping criterion gsupn <= 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. gsupn .lt. gsupnb ) ) then
            fbest = f
            gsupnb = gsupn
         end if
  
         if ( iter .ge. maxit ) then
            istop = 10
            return
         end if
      end if
  
      ! iterate
  
      iter = iter + 1
  
      extrecom = .false.
      
      fbtrial = lrgstreal  
      dim = n

  
      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 bgmalgp to at least ',dim,' and re-run.'
        stop
      end if
    
      lss%matrix%n = dim

      ! Set pointers to the diagonal elements

      idiag(1:n) = 0
        
      do k = 1,lss%matrix%ne
        if ( lss%matrix%row(k) .eq. lss%matrix%col(k) .and. lss%matrix%row(k) .le. n ) 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,n
        if ( idiag(i) .eq. 0 ) then
            if ( lss%matrix%ne + 1 .gt. lss%nnzmax ) then
            write(*,*) 'In algpuncb (algp.f90): 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        
  
      ! Compute search direction
  
      rhs(1:n) = - g(1:n)
      rhs(n+1:dim) = 0.0d0 
  
      mreps   = 10.0d0 ** ( amreps * log10( gsupn ) + bmreps )
      mreps   = max( mrepsf, min( mrepsi, mreps ) )
  
      write (*,*)
      write (*,*) "minreseps   = ", mreps
      write (*,*) "maxitminres = ", maxitmr
      call minres (dim,n,rhs(1:n),mreps,d(1:n),resmin(1:n),mriter,maxitmr,1,Dtype,lss)

      if (Dtype .eq. 1 .and. mriter .eq. 0) then 
         d(1:n) = resmin(1:n)
         write(*,*) 'Minus gradient accepted as direction.'
      else 
         normd1 = norm2( d(1:n) )
         if ( normd1 .gt. a1*gsupn ) then
            beta1 = ( a1 * gsupn ) / normd1
            d2(1:n) = beta1 * d(1:n)
         else
            d2(1:n) = d(1:n)
         end if
         
         gtd = dot_product( g(1:n), d2(1:n) )
         if ( gtd .gt. -a2 * ( gsupn**2 ) ) then 
            beta2 = ( 1 - a2 ) / ( 1 + ( gtd / ( gsupn**2 ) ) )
            d(1:n) = beta2 * d2(1:n) + ( 1.0d0 - beta2 ) * ( -g(1:n) )
         else
            d(1:n) = d2(1:n)
         end if
      end if 
      gtd = dot_product( g(1:n), d(1:n) )

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

      go to 300
     
      ! Backtracking
  300 continue

      write(*,*) 'gtd = ',gtd
        
      t = 1.0d0
  
  310 continue
        
      xtrial(1:n) = x(1:n) + t * d(1:n)
  
      call evalf(n,xtrial,ftrial,ierr,pdata)
      if ( ierr .ne. 0 ) return
      
      write(*,*) 'in backtracking, t = ',t,' ftrial = ',ftrial
  
      if ( ftrial .lt. fbtrial ) then
         fbtrial = ftrial
         xbtrial(1:n) = xtrial(1:n)
      end if
        
      if ( t .eq. 0.0d0 ) then
         write(*,*) 'null step in backtracking'
         xtrial(1:n) = x(1:n)
         ftrial = f
         go to 400
         
      else if ( ftrial .le. ftarget ) then
         write(*,*) 'ftrial <= ftarget'
         go to 400
         
      else if ( ftrial .le. f + t * gamma * gtd ) then
         write(*,*) 'Armijo holds'
         if ( t .eq. 1.0d0 ) then
            extrecom = .true.
         end if
         go to 400
  
      else if ( t .eq. 1.0d0 .and. sgmtrial .eq. 0.0d0 .and. &
           ftrial .le. f + macheps14 * max( 1.0d0, abs( f ) ) .and. &
           norm2( d(1:n) ) .le. macheps14 * norm2( x(1:n) ) ) then
         if ( navoidcy .lt. 2 ) then
            navoidcy = navoidcy + 1
            write(*,*) 'The Newton 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.'
            go to 400
         end if
      end if
  
      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
  
  400 continue
      
      if ( fbtrial .lt. ftrial ) then
         ftrial = fbtrial
         xtrial(1:n) = xbtrial(1:n)
      end if
         
      if ( sgmtrial .eq. 0.0d0 ) then
         sgmini = min( max( sgmmin, 0.5d0 * sgmini   ), sgmmax )
      else
         sgmini = min( max( sgmmin, 0.5d0 * sgmtrial ), sgmmax )
      end if
      
      call evalg(n,xtrial,gtrial,ierr,pdata)
      if ( ierr .ne. 0 ) return
  
      ! Extrapolate if recommended
      
      if ( extallowed .and. .not. extinhibited .and. extrecom ) then
         if ( .not. maxval( abs( gtrial(1:n) ) ) .le. eps ) then
            if ( dot_product( gtrial(1:n), d(1:n) ) .le. beta * gtd ) then     
               fref = ftrial
               xref(1:n) = xtrial(1:n)
      
               tex = 1
410 continue
               if ( tex .le. texmax ) then
                  if ( .not. fref .le. ftarget ) then
                     xext(1:n) = x(1:n) + 2.0d0 ** tex * d(1:n)

                     call evalf(n,xext,fext,ierr,pdata)
                     if ( ierr .ne. 0 ) return
                     
                     if ( fext .lt. fref ) then
                        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 410
                     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  
    end subroutine algpuncb
    
    ! *****************************************************************
    ! *****************************************************************
      
    subroutine algp(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 lssini, allocation error.'
      stop
    end if
    
    call algpb(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 algp, allocation error.'
        stop
    end if 
  end subroutine algp
  
  ! *****************************************************************
  ! *****************************************************************
  
  subroutine algpb(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
   integer :: extfails,i,navoidcy,nfnoprogr,nfr,tex
   real(kind=8) :: eps12,eps14,eps18,fbest,fbtrial,fext,fref,ftrial,gisupn,gtd, &
        gpsupnb,lspg,sts,sty,t,ttmp,amreps,bmreps,mreps,beta1,beta2,normd1,a1,a2,&
        xsupn,gsupn,tref
  
   ! LOCAL ARRAYS
   logical :: frvar(n),fvpre(n)
   integer :: countg(4),ifr(n)
   real(kind=8) :: d(n),gp(n),gprev(n),gptrial(n),gtrial(n),s(n),xbtrial(n), &
        xext(n),xprev(n),xref(n),xtrial(n),resmin(n),d2(n)

   !PARAMETERS
   real(kind=8), parameter :: mrepsi     =   1.0d-01
   real(kind=8), parameter :: mrepsf     =   1.0d-08
  
   ! ==========================================================================
   ! Initialization
   ! ==========================================================================
  
   iter = 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

   sgminiundef = .true.
  
   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

    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) ) )   

    if ( gpsupn .ne. 0.0d0 ) then
        amreps = log10( mrepsf / mrepsi ) / log10( sqrt(eps) / gpsupn )
        bmreps = log10( mrepsi ) - amreps * log10( gpsupn )
    else
        amreps = 0.0d0
        bmreps = mrepsf
    end if
  
   ! ==========================================================================
   ! Main loop
   ! ==========================================================================

  100 continue
  
   ! Check stopping criterion

    if(iter .gt. 0) then 
        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) ) )
    end if
    gsupn = maxval( abs( g(1:n) ) )
    xsupn = maxval( abs( x(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
  
   if ( gisupn .ge. r * gpsupn ) then
  
      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
   
      if ( fbtrial .lt. ftrial ) then
         ftrial = fbtrial
         xtrial(1:n) = xbtrial(1:n)
      end if
      
   else
  
      write(*,*) 'A leaving-face iteration will be done'
  
      if ( iter .eq. 1 .or. sty .le. 0.0d0 ) then
         lspg = max( 1.0d0, xsupn ) / gsupn
      else
         lspg = sts / sty
      end if
  
      lspg = max( lspgmin, min( lspg, lspgmax ) )
  
      t = 1.0d0
  
      xtrial(1:n) = x(1:n) - lspg * g(1:n)
      call project(n,lind,lbnd,uind,ubnd,xtrial)
  
      d(1:n) = xtrial(1:n) - x(1:n)
  
      gtd = dot_product( g(1:n), d(1:n) )
  
  200    continue
    
      call evalf(n,xtrial,ftrial,ierr,pdata)
      if ( ierr .ne. 0 ) return
  
      if ( t .eq. 0.0d0 ) then
         write(*,*) 'null step in backtracking'
         xtrial(1:n) = x(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
  
         xtrial(1:n) = x(1:n) + t * d(1:n)
  
         go to 200
      end if
  
   end if
  
   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
                     fref = fext
                     xref(1:n) = xext(1:n)
                     
                     write(*,*) 'Extrapolation with 2^tex = ',2.0d0 ** tex,' gave descent fext = ',fext, tref
                     
                     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)
   
   s(1:n) = x(1:n) - xprev(1:n)
   sty = dot_product( s(1:n), g(1:n) - gprev(1:n) )
   sts = sum( s(1:n) ** 2 )
  
   go to 100
  
  contains
   
   ! *****************************************************************
   ! *****************************************************************
   
   subroutine newtonls()
  
     implicit none
     
     ! LOCAL SCALARS
     integer :: dim,i,k,nkeepel,rbdnnz,mriter,maxitmr,Dtype
     real(kind=8) :: amax
     
     ! LOCAL ARRAYS
     integer :: idiag(nfr),rbdind(nfr),renumbered(dimmax)
     real(kind=8) :: amaxv(nfr),origv(nfr),rhs(dimmax),diagv(nfr)
     character :: amaxt(nfr),rbdtype(nfr)
  
     fbtrial = lrgstreal
     maxitmr = 2*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 bgmalgp 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 newtonls (algp.f90): 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
  
     ! Compute search direction

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

     mreps   = 10.0d0 ** ( amreps * log10( gpsupn ) + bmreps )
     mreps   = max( mrepsf, min( mrepsi, mreps ) )
     write (*,*)
     write (*,*) "minreseps   = ", mreps
     write (*,*) "maxitminres = ", maxitmr

     call minres (dim,nfr,rhs(1:dim),mreps,d(1:dim),resmin(1:dim),mriter,maxitmr,1,Dtype,lss)

     if (Dtype .eq. 1 .and. mriter .eq. 0) then 
       d(1:nfr) = resmin(1:nfr) ! residuo here is minus gradient, then d = d1
       write(*,*) 'Minus gradient accepted as direction.'
     else 
       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 - a2 ) / (1 + ( 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 .lt. fbtrial ) then
           fbtrial = ftrial
           xbtrial(1:n) = xtrial(1:n)
        end if
        
        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.
           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
             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 )
  
  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 ( ftrial .lt. fbtrial ) then
        fbtrial = ftrial
        xbtrial(1:n) = xtrial(1:n)
     end if
     
      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 ! Step 2.3
        write(*,*) 'Simple decrease was obtained and since we are in the frontier, this is enough.'
        extrecom   = .true.
        infrontier = .true.
        tref = amax
        return
  
     else if ( ftrial .le. f + t * gamma * gtd ) then
        write(*,*) 'Armijo holds'
        if ( t .eq. min( 1.0d0, amax ) ) then
           extrecom   = .true.
           infrontier = .false.
           tref = t ! = min(1.0d0, amax)
        end if
        return

     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
  
     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 subroutine newtonls
  
  end subroutine algpb
  
  
  ! *****************************************************************
  ! *****************************************************************
  
  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

  end module bgmalgp
  