C     *****************************************************************
C     *****************************************************************

      subroutine tnls(nind,x,l,u,m,lambda,rho,equatn,linear,f,g,amax,d,
     +rbdnnz,rbdind,rbdtype,xp,fp,gp,lsinfo,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer inform,lsinfo,m,nind,rbdnnz
      double precision amax,f,fp

C     ARRAY ARGUMENTS
      integer rbdind(nind)
      character rbdtype(nind)
      logical equatn(m),linear(m)
      double precision d(nind),g(*),gp(*),l(nind),lambda(m),rho(m),
     +        u(nind),x(*),xp(*)

C     This subroutine computes a line serach in direction d.
C
C     lsinfo:
C
C     At the first trial:
C
C     4: Interior first trial. Armijo and beta-condition hold.
C     5: x + amax d is at the boundary and f(x + amax d) is smaller than f

C     Extrapolation:
C
C     2: Unbounded objective function?
C     6: Maximum number of extrapolations reached
C     7: Similar consecutive projected points
C     8: Not-well-defined objective function
C     9: Functional value increases

C     Backtracking:
C
C     0: Armijo satisfied
C     1: Small step with functional value similar to the current one
C     2: Unbounded objective function?
C     3: Too small backtracking step. Wrong gradient?

#include "dim.par"
#include "machconst.com"
#include "algconst.par"
#include "counters.com"
#include "outtyp.com"

C     LOCAL SCALARS
      logical boundary,projected,samep
      integer extrap,i
      double precision alpha,atmp,dsupn,fbext,ftmp,gptd,gtd,xsupn

C     LOCAL ARRAYS
      double precision xbext(nmax),xtmp(nmax)

C     EXTERNAL SUBROUTINES
      external calcal

C     ==================================================================
C     ==================================================================

C     Test Armijo condition and beta condition. Decide between accept 
C     the first trial, extrapolate or to do backtracking. 

C     ==================================================================
C     ==================================================================

C     ------------------------------------------------------------------
C     Compute directional derivative, dsupn and xsupn
C     ------------------------------------------------------------------

      gtd = 0.0d0
      dsupn = 0.0d0
      xsupn = 0.0d0
      do i = 1,nind
          gtd = gtd + g(i) * d(i)
          dsupn = max( dsupn, abs( d(i) ) )
          xsupn = max( xsupn, abs( x(i) ) )
      end do

      if ( iprintinn .ge. 6 ) then
          write(*, 100) xsupn,amax,dsupn
          write(10,100) xsupn,amax,dsupn
      end if

C     ------------------------------------------------------------------
C     Compute first trial (projected point)
C     ------------------------------------------------------------------

      alpha = 1.0d0

      boundary = .false.
      do i = 1,nind
          xp(i) = x(i) + d(i)

          if ( xp(i) .lt. l(i) .or. xp(i) .gt. u(i) ) then
              boundary = .true.
              xp(i) = max( l(i), min( xp(i), u(i) ) )
          end if
      end do

      if ( amax .eq. 1.0d0 ) then
          boundary = .true.

          do i = 1,rbdnnz
              if ( rbdtype(i) .eq. 'L' ) then
                  xp(rbdind(i)) = l(rbdind(i))
              else if ( rbdtype(i) .eq. 'U' ) then
                  xp(rbdind(i)) = u(rbdind(i))
              end if
          end do
      end if

      call calcal(nind,xp,m,lambda,rho,equatn,linear,fp,inform)
      if ( inform .lt. 0 ) return

      if ( .not. boundary ) then

          if ( iprintinn .ge. 6 ) then
              write(*, 110) fp,fcnt
              write(10,110) fp,fcnt
          end if

      else
          if ( iprintinn .ge. 6 ) then
              write(*, 120) fp,fcnt
              write(10,120) fp,fcnt
          end if
      end if

C     ------------------------------------------------------------------
C     The first trial is an interior point.
C     ------------------------------------------------------------------

      if ( .not. boundary ) then

          if ( iprintinn .ge. 6 ) then
              write(*, 140)
              write(10,140)
          end if

C         Armijo condition holds.

          if ( fp .le. f + alpha * gamma * gtd ) then

              if ( iprintinn .ge. 6 ) then
                  write(*, 150)
                  write(10,150)
              end if

              call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform)
              if ( inform .lt. 0 ) return

              gptd = 0.0d0
              do i = 1,nind
                  gptd = gptd + gp(i) * d(i)
              end do

C             Beta condition holds. Line search is over.

              if ( gptd .ge. beta * gtd ) then

                  if ( iprintinn .ge. 6 ) then
                      write(*, 160)
                      write(10,160)
                  end if

                  lsinfo = 4

                  if ( iprintinn .ge. 6 ) then
                      write(*, 900)
                      write(10,900)
                  end if

                  return

              end if

