*     *************************************************************************

      subroutine lnsrch( alpha, alphas, alphab,  trial, dir, ds, x,
     +                   osol, rgra, rgnor,  dbd, newp, fcalc, lfc,
     +                   fuval, fx,  difval,  oldfi, nefval, elvar,
     +                   elptr, xst, bound,  ba, nba, su, nsu, nqa,
     +                   ftstem, bep, beptr, kdist, inform, ipdevc,
     +                   nflag, iflag )

*     *************************************************************************

*     Purpose :
*     ---------

*     Once the superbasic search direction DS is computed, this
*     routine computes a step of length ALPHA along this direction
*     such that the trial point is still feasible and the cost
*     function sufficiently reduced.

*     Since various bounds on the variables are present in the
*     calculation, the technique used in this linesearch procedure
*     is similar to that used by Bertsekas.

*     Parameters :
*     ------------

*     alpha   ( dble )

*       input  : meaningless.
*       output : the computed value of the steplength.

*     alphas  ( dble )

*       input  : meaningless.
*       output : The maximum steplength such that the feasibility
*                on the superbasic variables remains feasible.

*     alphab  ( dble )

*       input  : meaningless.
*       output : The maximum steplength such that the feasibility
*                on the basic variables remains feasible.

*     trial   ( int )

*       input  : meaningless.
*       output : the number of linesearch iterations performed
*                before finding a new and better feasible point.

*     dir     ( dble )

*       input  : the superbasic components of the direction.
*       output : the basic components of the projected direction
*                are computed.

*     ds      ( dble )

*       input  : vector used as workspace.
*       output : meaningless.

*     x       ( dble )

*       input  : the current iterate vector.
*       output : a new and better iterate is produced.

*     osol    ( dble )

*       input  : array used as workspace.
*       output : it contains the previous iterate.

*     rgra    ( dble )

*       input  : the reduced gradient vector.
*       output : unmodified.

*     rgnor   ( dble )

*       input  : the euclidean norm of the reduced gradient vector.
*       output : unmodified.

*     dbd     ( dble )

*       input  : array used as workspace.
*       output : vector containing for each superbasic variable
*                the absolute difference between its value and
*                its active bound value.

*     newp    ( log )

*       input  : meaningless.
*       output : .true.  iff a new and better feasible point
*                        has been found.
*                .false. otherwise.

*     fcalc   ( int )

*       input  : meaningless.
*       output : this vector contains the indices of the element
*                functions for which information is required.

*     lfc     ( int )

*       input  : meaningless.
*       output : the number of element functions for which
*                information is required.

*     fuval   ( dble )

*       input  : array used to store the function and derivative
*                values of the element functions.
*       output : the element functions whose indices are in FCALC(i)
*                for i=1 up to i=LFC are evaluated at new point X.

*     fx      ( dble )

*       input  : the objective function value.
*       output : the objective function value evaluated at
*                new point X.

*     difval  ( dble )

*       input  : meaningless.
*       output : the decrease in the objective function value.

*     oldfi   ( dble )

*       input  : array used as workspace.
*       output : it contains the element function values
*                evaluated at the previous iterate.

*     nefval  ( int )

*       input  : the number of element function evaluations.
*       output : this number is increased by one each time
*                an element function needs to be evaluated.

*     elvar   ( int )

*       input  : array containing the indices of the varaibles
*                in the first element, followed by those in the
*                second element, etc.
*       output : unmodified.

*     elptr   ( int )

*       input  : array whose kth value is the position of
*                the first variable of element k, in the
*                list ELVAR.
*       output : unmodified.

*     xst     ( int )

*       input  : vector containing the status of the variables.
*       output : unmodified.

*     bound   ( int )

*       input  : meaningless.
*       output : 0 iff no basic and superbasic variables hit a bound.
*                1 iff at least one basic variable hits a bound.
*                2 iff no basic variable hit a bound and a least one
*                  superbasic variable hits a bound.

*     ba      ( int )

*       input  : vector containing the indices of the basic
*                variables.
*       output : unmodified.

*     nba     ( int )

*       input  : the number of basic variables.
*       output : unmodified.

*     su      ( int )

*       input  : vector containing the indices of the superbasic
*                variables.
*       output : unmodified.

*     nsu     ( int )

*       input  : the number of superbasic variables.
*       output : unmodified.

*     nqa     ( int )

*       input  : the number of quasi-active variables.
*       output : unmodified.

*     ftstem  ( int )

*       input  : vector containing for each superbasic variable,
*                the length from its origine node to the stem
*                node, and the length of its flow augmenting path.
*       output : unmodified.

*     bep     ( int )

*       input  : array containing the flow augmenting paths of
*                the superbasic variables.
*       output : unmodified.

*     beptr   ( int )

*       input  : array whose kth value is the position of the
*                first element of the flow augmenting path of
*                superbasic arc k, in array BEP.
*       output : unmodified.

*     kdist   ( int )

*       input  : vector containing for each superbasic variable
*                the length of its flow augmenting path.
*       output : unmodified.

*     inform  ( int )

*       input  : must be set to zero the first time this routine
*                is called.
*                If it is not equal to zero, the routine is then
*                re-entered with some information that was needed
*                to proceed the calculation.
*       output : If it is not equal to zero, some element function
*                values are needed to proceed the calculation.
*                A return is then made to the main program.
*                Otherwise, no further information is needed.

*     ipdevc  ( int )

*       input  : output device unit number for printing messages.
*       output : unmodified.

*     nflag   ( log )

*       input  : meaningless.
*       output : .true. iff the difference between two successive
*                       function values along the current linesearch
*                       is insignificant compared to the noise on
*                       these function values.
*                .false. otherwise.

*     iflag   ( int )

*       input  : should be equal to zero when the routine is called.
*       output : = 18, if the directional derivative at the beginning
*                      of the linesearch procedure is non-negative.
*                Otherwise, it remains unmodified.

*     Routines used :
*     ---------------

*     sqrt, abs, min, max.
*     dcopy, dsetcd, active, xpdbs, dcopcd, bding.

*     Programming :
*     -------------

*     D. Tuyttens

*     ========================================================================

*     Routine parameters

      integer           ba(*),  nba,  su(*), nsu,  nqa,  ftstem(*),
     +                  bep(*), beptr(*),  nefval,  fcalc(*),  lfc,
     +                  trial, elvar(*), elptr(*), oldfi(*), bound,
     +                  xst(*),  kdist(*),  ipdevc, iflag,  inform
      logical           newp, nflag
      double precision  alpha, alphas, alphab, dir(*), ds(*), x(*),
     +                  osol(*),   rgra(*),  rgnor,   dbd(*),   fx,
     +                  difval,  fuval(*)

*     Internal variables

      integer           i, j, ik, k, kk,  kd, kdis, inoise, ifc
      logical           move,   bfeas
      double precision  sqrem,  alpham,  lambda,  rd,  rdini,  ard,
     +                  bdact,  f0, mindec,  calfa,  dsmax,
     +                  dsnor,  slope

*     Save variables

      save              i, j, ik, k, kk,  kd, kdis, inoise, ifc
      save              move,   bfeas
      save              sqrem,  alpham,  lambda,  rd,  rdini,  ard,
     +                  bdact,  f0, mindec,  calfa,  dsmax,
     +                  dsnor,  slope

      double precision  active
      external          active

*     Common specifications

      integer           arcs, nodes, elem
      double precision  epsmch, huge, tiny, tol
      common / prbdim / arcs, nodes, elem
      common / prbmch / epsmch, huge, tiny, tol

      double precision  zero, one, two, three, half, tenm1, tenm2, tenm4
      common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4
*     Parameter definition.

      integer ttrial
      parameter( ttrial = 10 )

*
*     Reverse communication test.
*
      if( inform.ne.0 ) goto 100
*
*     Some initializations.
*
      sqrem  = sqrt(epsmch)
      alphas = huge
      calfa  = half**ttrial
      rd     = zero
      inoise = 0
      bound  = 0
      newp   = .false.
      move   = .false.
      nflag  = .false.
*
      call dcopy( arcs, x, 1, osol, 1 )
      call dsetcd( nba, ba, dir, zero )
*
*     We compute the maximum steplength ALPHAS,
*     such that the feasibility on the superbasic
*     variables remains feasible.
*
      do 20 ik = 1 , nsu
        k = su(ik)
*
*       The direction is set to minus the reduced
*       gradient value for the quasi active variables.
*
        if( ik.gt.nsu-nqa ) dir(k) = -rgra(k)
*
*       We test if superbasic component K of the
*       direction is too small to be considered.
*
        if( abs(dir(k)).gt.sqrem ) then
          bdact  = active(k,dir)
          dbd(k) = abs( x(k)-bdact )
          alpha  = dbd(k) / abs(dir(k))
*
*         We test if variable K is near its bound.
*
          if( dbd(k).le.tol*(one+abs(bdact))
     +        .or. alpha.le.sqrem ) then
            dir(k) = zero
          else
            move   = .true.
            rd     = rd + rgra(k)*dir(k)
            rdini  = rd
            alphas = min( alphas , alpha )
            if( kd.eq.0 .or. abs(dir(k)).gt.abs(dir(kd)) ) kd = k
          endif
*
*         We computes the basic components of
*         the direction.
*
          kdis = kdist(ik)
          do 30 i = beptr(k) , beptr(k)+kdis-1
            kk = bep(i)
            if( kk.lt.0 ) then
              dir(-kk) = dir(-kk) - dir(k)
            else
              dir(kk)  = dir(kk)  + dir(k)
            endif
 30       continue
        endif
 20   continue
*
*     No move is possible. No
*     new point can be found.
*
      if( .not.move ) return