C             Beta condtion does not holds. We will extrapolate.

              if ( iprintinn .ge. 6 ) then
                  write(* ,170)
                  write(10,170)
              end if

              go to 1000

          end if

C         Armijo condition does not hold. We will do backtracking.  

          if ( iprintinn .ge. 6 ) then
              write(* ,180)
              write(10,180)
          end if

          go to 2000

      end if

C     ------------------------------------------------------------------
C     The first trial is at the boundary.
C     ------------------------------------------------------------------

      if ( iprintinn .ge. 6 ) then
          write(*, 190) 
          write(10,190) 
      end if

C     Function value is smaller than at the current point. We will 
C     extrapolate.

      if ( fp .lt. f ) then

          if ( iprintinn .ge. 6 ) then
              write(*, 200) 
              write(10,200) 
          end if

          go to 1000

      end if

C     Discard the projected point and consider x + amax d

      if ( iprintinn .ge. 6 ) then
          write(*, 210) 
          write(10,210) 
      end if

      alpha = amax

      do i = 1,nind
          xp(i) = x(i) + alpha * d(i)
      end do

      do i = 1,rbdnnz
          if ( rbdtype(i) .eq. 'L' ) then
              xp(rbdind(i)) = l(rbdind(i))
          else if ( rbdtype(i) .eq. 'U' ) then
              xp(rbdind(i)) = u(rbdind(i))
          end if
      end do

      call calcal(nind,xp,m,lambda,rho,equatn,linear,fp,inform)
      if ( inform .lt. 0 ) return

      if ( iprintinn .ge. 6 ) then
          write(*, 130) alpha,fp,fcnt
          write(10,130) alpha,fp,fcnt
      end if

C     Function value is smaller than or equal to (or even just a little 
C     bit greater than) at the current point. Line search is over.

      if ( fp .lt. f + macheps23 * abs(f) ) then

          if ( iprintinn .ge. 6 ) then
              write(*, 220) 
              write(10,220) 
          end if

          call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform)
          if ( inform .lt. 0 ) return

          lsinfo = 5

          if ( iprintinn .ge. 6 ) then
              write(*, 900)
              write(10,900)
          end if

          return

      end if

C     Function value is greater than at the current point. We will 
C     do backtracking.

      if ( iprintinn .ge. 6 ) then
          write(*, 230) 
          write(10,230) 
      end if

      go to 2000


C     ==================================================================
C     ==================================================================

C     Extrapolation

C     ==================================================================
C     ==================================================================

 1000 continue

      extrap = 0

C     Save f and x before extrapolation to return in case of a
C     not-well-defined objective function at an extrapolated point

      fbext = fp

      do i = 1,nind
          xbext(i) = xp(i)
      end do

 1010 continue

C     Test f going to -inf

      if ( fp .le. fmin ) then

C         Finish the extrapolation with the current point

          call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform)
          if ( inform .lt. 0 ) return

          lsinfo = 2

          if ( iprintinn .ge. 6 ) then
              write(*, 910)
              write(10,910)
          end if

          return

      end if

C     Test if the maximum number of extrapolations was exceeded

      if ( extrap .ge. maxextrap ) then

C         Finish the extrapolation with the current point

          call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform)
          if ( inform .lt. 0 ) return

          lsinfo = 6

          if ( iprintinn .ge. 6 ) then
              write(*, 930)
              write(10,930)
          end if

          return

      end if

C     Chose new step 

      extrap = extrap + 1

      if ( alpha .lt. amax .and. etaext * alpha .gt. amax ) then
          atmp = amax
      else
          atmp = etaext * alpha
      end if

C     Compute new trial point

      do i = 1,nind
          xtmp(i) = x(i) + atmp * d(i)
      end do

      if ( atmp .eq. amax ) then
          do i = 1,rbdnnz
              if ( rbdtype(i) .eq. 'L' ) then
                  xtmp(rbdind(i)) = l(rbdind(i))
              else if ( rbdtype(i) .eq. 'U' ) then
                  xtmp(rbdind(i)) = u(rbdind(i))
              end if
          end do
      end if

C     Project

      if ( atmp .gt. amax ) then

          projected = .false.
          do i = 1,nind
              if ( xtmp(i) .lt. l(i) .or. xtmp(i) .gt. u(i) ) then
                  projected = .true.
                  xtmp(i) = max( l(i), min( xtmp(i), u(i) ) )
              end if
          end do

      end if

C     Test if this is not the same point as the previous one. This test 
C     is performed only when xtmp is in fact a projected point.

      if ( projected ) then

          samep = .true.
          do i = 1,nind
              if ( abs( xtmp(i) - xp(i) ) .gt. 
     +             macheps * max( 1.0d0, abs( xp(i) ) ) ) then
                  samep = .false.
              end if
          end do

          if ( samep ) then

C             Finish the extrapolation with the current point

              call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform)
              if ( inform .lt. 0 ) return

              lsinfo = 7

              if ( iprintinn .ge. 6 ) then
                  write(*, 940)
                  write(10,940)
              end if

              return

          end if

      end if

C     Evaluate function

      call calcal(nind,xtmp,m,lambda,rho,equatn,linear,ftmp,inform)

C     If the objective function is not well defined in an extrapolated 
C     point, we discard all the extrapolated points and return to a 
C     safe region (to the last point before the extrapolation)

      if ( inform .lt. 0 ) then

          fp = fbext

          do i = 1,nind
              xp(i) = xbext(i)
          end do

          call csetp(nind,xp)

          call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform)
          if ( inform .lt. 0 ) return

          lsinfo = 8

          if ( iprintinn .ge. 6 ) then
              write(*, 950)
              write(10,950)
          end if

          return

      end if

C     Print information of this iteration

      if ( iprintinn .ge. 6 ) then
          write(*, 130) atmp,ftmp,fcnt
          write(10,130) atmp,ftmp,fcnt
      end if

C     If the functional value decreases then set the current point and 
C     continue the extrapolation

      if ( ftmp .lt. fp ) then

          alpha = atmp

          fp = ftmp

          do i = 1,nind
              xp(i) = xtmp(i)
          end do

          go to 1010

      end if

C     If the functional value does not decrease then discard the last 
C     trial and finish the extrapolation with the previous point

      call csetp(nind,xp)

      call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform)
      if ( inform .lt. 0 ) return

      lsinfo = 9

      if ( iprintinn .ge. 6 ) then
          write(*, 960)
          write(10,960)
      end if

      return

C     ==================================================================
C     End of extrapolation
C     ==================================================================

C     ==================================================================
C     ==================================================================

C     Backtracking

C     ==================================================================
C     ==================================================================

 2000 continue

      call backtracking(nind,x,m,lambda,rho,equatn,linear,f,xsupn,d,gtd,
     +dsupn,alpha,fp,xp,calcal,lsinfo,inform)
      if ( inform .lt. 0 ) return

      call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform)
      if ( inform .lt. 0 ) return

C     ==================================================================
C     End of backtracking
C     ==================================================================

C     NON-EXECUTABLE STATEMENTS

 100  format(/,5X,'TN Line search (xsupn = ',1P,D7.1,', amax = ',
     +             1P,D7.1,', dsupn = ',1P,D7.1,')')
 110  format(  5X,'Unitary step    F = ',1P,D24.16,' FE = ',I7)
 120  format(  5X,'Projected point F = ',1P,D24.16,' FE = ',I7)
 130  format(  5X,'Alpha = ',1P,D7.1,' F = ',1P,D24.16,' FE = ',I7)
 140  format(  5X,'The first trial is an interior point.')
 150  format(  5X,'Armijo condition holds.')
 160  format(  5X,'Beta condition also holds. Line search is over.')
 170  format(  5X,'Beta condition does not hold. We will extrapolate.')
 180  format(  5X,'Armijo condition does not hold. We will backtrack.')
 190  format(  5X,'The first trial is at the boundary.')
 200  format(  5X,'Function value at the boundary is smaller than at ',
     +            'the current point.',/,5X,'We will extrapolate.')
 210  format(  5X,'Discarding projected point. We will now consider x ',
     +            '+ amax d.')
 220  format(  5X,'Function value at the boundary is smaller than or ',
     +            'equal to than at the',/,5X,'current point. Line ',
     +            'search is over.')
 230  format(  5X,'Function value at the boundary is greater than at ',
     +            'the current point.',/,5X,'We will backtrack.')
 900  format(  5X,'Flag of TN Line search: First trial accepted.')
 910  format(  5X,'Flag of TN Line search: Unbounded objective ',
     +            'function?')
 930  format(  5X,'Flag of TN Line search: Maximum of consecutive ',
     +            'extrapolations reached.')
 940  format(  5X,'Flag of TN Line search: Very similar consecutive ',
     +            'projected points.')
 950  format(  5X,'Flag of TN Line search: Not-well-defined objective ',
     +            'function in an extrapolated point.')
 960  format(  5X,'Flag of TN Line search: Functional value increased ',
     +            'when extrapolating.')

      end