*
*     We compute the euclidean norm DSNOR
*     of the superbasic direction.
*
      dsnor = zero
      dsmax = abs(dir(kd))
      do 110 ik = 1 , nsu
        k = su(ik)
        if( k.ne.kd ) then
          dsnor = dsnor + (dir(k)/dsmax)**2
        endif
 110  continue
      dsnor = dsmax * sqrt(one+dsnor)
*
*     We test if the superbasic direction is descent ?
*
      slope = rd / (rgnor*dsnor)
      if( slope.ge.zero ) then
*
*       The superbasic direction is not descent.
*
        iflag = 18
        return
      endif
*
*     We compute the maximum steplength ALPHAB,
*     such that the feasibility on the basic
*     variables remains feasible.
*
      alphab = huge
      do 40 ik = 1 , nba
        k = ba(ik)
*
*       We test if basic component K of the
*       direction is too small to be considered.
*
        if( abs(dir(k)).gt.sqrem ) then
          bdact  = active(k,dir)
          alpha  = abs( (x(k)-bdact)/dir(k) )
*
*         We test if variable K is near its bound.
*
          if( abs(x(k)-bdact).le.tol*(one+abs(bdact))
     +        .or. alpha.le.epsmch ) then
*
*           No move is possible. No
*           new point can be found.
*
            alphab = zero
            return
          endif
          alphab = min( alpha , alphab )
        else
          dir(k) = zero
        endif
 40   continue
*
*     Depending on the computed values
*     ALPHAS and ALPHAB, we define the
*     first value of ALPHA.
*
      if( alphab.le.alphas ) then
        alpha = min( one , alphab )
        bfeas = .true.
      else
        alpha = one
        if( alphas.ge.one ) then
          bfeas = .true.
        else
          bfeas = .false.
        endif
      endif
      alpham = min( alphab , alphas )
*
*     The linesearch iteration loop.
*
      if ( ttrial .ge. 1 ) then
        trial = 1
 50     continue
*
*       We build the feasible basic-superbasic direction
*       corresponding to the value of steplength ALPHA.
*
        call xpdbs( alpha,  alpham,  x,  dir, ds, dbd, rgra, rd,
     +              bfeas, ba, nba, su, nsu, ftstem, bep, beptr,
     +              kdist )
*
*       We define the minimum function decrease required.
*
        mindec = tenm2 * rd * alpha
*
*       We select the element function that need
*       to be evaluated at new point X.
*
        lfc = 0
        do 60 i = 1 , elem
          do 70 j = elptr(i) , elptr(i+1)-1
            k = elvar(j)
            if( x(k).ne.osol(k) ) then
              lfc = lfc + 1
              fcalc(lfc) = i
              oldfi(i) = fuval(i)
              nefval   = nefval + 1
              go to 60
            endif
 70       continue
 60     continue
        inform = 29
        return
*
*       We compute the new objective function value.
*
 100    continue
        inform = 0
        f0 = fx
        fx = zero
        do 80 i = 1 , elem
          fx = fx + fuval(i)
 80     continue
*
*       We test if the decrease in the objective
*       function value is sufficient.
*
        difval = fx - f0
        if( difval.gt.mindec ) then
*
*         No sufficient decrease.
*
          fx = f0
          do 90 ifc = 1 , lfc
            i = fcalc(ifc)
            fuval(i) = oldfi(i)
 90       continue
        else
*
*         A new feasible point is found.
*
          newp = .true.
        endif
*
*       We test if the decrease in the objective function value
*       is insignificant compared to the noise on the function
*       value.
*
        if( abs(difval).le. elem*epsmch*(one+abs(f0)) ) then
          inoise = inoise + 1
          if( inoise.ge.2 ) then
            write(ipdevc,1001)
 1001       format('  Function value in the noise !!! ')
            nflag = .true.
            fx = f0
            call dcopcd( nba, ba, osol, x )
            call dcopcd( nsu, su, osol, x )
            return
          endif
        endif
*
*       We check if a variable hits a bound.
*
        if( newp .or. ( inoise.eq.1 .and.
     +      alpha.eq.alphab .and. alpha.le.calfa ) ) then
          call bding( bound, ba, nba, su, nsu, x, xst )
          return
        endif
*
*       No new point has been found, we reset
*       point X to its original values.
*
        call dcopcd( nba, ba, osol, x )
        call dcopcd( nsu, su, osol, x )
*
*       We update the steplength ALPHA.
*
        if( trial.eq.3 .and. alpha.gt.alphas ) then
*
*         We avoid a large number of halving steps.
*
          alpha = alphas
        else
          if( alpha.gt.alphas ) then
*
*           Bissection.
*
            alpha = half * alpha
          else
*
*           Safeguarded quadratic interpolation.
*
            ard    = alpha * rdini
            lambda = - half * ard / ( difval - ard )
            alpha  = min( max( tenm2 , lambda) , half ) * alpha
          endif
        endif
        if( alpha.le.alpham ) bfeas = .true.
*
        trial = trial + 1
        if ( trial .le. ttrial ) go to 50
      endif
*
*     The maximum number of linesearch iterations has been
*     performed without finding a better feasible point.
*
      return
      end