      subroutine ivmispg(n,m,x,y,arow,acol,aval,anze,b,epsdsn,epsden,
     +                   epsnlsn,maxoutit,maxinnit,maxfc,p,extrap,hess,
     +                   output,kprint,eta,theta,tau,tol,lmin,lmax,
     +                   gamma,sterel,steabs,epsmac,infty,tau1,tau2,f,
     +                   dsupn,deucn2,amax,outit,innit,fcnt,gcnt,hcnt,
     +                   flag)

C     Subroutine IVM implements the Inexact Variable Metric method
C     to find the local minimizers of a given function with convex 
C     constraints (linear inequality constraints in this implementation), 
C     described in:
C
C     R. Andreani, E. G. Birgin, J. M. Martínez and J. Yuan, "Spectral 
C     projected gradient and variable metric methods for optimization 
C     with linear inequalities", IMA Journal of Numerical Analysis 25, 
C     pp. 221-252, 2005.
C
C     The user must supply the external subroutines evalf, evalg, 
C     and evalH to evaluate the objective function, its gradient and
C     its hessian, respectively. Depending on parameter hess, evalH
C     subroutine may be an empty subroutine.
C
C     This version 22 SET 2003 by E.G.Birgin, J.M.Martinez and J.Yuan.
C     Final revision ... by E.G.Birgin, J.M.Martinez and J. Yuan.
C
C     ON ENTRY:
C
C     n     integer,
C           number of variables,
C
C     m     integer,
C           number of linear inequality constraints,
C
C     x     double precision x(n),
C           initial estimation of the solution,
C
C     y     double precision y(m),
C           initial estimation of the Lagrange multipliers,
C
C     lda   integer
C           leading dimension of matrix A (see below),
C
C     A     double precision A(m,n),
C           linear inequalities coefficients matrix,
C
C     b     double precision b(m),
C           linear inequalities right hand side,
C
C     epsdsn double precision,
C           tolerance for the stopping criterion || amax d ||_sup < epsdsn,
C           where d is the IVM search direction and amax is
C           the maximum step which preserves feasibility,
C
C     epsden double precision,
C           tolerance for the stopping criterion || amax d ||_2 < epsden,
C           where d is the IVM search direction and amax is
C           the maximum step which preserves feasibility,
C
C     epsnlsn double precision,
C           tolerance for the stopping criterion of the inner
C           solver based on the sup-norm of the projected gradient
C           of the Lagrangian function,
C     
C     maxoutit integer,
C           maximum number of outer iterations,
C
C     maxinnit integer,
C           maximum number of inner iterations,
C
C     maxfc integer,
C           maximum number of function evaluations,
C
C     p     integer,
C           number of previous function values to be considered 
C           in the nonmonotone line search,
C
C           RECOMMENDED = 10
C
C     extrap logical,
C           if extrap = TRUE then, when solving the subproblem (using
C           an active set strategy), extrapolations are tried to increase
C           the dimension of the active set (it may be recomended for high
C           dimensional subproblems which comes from original problems with
C           many constraints); else extrapolations are not used at all.
C
C           RECOMMENDED = .true.
C
C     hess  logical,
C           if hess = TRUE then the true Hessian is used in the
C           quadratic subproblems (the Identity is used if the function
C           is linear); else the spectral approximation is used. 
C
C           RECOMMENDED = try both!
C
C     output integer,
C           < 0 means no print at all;
C           = 0 just initial and final information;
C           = 1 information at each outer iteration;
C           = 2 initial and final information for the inner
C               jobs: inner subproblems and line searches;
C           = 3 information of the inner iterations and 
C               the line search iterations;
C           = 4 detailed information of the inner iterations.
C
C     kprint integer,
C           maximum number of components to be printed when an
C           array must be showed,
C
C           RECOMMENDED = 5
C
C     eta   double precision,
C           eta in (0,1] controls the inexactness for solving the subproblem. 
C           eta = 1 means "solve the subproblem exactly".
C
C           RECOMMENDED: 0.9
C
C     theta double precision,
C           theta in (eta,1) is used to warranty that the iterates are
C           interior enough.
C
C           RECOMMENDED: 0.95
C
C     tau   double precision,
C           tau > 0 is used to compute the search direction of a 
C           subproblem of type q(y) = 0.5 y^t B y + b^t y as 
C           (B + tau I) dk = - \nabla q(yk). See formulae (44) and
C           the end of Section 7 for details.
C
C           RECOMMENDED = 10^-6 
C
C     tol   double precision,
C           tol in (0,1) is used, i the algorithm to solve the subproblems
C           (which uses an active set strategy) to decide if the current 
C           face must be abandoned or not. Near 0 means to abandon the
C           face as fast as possible. Near 1 means try to find the solution
C           inside the current face.
C
C           RECOMMENDED: 0.9  
C
C     lmin  double precision,
C           lmin > 0, lower bound for the safeguarded spectral 
C           steplenght,
C
C           RECOMMENDED: 10^-3
C
C     lmax  double precision,
C           lmax >= lmin, upper bound for the safeguarded spectral 
C           steplenght,
C
C           RECOMMENDED: 10^3
C
C     gamma double precision,
C           gamma > 0, constant for the Armijo-like stooping criterion 
C           of the line search,
C
C           RECOMMENDED: 10^-4
C
C     sterel double precision,
C           epsrel > 0, a relative small number used for finite difference
C           approximation of derivatives,
C
C           RECOMMENDED: 10^-7
C
C     steabs double precision,
C           epsabs > 0, an absolute small number used for finite difference
C           approximation of derivatives,
C
C           RECOMMENDED: 10^-10
C
C     epsmac double precision,
C           epsmac > 0, machine epsilon, i.e., the smallest positive value 
C           such that 1 + epsmac > 1.
C
C           RECOMMENDED: In IEEE double precision floating point
C           arithmetic, the machine epsilon is something near 10^-16.
C           But the recomendation is to fix its value running the   
C           following subroutine:         
C
C           subroutine mcheps(epsmac)
C           double precision epsmac
C           epsmac = 1.0d0
C      20   if ( ( 1.0d0 + epsmac ) .ne. 1.0d0 ) then
C               epsmac = epsmac / 2.0d0
C               go to 20
C           end if
C           epsmac = 2.0d0 * epsmac
C           end
C
C     infty double precision,
C           a big number,
C
C           RECOMMENDED = 10^20
C
C     The next two parameters are parameters for subroutine modchl
C     which computes the "Eskow-Schnabel" modified Cholesky factorization
C     of a given symetric matrix, taken from:
C
C     ALGORITHM 695, COLLECTED ALGORITHMS FROM ACM.
C     PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C     VOL. 17, NO. 3, SEPTEMBER, 1991, PP. 306-312.
C
C     tau1  double precision,
C           tau1 > 0, tolerance used for determining when to switch phase 2
C
C           RECOMMENDED: 10^-6 (approx. epsmac ^ 1/3, where epsmac
C           is the machine epsilon (approx. 10^-16).
C
C     tau2  double precision,
C           tolerance used for determining the maximum condition 
C           number of the final 2X2 submatrix.
C
C           RECOMMENDED: 10^-6 (approx. epsmac ^ 1/3, where epsmac
C           tau2 > 0, is the machine epsilon (approx. 10^-16).
C
C     ON RETURN:
C
C     x     double precision x(n),
C           final estimation of to the local minimizer,
C
C     y     double precision y(m),
C           final estimation of the Lagrange multipliers,
C
C     f     double precision,
C           function value at the approximation to the local
C           minimizer,
C
C     dsupn double precision,
C           || d ||_sup at the final iteration,
C           where d is the IVM search direction,
C
C     deucn2 double precision,
C           || d ||_2^2 at the final iteration,
C           where d is the IVM search direction,
C
C     amax double precision,
C           maximum step to preserve feasibility in the last
C           computed direction,
C
C     outit  integer,
C           number of outer iterations,
C
C     innit  integer,
C           number of inner iterations,
C
C     fcnt  integer,
C           number of function evaluations,
C
C     gcnt  integer,
C           number of gradient evaluations,
C
C     hcnt  integer,
C           number of Hessian evaluations,
C
C     flag  integer,
C           termination parameter:
C           0= convergence with sup-norm of the IVM search direction,
C           1= convergence with euclidian-norm of the IVM search direction,

C           6= too many outer iterations,
C           7= too many inner iterations,
C           8= too many functional evaluations,

C          -1= error in evalf subroutine,
C          -2= error in evalg subroutine,
C          -3= error in evalH subroutine.

C     PARAMETERS
      integer mmax
      parameter (mmax    =   100 000)
      integer nmax
      parameter (nmax    =     4 000)
      integer pmax
      parameter (pmax    =       100)

C     SCALAR ARGUMENTS
      logical hess,extrap
      integer anze,fcnt,flag,gcnt,hcnt,innit,kprint,m,maxfc,maxinnit,
     +        maxoutit,n,outit,output,p
      double precision amax,dsupn,deucn2,epsdsn,epsden,epsnlsn,f,eta,
     +        theta,tau,tol,lmin,lmax,gamma,tau1,tau2,sterel,steabs,
     +        epsmac,infty

C     ARRAY ARGUMENTS
      integer arow(anze),acol(anze)
      double precision aval(anze),b(m),x(n),y(m)

C     COMMON SCALARS
      double precision lambda

C     COMMON ARRAYS
      double precision g(nmax),bbar(mmax)

C     LOCAL SCALARS
      integer i,inform,nprint,mprint
      double precision fbest,fnew,gtd,redc,L,Q,nLpsupn,alpha,atemp,fmax,
     +        sts,sty,gtg,tsmall,nLpi

C     LOCAL ARRAYS
      double precision xnew(nmax),xtmp(nmax),xbest(nmax),gnew(nmax),
     +        gtmp(nmax),d(nmax),ybest(mmax),nL(mmax),lastfv(0:pmax-1),
     +        w1(nmax),w3(nmax),w4(nmax),w5(mmax)

C     GENCAN VARIABLES
      logical Gnearlyq
      integer Gmaxitnfp,Gmaxitngp,Gmaxit,Gmaxfc,Gcgmaxit,Gcgscre,
     +        Gmaxitnqmp,Ggtype,Ghtvtype,Gtrtype,Giprint,Giter,Gfcnt,
     +        Ggcnt,Gcgcnt,Gspgiter,Gspgfcnt,Gtniter,Gtnfcnt,Gtnstpcnt,
     +        Gtnintcnt,Gtnexgcnt,Gtnexbcnt,Gtnintfe,Gtnexgfe,Gtnexbfe,
     +        Gflag,Gw2(mmax),Gmininterp,Gncomp
      double precision Gepsgpen,Gepsgpsn,Gepsnfp,Gfmin,Gdelta0,Gcggpnf,
     +        Gcgepsi,Gcgepsf,Gepsnqmp,Gf,Gg(mmax),Ggpeucn2,Ggpsupn,
     +        Gw11(mmax),Gw12(mmax),Gw13(mmax),Gw14(5*mmax),Geta,
     +        Gdelmin,Glammax,Glammin,Gtheta,Ggamma,Gbeta,Gsigma1,
     +        Gsigma2,Gsterel,Gsteabs,Gepsrel,Gepsabs,Ginfty,Gnint,
     +        Gnext,Gl(mmax),Gu(mmax)

C     COMMON BLOCKS
      common /ldata2/ g,bbar,lambda

C     EXTERNAL SUBROUTINES
      external evallnlispg,gencan

C     INTRINSIC FUNCTIONS
      intrinsic abs,max,min,mod,sqrt

C     ===============================================================
C     INEXACT VARIABLE METRIC METHOD
C     ===============================================================

C     ===============================================================
C     INITALIZATION
C     ===============================================================

C     PRINT PRESENTATION INFORMATION

      if ( output .ge. 0 ) then
          write(*, fmt=9010)
          write(10,fmt=9010)
      end if

C     CONSTANTS

      nprint = min(kprint,n)
      mprint = min(kprint,m)

C     COUNTERS

      outit = 0
      innit = 0
      fcnt  = 0
      gcnt  = 0
      hcnt  = 0

C     INITIALIZE LAST FUNCTIONAL VALUES VECTOR

      do i = 0,p - 1
          lastfv(i) = - infty
      end do

C     EVALUATE THE FUNCTION

      call evalf(n,x,f,inform)
      fcnt = fcnt + 1

      if ( inform .ne. 0 ) then
        ! ERROR IN EVALF SUBROUTINE, STOP
          flag = -1
          go to 900
      end if

C     EVALUATE THE GRADIENT

      call evalg(n,x,g,inform)
      gcnt = gcnt + 1

      if ( inform .ne. 0 ) then
        ! ERROR IN EVALG SUBROUTINE, STOP
          flag = -2
          go to 900
      end if

C     STORE FUNCTION VALUE FOR THE NONMONOTONE LINE SEARCH

      lastfv(0) = f

C     INITIALIZE BEST SOLUTION, BEST FUNCTIONAL VALUE AND ITS
C     ASSOCIATED LAGRANGE MULTIPLIERS

      fbest = f

      do i = 1,n
          xbest(i) = x(i)
      end do

      do i = 1,m
          ybest(i) = y(i)
      end do

C     DEFINE INITIAL SPECTRAL STEPLENGHT

      tsmall = 0.0d0
      do i = 1,n
          tsmall = max( tsmall, abs( x(i) ) )
      end do

      tsmall = max( sterel * tsmall, steabs )          

      do i = 1,n
          xtmp(i) = x(i) - tsmall * g(i)
      end do

      call evalg(n,xtmp,gtmp,inform)
      gcnt = gcnt + 1

      gtg = 0.0d0
      do i = 1,n
          gtg = gtg + g(i) ** 2
      end do

      sts = tsmall ** 2 * gtg

      sty = 0.0d0
      do i = 1,n
          sty = sty + g(i) * ( gtmp(i) - g(i) )
      end do
      sty = - tsmall * sty

      if ( sty .le. 0.0d0 ) then
          lambda = lmax

      else
          lambda = min( lmax, max( lmin, sts / sty ) )
      end if

C     ===============================================================
C     MAIN OUTER LOOP
C     ===============================================================

 100  continue

      outit = outit + 1

C     ===============================================================
C     STEP 1: COMPUTE THE SEARCH DIRECTION
C     ===============================================================

C     ---------------------------------------------------------------
C     Inner subproblem
C     ---------------------------------------------------------------

C     To compute the search direction we want to 
C
C     Minimize (approx) Qk(d) = 0.5 d^t Hk d + gk^t d 
C     subject to A(xk + d) <= b.
C
C     ---> In the case of the ISPG, Hk = (1 / lambda) I

C     Equivalently (via duality), we want to 
C
C     Minimize (approx) -Lk(y) = 0.5 (gk + A^t y) Hk^-1 (gk + A^t y) + 
C                                y^t (b - A xk)
C     subject to y >= 0.
C
C     ---> In the case of the ISPG,
C
C     Minimize (approx) -Lk(y) = 0.5 lambda || (gk + A^t y) ||_2^2 + 
C                                y^t (b - A xk)
C     subject to y >= 0.

C     The primal solution is then obtained by: 
C
C     d = - Hk^-1 (gk + A^t y)
C
C     ---> In the case of the ISPG,
C
C     d = - lambda (gk + A^t y)

C     ---------------------------------------------------------------
C     Initialization
C     ---------------------------------------------------------------

      if ( output .ge. 2 ) then
          write(*, fmt=9200)
          write(*, fmt=9203) lambda

          write(10,fmt=9200)
          write(10,fmt=9203) lambda
      end if

C     Compute bbar = b - A xk

      do i = 1,m
          bbar(i) = b(i)
      end do

      do i = 1,anze
          bbar(arow(i)) = bbar(arow(i)) - aval(i) * x(acol(i))
      end do

c     Just to avoid bbar(i) = -10^-16 due to rounding errors
      do i = 1,m
          bbar(i) = max( 0.0d0, bbar(i) )
      end do

C     ---------------------------------------------------------------
C     Inner main loop
C     ---------------------------------------------------------------

 200  continue

C     ---------------------------------------------------------------
C     Given y, this subroutine computes:
C
C     (a) The corresponding primal point: d = - Hk^-1 (gk + A^t y);
C
C     (b) The maximum step to mantain feasibility of xk + d in the
C         primal feasible set A( xk + alpha d ) <= b;
C
C     (c) The reduction constant to warranty that xk + redc * d is
C         an (enough) interior point;
C
C     (d) Lk(d,y) and Qk(redc d);
C
C     (e) The gradient of -L as a function of y.
C     ---------------------------------------------------------------

      call evallnlispg(m,n,y,g,arow,acol,aval,anze,bbar,lambda,d,amax,
     +                 redc,L,Q,nL,w1,w5,w3,w4,theta,infty)

C     Compute chopped antigradient and projected antigradient norms; 
C     and store in nind the number of free variables and in array nind 
C     its identifiers.

      nLpsupn = 0.0d0
      do i = 1,m
          if ( y(i) .gt. 0.0d0 ) then
              nLpi   = - nL(i)
          else
              nLpi   = max( 0.0d0, - nL(i) )
          end if
          nLpsupn = max( nLpsupn, abs( nLpi   ) )
      end do

C     Print inner iteration information

      if ( output .ge. 3 ) then
          if ( output .ge. 4 ) then
              write(*, fmt=9210) innit,mprint,m,(y(i),i=1,mprint)
              write(*, fmt=9240) -L,nLpsupn
              write(*, fmt=9250) nprint,n,(d(i),i=1,nprint)
              write(*, fmt=9260) amax,redc,Q

              write(10,fmt=9210) innit,mprint,m,(y(i),i=1,mprint)
              write(10,fmt=9240) -L,nLpsupn
              write(10,fmt=9250) nprint,n,(d(i),i=1,nprint)
              write(10,fmt=9260) amax,redc,Q
          else
              write(*, fmt=9205) innit,-L,nLpsupn
              write(10,fmt=9205) innit,-L,nLpsupn
          end if
      end if

C     ---------------------------------------------------------------
C     Inner stooping criteria
C     ---------------------------------------------------------------

C     Testing precognized inner stooping criterion: 
C     Qk(redc d) <= eta L(d,y)

      if ( Q .le. eta * L ) then
        ! Precongnized stooping criterion satisfied, 
        ! stop the inner algorithm
          if ( output .ge. 2 ) then
              write(*, fmt=9270)
              write(10,fmt=9270)
          end if

          go to 300

      end if

C     Testing projected-gradient norm inner stooping criterion

      if ( nLpsupn .le. epsnlsn ) then
        ! Stationary point of the subproblem, 
        ! stop the inner algorithm
          if ( output .ge. 2 ) then
              write(*, fmt=9280)
              write(10,fmt=9280)
          end if

          go to 300

      end if

C     Testing number of inner iterations

      if ( innit .ge. maxinnit ) then
        ! Maximum number of inner iterations exceeded, 
        ! stop the algorithm
          flag = 7
          go to 900
      end if

C     ---------------------------------------------------------------
C     Doing an inner iteration
C     ---------------------------------------------------------------

c     innit = innit + 1

C     GENCAN ITERATION, STARTING FROM y TO MINIMZE
C
C     -L(y) = 0.5 lambda || (gk + A^t y) ||_2^2 + bbar^t y
C
C     SUBJECT TO y >= 0.

      do i = 1,m
          Gl(i) = 0.0d0
          Gu(i) = infty
      end do

      Gepsgpen   =  0.0d0
      Gepsgpsn   =  epsnlsn
      Gepsnfp    =  0.0d0
      Gmaxitnfp  =  10000
      Gmaxitngp  =  10000
      Gfmin      = -1.0d+99
      Gmaxit     =  100
      Gmaxfc     =  10000
      Gdelta0    = -1.0d0
      Gcgmaxit   = -1
      Gcgscre    =  2
      Gcggpnf    =  Gepsgpsn
      Gcgepsi    =  1.0d-1
      Gcgepsf    =  1.0d-5  
      Gepsnqmp   =  1.0d-4
      Gmaxitnqmp =  5
      Gnearlyq   = .true.
      Ggtype     =  0
      Ghtvtype   =  0
      Gtrtype    =  0
      Giprint    = -1
      Geta       =  0.9d0
      Gdelmin    =  0.1d0
      Glammax    =  1.0d+03
      Glammin    =  1.0d-03
      Gtheta     =  1.0d-06
      Ggamma     =  1.0d-04
      Gbeta      =  0.5d0
      Gsigma1    =  0.1d0
      Gsigma2    =  0.9d0
      Gnint      =  2
      Gnext      =  2
      Gmininterp =  4
      Gncomp     =  5
      Gsterel    =  1.0d-07
      Gsteabs    =  1.0d-10
      Gepsrel    =  1.0d-10
      Gepsabs    =  1.0d-20
      Ginfty     =  1.0d+20

      call gencan(m,y,Gl,Gu,Gepsgpen,Gepsgpsn,Gepsnfp,Gmaxitnfp,
     +Gmaxitngp,Gfmin,Gmaxit,Gmaxfc,Gdelta0,Gcgmaxit,Gcgscre,Gcggpnf,
     +Gcgepsi,Gcgepsf,Gepsnqmp,Gmaxitnqmp,Gnearlyq,Ggtype,Ghtvtype,
     +Gtrtype,Giprint,Gf,Gg,Ggpeucn2,Ggpsupn,Giter,Gfcnt,Ggcnt,Gcgcnt,
     +Gspgiter,Gspgfcnt,Gtniter,Gtnfcnt,Gtnstpcnt,Gtnintcnt,Gtnexgcnt,
     +Gtnexbcnt,Gtnintfe,Gtnexgfe,Gtnexbfe,Gflag,Gw11,Gw12,
     +Gw13,Gw2(1),Gw14,Geta,Gdelmin,Glammax,Glammin,Gtheta,
     +Ggamma,Gbeta,Gsigma1,Gsigma2,Gnint,Gnext,Gmininterp,Gncomp,
     +Gsterel,Gsteabs,Gepsrel,Gepsabs,Ginfty)

      do i = 1,m
          if ( y(i) .le. epsmac ) then
              y(i) = 0.0d0
          end if
      end do

      innit = innit + Giter

C     ---------------------------------------------------------------
C     End of the inner iteration
C     ---------------------------------------------------------------

      go to 200

C     ---------------------------------------------------------------
C     End of the inner main loop
C     ---------------------------------------------------------------

C     ===============================================================
C     PRINT OUTER ITERATION INFORMATION AND TEST STOPPING CRITERIA
C     ===============================================================

 300  continue

C     ---------------------------------------------------------------
C     Compute d norms
C     ---------------------------------------------------------------

      dsupn  = 0.0d0
      deucn2 = 0.0d0
      do i = 1,n
          dsupn  = max( dsupn, abs( d(i) ) )
          deucn2 = deucn2 + d(i) ** 2
      end do

C     ---------------------------------------------------------------
C     Print outer iteration information
C     ---------------------------------------------------------------

      if ( output .ge. 1 ) then
          write(*, fmt=9020) outit,f,dsupn,amax
          write(10,fmt=9020) outit,f,dsupn,amax
      end if

C     ---------------------------------------------------------------
C     Stooping criterion
C     ---------------------------------------------------------------

      if ( abs( redc ) * dsupn .le. epsdsn ) then
        ! IVM SEARCH DIRECTION SUP-NORM STOPPING CRITERION SATISFIED, STOP
          flag = 0
          go to 900
      end if

      if ( redc ** 2 * deucn2 .le. epsden ** 2 ) then
        ! IVM SEARCH DIRECTION 2-NORM STOPPING CRITERION SATISFIED, STOP
          flag = 1
          go to 900
      end if

      if ( outit .gt. maxoutit ) then
        ! MAXIMUM NUMBER OF OUTER ITERATIONS EXCEEDED, STOP
          flag = 6
          go to 900
      end if

C     ===============================================================
C     STEP 2: NONMONOTONE LINE SEARCH
C     ===============================================================

C     ---------------------------------------------------------------
C     Line search initialization
C     ---------------------------------------------------------------

      if ( output .ge. 2 ) then
          write(*, fmt=9500)
          write(10,fmt=9500)
      end if

C     Compute g^t d

      gtd = 0.0d0
      do i = 1,n
          gtd    = gtd + g(i) * d(i)
      end do

C     Compute the maximum functional value of the last p iterations

      fmax = lastfv(0)
      do i = 1,p - 1
          fmax = max(fmax,lastfv(i))
      end do

C     Compute first trial

      alpha = redc

      do i = 1,n
          xnew(i) = x(i) + alpha * d(i)
      end do

      call evalf(n,xnew,fnew,inform)
      fcnt = fcnt + 1

      if (inform .ne. 0) then
        ! Error in evalf subroutine, stop
          flag = -1
          go to 900
      end if

C     ---------------------------------------------------------------
C     Line search main loop
C     ---------------------------------------------------------------

 400  continue

C     Print iteration information

      if ( output .ge. 3 ) then
          write(*, fmt=9510) alpha,fnew
          write(10,fmt=9510) alpha,fnew
      end if

C     ---------------------------------------------------------------
C     Line search stooping criteria
C     ---------------------------------------------------------------

      if ( fnew .le. fmax + gamma * alpha * gtd ) then
        ! Nonmonotone Armijo-like stopping criterion satisfied, 
        ! stop the line search

          if ( output .ge. 2 ) then
              write(*, fmt=9520)
              write(10,fmt=9520)
          end if

          go to 500
      end if

      if ( fcnt .ge. maxfc ) then
        ! Maximum number of functional evaluations exceeded,
        ! stop the algorithm
          flag = 8
          go to 900
      end if 

C     ---------------------------------------------------------------
C     Line search iteration
C     ---------------------------------------------------------------

C     Safeguarded quadratic interpolation

      if ( alpha .le. 0.1d0 ) then
          alpha = 0.5d0 * alpha

      else
          atemp = ( - gtd * alpha ** 2 ) / 
     +            ( 2.0d0 * ( fnew - f - alpha * gtd ) )
          if ( atemp .lt. 0.1d0 .or. atemp .gt. 0.9d0 * alpha ) then
              atemp = 0.5d0 * alpha
          end if
          alpha = atemp
      end if

C     Compute trial point

      do i = 1,n
          xnew(i) = x(i) + alpha * d(i)
      end do

C     Evaluate the function

      call evalf(n,xnew,fnew,inform)
      fcnt = fcnt + 1

      if (inform .ne. 0) then
        ! Error in evalf subroutine, stop
          flag = -1
          go to 900
      end if

C     ---------------------------------------------------------------
C     End of the line search iteration
C     ---------------------------------------------------------------

      go to 400

C     ---------------------------------------------------------------
C     End of the line search main loop
C     ---------------------------------------------------------------

C     ===============================================================
C     PREPARING FOR THE NEXT OUTER ITERATION
C     ===============================================================

 500  continue

C     EVALUATE THE GRADIENT AT THE NEW ITERATE

      call evalg(n,xnew,gnew,inform)
      gcnt = gcnt + 1

      if ( inform .ne. 0 ) then
        ! ERROR IN EVALG SUBROUTINE, STOP
          flag = -2
          go to 900
      end if

C     COMPUTE THE SPECTRAL STEPLENGHT

      sts = 0.0d0
      sty = 0.0d0
      do i = 1,n
          sts = sts + ( xnew(i) - x(i) ) ** 2
          sty = sty + ( xnew(i) - x(i) ) * ( gnew(i) - g(i) )
      end do

      if ( sty .le. 0.0d0 ) then
          lambda = lmax

      else
          lambda = min( lmax, max( lmin, sts / sty ) )
      end if

C     SET THE NEW CURRENT POINT AND ITS GRADIENT

      do i = 1,n
          x(i) = xnew(i)
          g(i) = gnew(i)
      end do

C     SET NEW FUNCTION VALUE AND SAVE IT FOR THE NONMONOTONE
C     LINE SEARCH

      f = fnew
      lastfv(mod(outit,p)) = f

C     COMPARE THE NEW FUNCTION VALUE AGAINST THE BEST FUNCTION 
C     VALUE  AND, IF SMALLER, UPDATE THE BEST FUNCTION 
C     VALUE, THE CORRESPONDING BEST POINT AND ITS ASSOCIATED
C     LAGRANGE MULTIPLIERS

      if ( f .lt. fbest ) then
          fbest = f
          do i = 1,n
              xbest(i) = x(i)
          end do
          do i = 1,m
              ybest(i) = y(i)
          end do
      end if

C     ===============================================================
C     END OF THE OUTER ITERATION
C     ===============================================================

      go to 100

C     ===============================================================
C     END OF THE MAIN LOOP
C     ===============================================================

C     ===============================================================
C     RETURN
C     ===============================================================

 900  continue

C     SET X, Y AND F WITH THE BEST SOLUTION, ITS ASSOCIATED LAGRANGE
C     MULTIPLIERS AND ITS FUNCTIONAL VALUE

      f = fbest

      do i = 1,n
          x(i) = xbest(i)
      end do

      do i = 1,m
          y(i) = ybest(i)
      end do

C     WRITE OUTPUT STATUS

      if ( output .ge. 0 ) then

          if ( flag .eq. 0 ) then
              write(*, fmt=9600) flag,epsdsn
              write(10,fmt=9600) flag,epsdsn

          else if ( flag .eq. 1 ) then
              write(*, fmt=9610) flag,epsden
              write(10,fmt=9610) flag,epsden

          else if ( flag .eq. 6 ) then
              write(*, fmt=9620) flag,maxoutit
              write(10,fmt=9620) flag,maxoutit

          else if ( flag .eq. 7 ) then
              write(*, fmt=9630) flag,maxinnit
              write(10,fmt=9630) flag,maxinnit

          else if ( flag .eq. 8 ) then
              write(*, fmt=9640) flag,maxfc
              write(10,fmt=9640) flag,maxfc

          else if ( flag .eq. -1 ) then
              write(*, fmt=9650) flag
              write(10,fmt=9650) flag

          else if ( flag .eq. -2 ) then
              write(*, fmt=9660) flag
              write(10,fmt=9660) flag

          else if ( flag .eq. -3 ) then
              write(*, fmt=9670) flag
              write(10,fmt=9670) flag

          end if
      end if

C     WRITE STATISTICS

      if ( output .ge. 0 ) then
          write(*, fmt=9700) f,dsupn,sqrt(deucn2),amax,flag
          write(*, fmt=9710) outit,fcnt,gcnt,hcnt,innit

          write(10,fmt=9700) f,dsupn,sqrt(deucn2),amax,flag
          write(10,fmt=9710) outit,fcnt,gcnt,hcnt,innit
      end if

      return

C     NON-EXECUTABLE STATEMENTS

 9010 format(/,1X,'INEXACT VARIABLE METRIC METHOD')
 9020 format(/,1X,'OUTIT= ',I10,' F= ',1P,D14.7,' DSUPN= ',1P,D14.7,
     +            ' AMAX= ',1P,D14.7)

 9200 format(/,6X,'INNER SOLVER (COMPUTING THE SEARCH DIRECTION FOR ',
     +             'THE NEXT ITERATION)')
 9201 format(  6X,'(The function is linear, using Gk = I)')
 9202 format(  6X,'(Using Gk = the true Hessian)')
 9203 format(  6X,'(Using Gk = (1 / lambda) I, where lambda = ',
     +         1P,D14.7,' is the spectral',
     +       /,6X,'steplenght)')
 9205 format(  6X,'INNIT = ',I6,' -L =',1P,D14.7,' NLPSUPN = ',
     +         1P,D14.7)
 9210 format(/,6X,'INNER ITERATION: ',I6,
     +       /,6X,'Current y (first ',I6,' components over a total of ',
     +         I6,'):',
     +       /,1(6X,6(1X,1PD11.4)))
 9240 format(  6X,'Functional value                = ',1PD14.7,
     +       /,6X,'Projected gradient sup-norm     = ',1PD14.7)
 9250 format(  6X,'Associated primal point (first ',I6,' components ',
     +            'over a total of ',I6,'):',
     +       /,1(6X,6(1X,1PD11.4)))
 9260 format(  6X,'Maximum primal step             = ',1PD14.7,
     +       /,6X,'Including constant of reduction = ',1PD14.7,
     +       /,6X,'Q(redc d)                       = ',1PD14.7)
 9270 format(  6X,'Precognized stooping criterion satisfied')
 9280 format(  6X,'Projected gradient norm stooping criterion ',
     +            'satisfied')
 9290 format(  6X,'Leaving the face')
 9300 format(  6X,'Remaining into the face')
 9310 format(  6X,'dy direction (first ',I6, ' components): ',
     +       /,1(6X,6(1X,1PD11.4)))
 9320 format(  6X,'Maximum step: ',1PD14.7,' Variable index: ',I6)
 9330 format(  6X,'The one-dimensional minimizer does not exist')
 9340 format(  6X,'One-dimensional minimizer: ',1PD14.7)

 9500 format(/,6X,'LINE SEARCH')
 9510 format(  6X,'Step = ',1P,D14.7,' F = ',1P,D14.7)
 9520 format(  6X,'Armijo-like stooping criterion satisfied')
 9600 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'convergence with the sup-norm of the search ',
     +       /,1X,'direction smaller than ',1P,D14.7,
     +            ' (parameter epsdsn)')
 9610 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'convergence with the euclidian-norm of the ',
     +            'search',
     +       /,1X,'direction smaller than ',1P,D14.7,
     +            ' (parameter epsden)')
 9620 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'too many outer iterations ',
     +       /,1X,'(more than parameter maxoutit = ',I6,')')
 9630 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'too many inner iterations ',
     +       /,1X,'(more than parameter maxinnit = ',I6,')')
 9640 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'too many functional evaluations ',
     +       /,1X,'(more than parameter maxfc = ',I6,')')
 9650 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'Error in evalf subroutine')
 9660 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'Error in evalg subroutine')
 9670 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'Error in evalH subroutine')
 9700 format(/,1X,'F        = ',   1P,D17.10,
     +       /,1X,'DSUPNORM = ',1X,1P,D16.10,
     +       /,1X,'DEUCNORM = ',1X,1P,D16.10,
     +       /,1X,'AMAX     = ',1X,1P,D16.10,
     +       /,1X,'FLAG     = ',15X,I2)
 9710 format(/,1X,'OUTIT    = ',1X,I16,
     +       /,1X,'FCNT     = ',1X,I16,
     +       /,1X,'GCNT     = ',1X,I16,
     +       /,1X,'HCNT     = ',1X,I16,
     +       /,1X,'INNIT    = ',1X,I16)

      end

c     ******************************************************
c     ******************************************************

      subroutine evallnlispg(m,n,y,g,arow,acol,aval,anze,bbar,lambda,d,
     +                       amax,redc,L,Q,nL,w1,w2,w3,w4,theta,infty)

C     SCALAR ARGUMENTS
      integer m,n,anze
      double precision lambda,amax,redc,L,Q,theta,infty

C     ARRAY ARGUMENTS
      integer arow(anze),acol(anze)
      double precision y(m),g(n),aval(anze),d(n),bbar(m),nL(m),w1(n),
     +        w2(m),w3(n),w4(n)

C     LOCAL SCALARS
      integer i,j
      double precision tmp

C     ---------------------------------------------------------------
C     Computing corresponding primal point d = - lambda ( gk + A^t y )
C     ---------------------------------------------------------------

C     w1 = - ( gk + A^t y )

      do j = 1,n
          w1(j) = - g(j)
      end do

      do i = 1,anze
          w1(acol(i)) = w1(acol(i)) - aval(i) * y(arow(i))
      end do

C     d = - lambda w1

      do i = 1,n
          d(i) = lambda * w1(i)
      end do

C     ---------------------------------------------------------------
C     Computing maximum step to mantain feasibility of xk + d
C     in the primal feasible set A( xk + alpha d ) <= b
C     ---------------------------------------------------------------

C     w2  = A d

      do i = 1,m
          w2(i) = 0.0d0
      end do

      do i = 1,anze
          w2(arow(i)) = w2(arow(i)) + aval(i) * d(acol(i))
      end do

C     amax = min_{i | [A d]_i > 0} { [bbar]_i / [A d]_i }

      amax = infty
      do i = 1,m
          if ( w2(i) .gt. 0.0d0 ) then
              amax = min( amax, bbar(i) / w2(i) )
          end if
      end do

      if ( amax .lt. 0.0d0 ) then
          amax = infty
          do i = 1,m
              if ( w2(i) .gt. 0.0d0 ) then
                  amax = min( amax, bbar(i) / w2(i) )
                  write(*,*) amax,bbar(i),w2(i)
              end if
          end do
          stop
      end if

C     ---------------------------------------------------------------
C     Computing the reduction constant to warranty that xk + redc * d 
C     is an (enough) interior point
C     ---------------------------------------------------------------

      redc = min( 1.0d0, theta * amax )

C     ---------------------------------------------------------------
C     Computing Lk(d,y) and Qk(redc d) 
C     ---------------------------------------------------------------

C     Computing 0.5 d^t Hk d = - d^t (gk + A^t y) = d^t w1

      tmp = 0.0d0
      do i = 1,n
          tmp = tmp + d(i) * w1(i)
      end do

      L = 0.5d0 * tmp
      Q = 0.5d0 * redc ** 2 * tmp

C     Computing gk^t d

      tmp = 0.0d0
      do i = 1,n
          tmp = tmp + g(i) * d(i)
      end do

      L = L + tmp
      Q = Q + redc * tmp

C     Computing y^t ( A xk - b + A d  ) = y^t ( w2 - bbar )

      tmp = 0.0d0
      do i = 1,m
          tmp = tmp + y(i) * ( w2(i) - bbar(i) )
      end do

      L = L + tmp

C     ---------------------------------------------------------------
C     Computing - \nabla L(d,y) = A lambda (gk + A^t y) + bbar = bbar - A d
C     ---------------------------------------------------------------

C     nL is the antigradient of L. As we are minimizing -L, nL is the
C     positive gradient of the function being minimized.

      do i = 1,m
          nL(i) = bbar(i) - w2(i)
      end do

      end

c ****************************************************
c ****************************************************

c     Last update: March 26th, 2003

c     See report of last modifications at the end of this file.

      subroutine gencan(n,x,l,u,epsgpen,epsgpsn,epsnfp,maxitnfp,
     *maxitngp,fmin,maxit,maxfc,udelta0,ucgmaxit,cgscre,cggpnf,cgepsi,
     *cgepsf,epsnqmp,maxitnqmp,nearlyq,gtype,htvtype,trtype,iprint,
     *f,g,gpeucn2,gpsupn,iter,fcnt,gcnt,cgcnt,spgiter,spgfcnt,
     *tniter,tnfcnt,tnstpcnt,tnintcnt,tnexgcnt,tnexbcnt,tnintfe,tnexgfe,
     *tnexbfe,flag,s,y,d,ind,w,eta,delmin,lammax,lammin,theta,gamma,
     *beta,sigma1,sigma2,nint,next,mininterp,ncomp,sterel,steabs,epsrel,
     *epsabs,infty)

      logical nearlyq

      integer n,maxitnfp,maxitngp,maxit,maxfc,maxitnqmp,ucgmaxit,
     *cgscre,gtype,htvtype,trtype,iprint,iter,fcnt,gcnt,cgcnt,
     *spgiter,spgfcnt,tniter,tnfcnt,tnstpcnt,tnintcnt,tnexgcnt,tnexbcnt,
     *tnintfe,tnexgfe,tnexbfe,flag,ind(n),mininterp,ncomp

      double precision x(n),l(n),u(n),epsgpen,epsgpsn,epsnfp,fmin,
     *udelta0,cggpnf,cgepsi,cgepsf,epsnqmp,f,g(n),gpeucn2,gpsupn,s(n),
     *y(n),d(n),w(5*n),eta,delmin,lammax,lammin,theta,gamma,beta,sigma1,
     *sigma2,nint,next,sterel,steabs,epsrel,epsabs,infty

c     Solves the box-constrained minimization problem
c
c                  Minimize f(x)
c                  subject to l \leq x \leq u
c     
c     using a method described in 
c
c     E. G. Birgin and J. M. Martinez, "Large-scale active-set 
c     box-constrained optimization method with spectral projected 
c     gradients", Computational Optimization and Applications, 
c     23, pp. 101-125, 2002.  
c
c     Subroutines evalf and evalg must be supplied by the user to 
c     evaluate the function f and its gradient, respectively. The 
c     calling sequences are
c
c                  call evalf(n, x, f, inform)
c
c                  call evalg(n, x, g, inform)
c
c     where x is the point where the function (the gradient) must be 
c     evaluated, n is the number of variables and f (g) is the 
c     functional value (the gradient vector). The real parameters
c     x, f, g must be double precision.
c
c     A subroutine evalhd to compute the Hessian times vector products
c     is optional. If this subroutine is not provided an incremental
c     quotients version will be used instead. The calling sequence of
c     this subroutine should be
c
c                  call evalhd(nind,ind,n,x,u,hu,inform)
c
c     where x is the point where the approx-Hessian is being considered, 
c     n is the number of variables, u is the vector which should be 
c     multiplied by the approx-Hessian H and hu is the vector where the 
c     product should be placed. The information about the matrix H must 
c     be passed to evalhd by means of common declarations. The necessary 
c     computations must be done in evalg. The real parameters x, u, hu
c     must be double precision.
c
c     This subroutine must be coded by the user, taking into account 
c     that n is the number of variables of the problem and that hu must 
c     be the product H u. Moreover, you must assume, when you code evalhd,
c     that only nind components of u are nonnull and that ind is the set
c     of indices of those components. In other words, you must write
c     evalhd in such a way that hu is the vector whose i-th entry is
c 
c               hu(i) = \Sum_{j=1}^{nind} H_{i,ind(j)} u_ind(j)
c
c     Moreover, the only components of hu that you need to compute are
c     those which corresponds to the indices ind(1),...,ind(nind). 
c     However, observe that you must assume that, in u, the whole 
c     vector is present, with its n components, even the zeroes. So, 
c     if you decide to code evalhd without taking into account the 
c     presence of ind and nind, you can do it. A final observation: 
c     probably, if nind is close to n, it is not worthwhile to use ind, 
c     due to the cost of accessing the correct indices. If you want, 
c     you can test, within your evalhd, if (say) nind > n/2, and, in 
c     this case use a straightforward scalar product for the components 
c     of hu.
c
c     Example: Suppose that the matrix H is full. The main steps of 
c     evalhd could be:
c
c          do i= 1, nind
c              indi= ind(i)
c              hu(indi)= 0.0d0
c              do j= 1, nind
c                  indj= ind(j)
c                  hu(indi)= hu(indi) + H(indi,indj) * u(indj)
c              end do
c          end do
c
c
c     On Entry
c
c     n    integer
c          the space dimension
c
c     x    double precision x(n)
c          initial estimate to the solution
c
c     l    double precision l(n)
c          lower bounds
c
c     u    double precision u(n)
c          upper bounds
c
c     epsgpen double precision
c          small positive number for declaring convergence when the 
c          euclidian norm of the projected gradient is less than
c          or equal to epsgpen
c
c          RECOMMENDED: epsgpen = 1.0d-5
c
c     epsgpsn double precision
c          small positive number for declaring convergence when the 
c          infinite norm of the projected gradient is less than
c          or equal to epsgpsn
c
c          RECOMMENDED: epsgpsn = 1.0d-5
c
c     epsnfp double precision
c          ``lack of enough progress'' measure. The algorithm stops by
c          ``lack of enough progress'' when f(x_k) - f(x_{k+1}) <=  
c          epsnfp * max { f(x_j)-f(x_{j+1}, j<k} during maxitnfp 
c          consecutive iterations. This stopping criterion may be 
c          inhibited setting epsnfp = 0. We recommend, preliminary, to 
c          set epsnfp = 0.01 and maxitnfp = 5
c
c          RECOMMENDED: epsnfp = 1.0d-2
c
c     maxitnfp integer
c          see the meaning of epsnfp, above
c
c          RECOMMENDED: maxitnfp = 5
c
c     maxitngp integer
c          If the order of the euclidian-norm of the continuous projected 
c          gradient did not change during maxitngp consecutive iterations 
c          the execution stops. Recommended: maxitngp= 10. In any case 
c          maxitngp must be greater than or equal to 1
c
c          RECOMMENDED: maxitngp = 10
c
c     fmin double precision
c          function value for the stopping criteria f <= fmin
c
c          RECOMMENDED: fmin = -1.0d+99 (inhibited)
c
c     maxit integer
c          maximum number of iterations allowed
c
c          RECOMMENDED: maxit = 1000
c
c     maxfc integer
c          maximum number of funtion evaluations allowed
c
c          RECOMMENDED: maxfc = 5000
c
c     udelta0 double precision
c          initial trust-region radius. Default max{0.1||x||,0.1} is set
c          if you set udelta0 < 0. Otherwise, the parameters udelta0
c          will be the ones set by the user.
c
c          RECOMMENDED: udelta0 = -1
c
c     ucgmaxit integer
c          maximum number of iterations allowed for the cg subalgorithm
c
c          Default values for this parameter and the previous one are 
c          0.1 and 10 * log (number of free variables). Default values 
c          are taken if you set ucgeps < 0 and ucgmaxit < 0, 
c          respectively. Otherwise, the parameters ucgeps and ucgmaxit 
c          will be the ones set by the user
c
c          RECOMMENDED: ucgmaxit = -1
c
c     cgscre, cggpnf double precision
c          cgscre means cunjugate gradient stopping criterion relation, and
c          cggpnf means conjugate gradients projected gradient final norm.
c          Both are related to a stopping criterion of conjugate gradients.
c          This stopping criterion depends on the norm of the residual
c          of the linear system. The norm of the this residual should be
c          less or equal than a ``small''quantity which decreases as we are
c          approximating the solution of the minimization problem (near the
c          solution, better the truncated-Newton direction we aim). Then, the
c          log of the required precision requested to conjugate gradient has 
c          a linear dependence on the log of the norm of the projected 
c          gradient. This linear relation uses the squared euclidian-norm
c          of the projected gradient if cgscre = 1 and uses the sup-norm if
c          cgscre = 2. In adition, the precision required to CG is equal to
c          cgepsi (conjugate gradient initial epsilon) at x0 and cgepsf
c          (conjugate gradient final epsilon) when the euclidian- or sup-norm
c          of the projected gradient is equal to cggpnf (conjugate gradients 
c          projected gradient final norm) which is an estimation of the value
c          of the euclidian- or sup-norm of the projected gradient at the
c          solution.
c
c          RECOMMENDED: cgscre = 1, cggpnf = epsgpen; or
c                       cgscre = 2, cggpnf = epsgpsn.
c
c     cgepsi, cgepsf double precision
c          small positive numbers for declaring convergence of the
c          conjugate gradient subalgorithm when 
c          ||r||_2 < cgeps * ||rhs||_2, where r is the residual and
c          rhs is the right hand side of the linear system, i.e., cg
c          stops when the relative error of the solution is smaller 
c          that cgeps. 
c
c          cgeps varies from cgepsi to cgepsf in such a way that, depending
c          on cgscre (see above),
c
c          i) log10(cgeps^2) depends linearly on log10(||g_P(x)||_2^2)
c          which varies from ||g_P(x_0)||_2^2 to epsgpen^2; or
c
c          ii) log10(cgeps) depends linearly on log10(||g_P(x)||_inf)
c          which varies from ||g_P(x_0)||_inf to epsgpsn.
c
c          RECOMMENDED: cgepsi = 1.0d-1, cgepsf = 1.0d-5
c
c     epsnqmp double precision
c          see below
c
c     maxitnqmp integer
c          This and the previous one parameter are used for a stopping
c          criterion of the conjugate gradient subalgorithm. If the
c          progress in the quadratic model is less or equal than a 
c          fraction of the best progress ( epsnqmp * bestprog ) during
c          maxitnqmp consecutive iterations then CG is stopped by not
c          enough progress of the quadratic model.
c
c          RECOMMENDED: epsnqmp = 1.0d-4, maxitnqmp = 5
c
c     nearlyq logical
c          if function f is (nearly) quadratic, use the option 
c          nearlyq = TRUE. Otherwise, keep the default option.
c
c          if, at an iteration of CG we find a direction d such
c          that d^T H d <= 0 then we take the following decision:
c
c          (i) if nearlyq = TRUE then take direction d and try to 
c          go to the boundary chosing the best point among the two
c          point at the boundary and the current point. 
c
c          (ii) if nearlyq = FALSE then we stop at the current point.
c
c          RECOMMENDED: nearlyq = .FALSE.
c
c     gtype integer
c          type of gradient calculation
c          gtype = 0 means user suplied evalg subroutine,
c          gtype = 1 means central diference approximation.
c
c          RECOMMENDED: gtype = 0
c
c          (provided you have the evalg subroutine)
c
c     htvtype integer
c          type of gradient calculation
c          htvtype = 0 means user suplied evalhd subroutine,
c          htvtype = 1 means incremental quotients approximation.
c
c          RECOMMENDED: htvtype = 1
c
c          (you take some risk using this option but, unless you have
c          a good evalhd subroutine, incremental quotients is a very 
c          cheap option)
c
c     trtype integer
c          type of trust-region radius
c          trtype = 0 means 2-norm trust-region
c          trtype = 1 means infinite-norm trust-region
c
c          RECOMMENDED: trtype = 0
c
c     iprint integer
c          commands printing. Nothing is printed if iprint < 0.
c          If iprint >= 0, gencan iterations information is printed.
c          If iprint >= 1, line searches and conjugate gradient information 
c          is printed.
c
c          RECOMMENDED: iprint = 1
c
c     s,y,d,w double precision s(n),y(n),d(n),w(5*n)
c          working vectors
c
c     ind  integer ind(n)
c          working vector
c
c     eta  double precision
c          constant for deciding abandon the current face or not
c          We abandon the current face if the norm of the internal 
c          gradient (here, internal components of the continuous 
c          projected gradient) is smaller than (1-eta) times the 
c          norm of the continuous projected gradient. Using eta=0.9 
c          is a rather conservative strategy in the sense that 
c          internal iterations are preferred over SPG iterations. 
c
c          RECOMMENDED: eta = 0.9
c
c     delmin double precision
c          minimum ``trust region'' to compute the Truncated Newton
c          direction
c
c          RECOMMENDED: delmin = 0.1
c
c     lammin, lammax double precision
c          The spectral steplength, called lambda, is projected
c          inside the box [lammin,lammax] 
c
c          RECOMMENDED: lammin = 10^{-10} and lammax = 10^{10}
c 
c     theta double precision
c          constant for the angle condition, i.e., at iteration k
c          we need a direction d_k such that
c          <g_k,d_k> <= -theta ||g||_2 ||d_k||_2,
c          where g_k is \nabla f(x_k)
c
c          RECOMMENDED: theta = 10^{-6}
c
c     gamma double precision
c          constant for the Armijo crtierion
c          f(x + alpha d) <= f(x) + gamma * alpha * <\nabla f(x),d>
c
c          RECOMMENDED: gamma = 10^{-4}
c
c     beta double precision
c          constant for the beta condition
c          <d_k, g(x_k + d_k)>  .ge.  beta * <d_k,g_k>
c          if (x_k + d_k) satisfies the Armijo condition but does not
c          satisfy the beta condition then the point is accepted, but
c          if it satisfied the Armijo condition and also satisfies the
c          beta condition then we know that there is the possibility
c          for a succesful extrapolation
c
c          RECOMMENDED: beta = 0.5
c
c     sigma1, sigma2 double precision
c          constant for the safeguarded interpolation
c          if alpha_new \notin [sigma1, sigma*alpha] then we take
c          alpha_new = alpha / nint
c
c          RECOMMENDED: sigma1 = 0.1 and sigma2 = 0.9
c
c     nint double precision
c          constant for the interpolation. See the description of
c          sigma1 and sigma2 above. Sometimes we take as a new trial
c          step the previous one divided by nint
c
c          RECOMMENDED: nint = 2.0
c
c     next double precision
c          constant for the extrapolation
c          when extrapolating we try alpha_new = alpha * next
c
c          RECOMMENDED: next = 2.0
c
c     mininterp integer
c          constant for testing if, after having made at least mininterp
c          interpolations, the steplength is too small. In that case 
c          failure of the line search is declared (may be the direction
c          is not a descent direction due to an error in the gradient 
c          calculations) 
c
c          RECOMMENDED: mininterp = 4 
c
c          (use mininterp > maxfc for inhibit this stopping criterion)
c
c     ncomp integer
c          this constant is just for printing. In a detailed printing
c          option, ncomp component of the current point will be printed
c
c          RECOMMENDED: ncomp = 5
c
c     sterel, steabs double precision
c          this constants mean a ``relative small number'' and ``an 
c          absolute small number'' for the increments in finite
c          difference approximations of derivatives
c
c          RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} 
c
c     epsrel, epsabs, infty  double precision
c          this constants mean a ``relative small number'', ``an 
c          absolute small number'', and ``infinite or a very big
c          number''. Basically, a quantity A is considered negligeble 
c          with respect to another quantity B if 
c          |A| < max ( epsrel * |B|, epsabs ) 
c
c          RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20} and 
c          infty = 10^{+20}
c
c
c     On Return
c
c     x    double precision x(n)
c          final estimation to the solution
c
c     f    double precision
c          function value at the final estimation 
c
c     g    double precision g(n)
c          gradient at the final estimation
c
c     gpeucn2 double precision
c          squared 2-norm of the continuous projected
c          gradient g_p at the final estimation (||g_p||_2^2)
c
c     gpsupn double precision
c          ||g_p||_inf at the final estimation
c
c     iter integer
c          number of iterations
c
c     fcnt integer
c          number of function evaluations   
c
c     gcnt integer
c          number of gradient evaluations   
c
c     cgcnt integer
c          number of conjugate gradient iterations   
c
c     spgiter integer
c          number of SPG iterations
c
c     spgfcnt integer
c          number of function evaluations in SPG-directions line searches
c
c     tniter integer
c          number of Truncated Newton iterations
c
c     tnfcnt integer
c          number of function evaluations in TN-directions line searches
c
c     tnintcnt integer
c          number of times a backtracking in a TN-drection was needed
c
c     tnexgcnt integer
c          number of times an extrapolation in a TN-direction was
c          succesful in decreass the function value
c
c     tnexbcnt integer
c          number of times an extrapolation was aborted in the first
c          extrapolated point by augment of the function value
c
c     flag integer
c          This output parameter tells what happened in this 
c          subroutine, according to the following conventions:
c 
c          0= convergence with small euclidian-norm of the
c             projected gradient (smaller than epsgpen);
c
c          1= convergence with small infinite-norm of the
c             projected gradient (smaller than epsgpsn);
c
c          2= the algorithm stopped by ``lack of enough progress'',
c             that means that f(x_k) - f(x_{k+1}) <=  
c             epsnfp * max { f(x_j)-f(x_{j+1}, j<k} during maxitnfp 
c             consecutive iterations;
c
c          3= the algorithm stopped because the order of the euclidian-
c             norm of the continuous projected gradient did not change 
c             during maxitngp consecutive iterations. Probably, we 
c             are asking for an exagerately small norm of continuous 
c             projected gradient for declaring convergence;
c
c          4= the algorithm stopped because the functional value
c             is very small (f <= fmin);
c
c          6= too small step in a line search. After having made at 
c             least mininterp interpolations, the steplength becames 
c             small. ``small steplength'' means that we are at point
c             x with direction d and step alpha, and 
c
c             alpha * ||d||_infty < max(epsabs, epsrel * ||x||_infty). 
c 
c             In that case failure of the line search is declared 
c             (may be the direction is not a descent direction 
c             due to an error in the gradient calculations). Use
c             mininterp > maxfc for inhibit this criterion;
c
c          7= it was achieved the maximum allowed number of
c             iterations (maxit);
c
c          8= it was achieved the maximum allowed number of
c             function evaluations (maxfc);
c
c          9= error in evalf subroutine;
c
c         10= error in evalg subroutine;
c
c         11= error in evalhd subroutine.

      character*3 ittype
      integer i,itnfp,itngp,cgflag,lsflag,nind,nprint,cgmaxit,cgiter,
     *fcntprev,tnintprev,tnexgprev,tnexbprev,rbdtype,rbdind,inform
      double precision lambda,gieucn2,gpi,fprev,sts,sty,epsgpen2,delta,
     *currprog,bestprog,gexp,gexpprev,cgeps,xnorm,ometa2,amax,amaxx,
     *acgeps,bcgeps,kappa,gpeucn20,gpsupn0

c     =======================================================
c     Initialization
c     =======================================================

c     Set some initial values:

c     counters,
      iter= 0
      fcnt= 0
      gcnt= 0
      cgcnt= 0

      spgiter= 0
      spgfcnt= 0

      tniter= 0
      tnfcnt= 0

      tnstpcnt= 0
      tnintcnt= 0
      tnexgcnt= 0
      tnexbcnt= 0

      tnintfe= 0
      tnexgfe= 0
      tnexbfe= 0

c     just for printing,
      nprint= min0(n,ncomp)

c     for testing convergence,
      epsgpen2= epsgpen**2

c     for testing wether abandon the current face or not,
c     (ometa2 means '(one minus eta) squared')
      ometa2= (1.0d0-eta)**2

c     for testing progress in f,
      fprev= infty
      bestprog= 0.0d0
      itnfp= 0

c     for testing progress in the order of the gradient norm, and   
      gexpprev= infty
      itngp= 0

c     Print problem information

      if(iprint.ge.0) then
          write(*,977) n
          write(*,978) nprint,(l(i),i=1,nprint)
          write(*,979) nprint,(u(i),i=1,nprint)
          write(*,980) nprint,(x(i),i=1,nprint)

          write(10,977) n
          write(10,978) nprint,(l(i),i=1,nprint)
          write(10,979) nprint,(u(i),i=1,nprint)
          write(10,980) nprint,(x(i),i=1,nprint)
      end if

c     Project initial guess. If the initial guess is infeasible, 
c     projection puts it in the box.

      do i= 1, n
          x(i)= dmax1(l(i),dmin1(x(i),u(i)))
      end do

c     Compute x euclidian-norm

      xnorm= 0.0d0
      do i= 1, n
          xnorm= xnorm + x(i)**2
      end do
      xnorm= dsqrt(xnorm)

c     Compute function and gradient at the initial point

      call evalq(n,x,f,inform)
      fcnt= fcnt + 1

      if (inform.ne.0) then
          flag = 9

          if (iprint.ge.0) then
              write(*,999)  flag
              write(10,999) flag
          end if

          go to 500
      end if

      if (gtype.eq.0) then
          call evalnq(n,x,g,inform)
      else if (gtype.eq.1) then
          call evalgdiff(n,x,g,sterel,steabs,inform)
      end if
      gcnt= gcnt + 1

      if (inform.ne.0) then
          flag = 10

          if (iprint.ge.0) then
              write(*,1000)  flag
              write(10,1000) flag
          end if

          go to 500
      end if

c     Compute continuous project gradient infinite- and euclidian-
c     norm, internal gradient euclidian-norm, and store in nind the
c     number of free variables and in array p their identifiers.

      nind= 0
      gpsupn= 0.0d0
      gpeucn2= 0.0d0
      gieucn2= 0.0d0
      do i= 1, n
          gpi= dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i)
          gpsupn= dmax1(gpsupn, dabs(gpi))
          gpeucn2= gpeucn2 + gpi**2
          if (x(i).gt.l(i).and.x(i).lt.u(i)) then
              gieucn2= gieucn2 + gpi**2
              nind= nind + 1
              ind(nind)= i
          end if
      end do

c     Compute a linear relation between gpeucn2 and cgeps2, i.e.,
c     scalars a and b such that 
c
c         a * log10(||g_P(x_0)||_2^2) + b = log10(cgeps_0^2) and
c
c         a * log10(||g_P(x_f)||_2^2) + b = log10(cgeps_f^2),
c
c     where cgeps_0 and cgeps_f are provided. Note that if 
c     cgeps_0 is equal to cgeps_f then cgeps will be allways 
c     equal to cgeps_0 and cgeps_f.

c     We introduce now a linear relation between gpsupn and cgeps also.

      if (cgscre.eq.1) then
          acgeps= 2.0d0*dlog10(cgepsf/cgepsi)/dlog10(cggpnf**2/gpeucn2)
          bcgeps= 2.0d0*dlog10(cgepsi)-acgeps*dlog10(gpeucn2)
      else ! if (cgscre.eq.2) then
          acgeps= dlog10(cgepsf/cgepsi)/dlog10(cggpnf/gpsupn)
          bcgeps= dlog10(cgepsi)-acgeps*dlog10(gpsupn)
      end if 

c     And it will be used for the linear relation of cgmaxit

      gpsupn0= gpsupn
      gpeucn20= gpeucn2

c     Print initial information

      if(iprint.ge.0) then
          write(*,981) iter
          write(*,985) nprint,(x(i),i=1,nprint)
          write(*,986) nprint,(g(i),i=1,nprint)
          write(*,987) nprint,
     *    (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint)
          write(*,988) min0(nprint,nind),nind,
     *    (ind(i),i=1,min0(nprint,nind))
          write(*,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n,
     *    spgiter,tniter,fcnt,gcnt,cgcnt

          write(10,981) iter
          write(10,985) nprint,(x(i),i=1,nprint)
          write(10,986) nprint,(g(i),i=1,nprint)
          write(10,987) nprint,
     *    (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint)
          write(10,988) min0(nprint,nind),nind,
     *    (ind(i),i=1,min0(nprint,nind))
          write(10,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n,
     *    spgiter,tniter,fcnt,gcnt,cgcnt
      end if

c     Test whether the initial functional value is very small

      if (f.le.fmin) then
          flag= 4

          if (iprint.ge.0) then
              write(*,994)  flag,fmin
              write(10,994) flag,fmin
          end if

          go to 500
      end if

c     Test whether the number of functional evaluations is 
c     exhausted (i.e., maxfc = 1)

      if (fcnt.ge.maxfc) then
          flag= 8

          if (iprint.ge.0) then
              write(*,998)  flag,maxfc
              write(10,998) flag,maxfc
          end if

          go to 500
      end if

c     =======================================================
c     Main loop
c     =======================================================

 100  continue

c     =======================================================
c     Test stopping criteria
c     =======================================================

c     Test whether the continuous projected gradient euclidian-norm
c     is small enough to declare convergence

      if (gpeucn2.le.epsgpen2) then
          flag= 0

          if (iprint.ge.0) then
              write(*,990)  flag,epsgpen
              write(10,990) flag,epsgpen
          end if

          go to 500
      end if

c     Test whether the continuous projected gradient infinite-norm
c     is small enough to declare convergence

      if (gpsupn.le.epsgpsn) then
          flag= 1

          if (iprint.ge.0) then
              write(*,991)  flag,epsgpsn
              write(10,991) flag,epsgpsn
          end if

          go to 500
      end if

c     Test whether we performed many iterations without good progress
c     of the functional value

      currprog= fprev - f
      bestprog= dmax1(currprog, bestprog)

      if (currprog.le.epsnfp*bestprog) then

          itnfp= itnfp + 1

          if (itnfp.ge.maxitnfp) then
              flag= 2

              if (iprint.ge.0) then
                  write(*,992)  flag,epsnfp,maxitnfp
                  write(10,992) flag,epsnfp,maxitnfp
              end if

              go to 500
          endif

      else
          itnfp= 0
      endif

c     Test whether we have performed many iterations without good 
c     reduction of the euclidian-norm of the projected gradient

      gexp= dlog10(gpeucn2)
   
      if (gexp.ge.gexpprev) then

          itngp= itngp + 1

          if(itngp.ge.maxitngp) then
              flag= 3

              if (iprint.ge.0) then
                  write(*,993)  flag,maxitngp
                  write(10,993) flag,maxitngp
              end if

              go to 500
          endif

      else
          itngp= 0
      endif

      gexpprev= gexp

c     Test whether the number of iterations is exhausted

      if (iter.ge.maxit) then
          flag= 7

          if (iprint.ge.0) then
              write(*,997)  flag,maxit
              write(10,997) flag,maxit
          end if

          go to 500
      end if

c     Note that the stopping criteria related to small functional
c     values ( f <= fmin ) and number of functional evaluations
c     are detected inside the line searches and do not need to
c     be tested here

c     =======================================================
c     The stopping criteria were not satisfied, a new 
c     iteration will be made
c     =======================================================

      iter= iter + 1

c     =======================================================
c     Save current values, f, x and g
c     =======================================================

      fprev= f

      do i= 1, n
          s(i)= x(i)
          y(i)= g(i)
      end do

c     =======================================================
c     Compute new iterate
c     =======================================================

c     We abandon the current face if the norm of the internal gradient
c     (here, internal components of the continuous projected gradient)
c     is smaller than (1-eta) times the norm of the continuous projected 
c     gradient. Using eta=0.9 is a rather conservative strategy in the 
c     sense that internal iterations are preferred over SPG iterations. 
c     Replace eta=0.9 by other tolerance in (0,1) if you find it convenient. 

      if (gieucn2.le.ometa2*gpeucn2) then

c         ===================================================
c         Some constraints should be abandoned. Compute
c         the new iterate using an SPG iteration
c         ===================================================

          ittype= 'SPG'
          spgiter= spgiter + 1

c         Compute spectral steplength

          if (iter.eq.1.or.sty.le.0.0d0) then
              lambda= dmax1(1.0d0,xnorm)/dsqrt(gpeucn2)
          else
              lambda= sts/sty
          end if
          lambda= dmin1(lammax,dmax1(lammin,lambda))

c         Perform a line search with safeguarded quadratic 
c         interpolation along the direction of the spectral 
c         continuous projected gradient

          fcntprev= fcnt

          call spgls(n,x,f,g,l,u,lambda,fmin,maxfc,iprint,fcnt,lsflag,
     *    w(1),w(n+1),gamma,sigma1,sigma2,nint,mininterp,sterel,steabs,
     *    epsrel,epsabs,infty) 

          spgfcnt= spgfcnt + (fcnt-fcntprev) 

          if (lsflag.eq.9) then
              flag = 9

              if (iprint.ge.0) then
                  write(*,999)  flag
                  write(10,999) flag
              end if

              go to 500
          end if

c         Compute the gradient of the new iterate

          if (gtype.eq.0) then
              call evalnq(n,x,g,inform)
          else if (gtype.eq.1) then
              call evalgdiff(n,x,g,sterel,steabs,inform)
          end if
          gcnt= gcnt + 1
 
          if (inform.ne.0) then
              flag = 10

              if (iprint.ge.0) then
                  write(*,1000)  flag
                  write(10,1000) flag
              end if

              go to 500
          end if

      else

c         ===================================================
c         The new iterate will belong to the closure of 
c         the current face
c         ===================================================

          ittype= 'TN '
          tniter= tniter + 1

c         Compute trust-region radius

          if (iter.eq.1) then
              if(udelta0.lt.0.d0) then
                  delta= dmax1(delmin,100.0d0*dmax1(1.0d0,xnorm))
              else
                  delta= udelta0
              end if
          else
              delta= dmax1(delmin,10.0d0*dsqrt(sts))
          end if

c         Shrink the point, its gradient and the bounds

          call Gshrink(nind,ind,n,x,g,l,u)

c         Compute the descent direction solving the newtonian 
c         system by conjugate gradients

c         Set conjugate gradient stopping criteria. Default values are 
c         taken if you set ucgeps < 0 and ucgmaxit < 0, respectively. 
c         Otherwise, the parameters cgeps and cgmaxit will be the ones 
c         set by the user.

          if(ucgmaxit.lt.0) then
              if (nearlyq) then
                  cgmaxit= nind
              else
                  if (cgscre.eq.1) then
                      kappa= dlog10(gpeucn2/gpeucn20)/
     *                       dlog10(epsgpen2/gpeucn20)
                  else ! if (cgscre.eq.2) then
                      kappa= dlog10(gpsupn/gpsupn0)/ 
     *                       dlog10(epsgpsn/gpsupn0)
                  end if
                  kappa= dmax1(0.0d0,dmin1(1.0d0,kappa))
                  cgmaxit= int(
     *            (1.0d0-kappa)*dmax1(1.0d0,10.0d0*dlog10(dfloat(nind)))
     *            + kappa*dfloat(nind) )
              end if
c             cgmaxit = 2 * nind
          else
              cgmaxit= ucgmaxit
          end if

          if (cgscre.eq.1) then
              cgeps= dsqrt(10.0d0**(acgeps*dlog10(gpeucn2)+bcgeps))
          else ! if (cgscre.eq.2) then
              cgeps= 10.0d0**(acgeps*dlog10(gpsupn)+bcgeps)
          end if
          cgeps= dmax1(cgepsf,dmin1(cgepsi,cgeps))

c         Call conjugate gradients

          call cg(nind,ind,n,x,g,delta,l,u,cgeps,epsnqmp,maxitnqmp,
     *    cgmaxit,nearlyq,gtype,htvtype,trtype,iprint,ncomp,d,cgiter,
     *    rbdtype,rbdind,cgflag,w(1),w(n+1),w(2*n+1),w(3*n+1),w(4*n+1),
     *    theta,sterel,steabs,epsrel,epsabs,infty)

	  cgcnt= cgcnt + cgiter

          if (cgflag.eq.10) then
              flag = 10

              if (iprint.ge.0) then
                  write(*,1000)  flag
                  write(10,1000) flag
              end if

              go to 500

          else if (cgflag.eq.11) then
              flag = 11

              if (iprint.ge.0) then
                  write(*,1001)  flag
                  write(10,1001) flag
              end if

              go to 500
          end if

c         Compute maximum step

          if (cgflag.eq.2) then
              amax= 1.0d0
          else
              amax= infty
              do i= 1, nind
                  if (d(i).gt.0.0d0.and.u(i).lt.infty) then
                      amaxx= (u(i)-x(i))/d(i)
                      if (amaxx.lt.amax) then
                          amax= amaxx
                          rbdind= i
                          rbdtype= 2
                      end if
                  else if (d(i).lt.0.0d0.and.l(i).gt.-infty) then
                      amaxx= (l(i)-x(i))/d(i)
                      if (amaxx.lt.amax) then
                          amax= amaxx
                          rbdind= i
                          rbdtype= 1
                      end if
                  end if
               end do
          end if

c         Perform the line search

          tnintprev= tnintcnt
          tnexgprev= tnexgcnt
          tnexbprev= tnexbcnt

          fcntprev= fcnt

          call tnls(nind,ind,n,x,f,g,d,l,u,amax,rbdtype,rbdind,fmin,
     *    maxfc,gtype,iprint,fcnt,gcnt,tnintcnt,tnexgcnt,tnexbcnt,
     *    lsflag,w(1),w(n+1),gamma,beta,sigma1,sigma2,nint,next,
     *    mininterp,sterel,steabs,epsrel,epsabs,infty)

          if (lsflag.eq.9) then
              flag = 9

              if (iprint.ge.0) then
                  write(*,999)  flag
                  write(10,999) flag
              end if

              go to 500

          else if (lsflag.eq.10) then
              flag = 10

              if (iprint.ge.0) then
                  write(*,1000)  flag
                  write(10,1000) flag
              end if

              go to 500
          end if

          tnfcnt= tnfcnt + (fcnt-fcntprev)

          if (tnintcnt.gt.tnintprev) then
              tnintfe= tnintfe + (fcnt-fcntprev)
          else if (tnexgcnt.gt.tnexgprev) then
              tnexgfe= tnexgfe + (fcnt-fcntprev)
          else if (tnexbcnt.gt.tnexbprev) then
              tnexbfe= tnexbfe + (fcnt-fcntprev)
          else
              tnstpcnt= tnstpcnt + 1
          end if

c         Expand the point, its gradient and the bounds

          call Gexpand(nind,ind,n,x,g,l,u)

c         If the line search (interpolation) in the Truncated Newton
c         direction stopped due to a very small step (lsflag=6) we 
c         will discard this iteration and force a SPG iteration

c         Note that tnls subroutine was coded in such a way that in
c         case of lsflag=6 termination the subroutine discards all
c         what was made and returns with the same point it started

          if (lsflag.eq.6) then

              if (iprint.ge.0) then
                  write(*,*)  
                  write(*,*)  
     *            '     The previous NT iteration was discarded due to',
     *            '     a termination for very small step in the line ',
     *            '     search. A SPG iteration will be forced now.   '

                  write(10,*)  
                  write(10,*)  
     *            '     The previous NT iteration was discarded due to',
     *            '     a termination for very small step in the line ',
     *            '     search. A SPG iteration will be forced now.   '
              end if

              ittype= 'SPG'
              spgiter= spgiter + 1

c             Compute spectral steplength

              if (iter.eq.1.or.sty.le.0.0d0) then
                  lambda= dmax1(1.0d0,xnorm)/dsqrt(gpeucn2)
              else
                  lambda= sts/sty
              end if
              lambda= dmin1(lammax,dmax1(lammin,lambda))

c             Perform a line search with safeguarded quadratic 
c             interpolation along the direction of the spectral 
c             continuous projected gradient

              fcntprev= fcnt

              call spgls(n,x,f,g,l,u,lambda,fmin,maxfc,iprint,fcnt,
     *        lsflag,w(1),w(n+1),gamma,sigma1,sigma2,nint,mininterp,
     *        sterel,steabs,epsrel,epsabs,infty) 

              spgfcnt= spgfcnt + (fcnt-fcntprev)

              if (lsflag.eq.9) then
                  flag = 9

                  if (iprint.ge.0) then
                      write(*,999)  flag
                      write(10,999) flag
                  end if

                  go to 500
              end if

c             Compute the gradient in the new iterate

              if (gtype.eq.0) then
                  call evalnq(n,x,g,inform)
              else if (gtype.eq.1) then
                  call evalgdiff(n,x,g,sterel,steabs,inform)
              end if
              gcnt= gcnt + 1
 
              if (inform.ne.0) then
                  flag = 10

                  if (iprint.ge.0) then
                      write(*,1000)  flag
                      write(10,1000) flag
                  end if

                  go to 500
              end if

          end if

      end if

c     =======================================================
c     Prepare for the next iteration 
c     =======================================================

c     This adjustment/projection is "por lo que las putas pudiera"

      do i= 1, n
          if (x(i).le.l(i)+dmax1(epsrel*dabs(l(i)),epsabs)) then
              x(i)= l(i)
          else if (x(i).ge.u(i)-dmax1(epsrel*dabs(u(i)),epsabs)) then  
              x(i)= u(i)
          end if
      end do

c     Compute x euclidian-norm

      xnorm= 0.0d0
      do i= 1, n
          xnorm= xnorm + x(i)**2
      end do
      xnorm= dsqrt(xnorm)

c     Compute s = x_{k+1} - x_k, y = g_{k+1} - g_k, <s,s> and <s,y>

      sts= 0.0d0
      sty= 0.0d0
      do i= 1, n
          s(i)= x(i)-s(i)
          y(i)= g(i)-y(i)
          sts= sts + s(i)*s(i)
          sty= sty + s(i)*y(i)
      end do

c     Compute continuous project gradient infinite- and euclidian-
c     norm, internal gradient euclidian-norm, and store in nind the
c     number of free variables and in array p its identifiers.
      
      nind= 0
      gpsupn= 0.0d0
      gpeucn2= 0.0d0
      gieucn2= 0.0d0
      do i= 1, n
          gpi= dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i)
          gpsupn= dmax1(gpsupn, dabs(gpi))
          gpeucn2= gpeucn2 + gpi**2
          if (x(i).gt.l(i).and.x(i).lt.u(i)) then
              gieucn2= gieucn2 + gpi**2
              nind= nind + 1
              ind(nind)= i
          end if
      end do

c     Print information of this iteration

      if (iprint.ge.0) then 
          write(*,983) iter,ittype
          write(*,985) nprint,(x(i),i=1,nprint)
          write(*,986) nprint,(g(i),i=1,nprint)
          write(*,987) nprint,
     *    (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint)
          write(*,988) min0(nprint,nind),nind,
     *    (ind(i),i=1,min0(nprint,nind))
          write(*,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n,
     *    spgiter,tniter,fcnt,gcnt,cgcnt

          write(10,983) iter,ittype
          write(10,985) nprint,(x(i),i=1,nprint)
          write(10,986) nprint,(g(i),i=1,nprint)
          write(10,987) nprint,
     *    (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint)
          write(10,988) min0(nprint,nind),nind,
     *    (ind(i),i=1,min0(nprint,nind))
          write(10,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n, 
     *    spgiter,tniter,fcnt,gcnt,cgcnt
      end if

c     =======================================================
c     Test some stopping criteria that may occur inside the
c     line searches 
c     =======================================================

      if (lsflag.eq.4) then

          flag= 4

          if (iprint.ge.0) then
              write(*,994)  flag,fmin
              write(10,994) flag,fmin
          end if

          go to 500

      else if (lsflag.eq.6) then

          flag= 6

          if (iprint.ge.0) then
              write(*,996)  flag,mininterp,epsrel,epsabs
              write(10,996) flag,mininterp,epsrel,epsabs
          end if

          go to 500

      else if (lsflag.eq.8) then

          flag= 8

          if (iprint.ge.0) then
              write(*,998)  flag,maxfc
              write(10,998) flag,maxfc
          end if

          go to 500

      end if

c     =======================================================
c     Iterate 
c     =======================================================

      go to 100

c     =======================================================
c     End of main loop
c     =======================================================

c     =======================================================
c     Report output status and return
c     =======================================================

 500  continue

c     Print final information

      if (iprint.ge.0) then

          write(*,982) iter
          write(*,985) nprint,(x(i),i=1,nprint)
          write(*,986) nprint,(g(i),i=1,nprint)
          write(*,987) nprint,
     *    (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint)
          write(*,988) min0(nprint,nind),nind,
     *    (ind(i),i=1,min0(nprint,nind))
          write(*,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n,
     *    spgiter,tniter,fcnt,gcnt,cgcnt

          write(10,982) iter
          write(10,985) nprint,(x(i),i=1,nprint)
          write(10,986) nprint,(g(i),i=1,nprint)
          write(10,987) nprint,
     *    (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint)
          write(10,988) min0(nprint,nind),nind,
     *    (ind(i),i=1,min0(nprint,nind))
          write(10,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n,
     *    spgiter,tniter,fcnt,gcnt,cgcnt

      end if
  
      return 

c     Non-executable statements

 977  format(/1X, 'Entry to GENCAN. Number of variables: ',I7)
 978  format(/1X,'Lower bounds (first ',I6, ' components): ',
     */,6(1X,1PD11.4))
 979  format(/1X,'Upper bounds (first ',I6, ' components): ',
     */,6(1X,1PD11.4))
 980  format(/1X,'Initial point (first ',I6, ' components): ',
     */,6(1X,1PD11.4))
 981  format(/1X,'GENCAN iteration: ',I6, ' (Initial point)')
 982  format(/1X,'GENCAN iteration: ',I6, ' (Final point)')
 983  format(/,1X,'GENCAN iteration: ',I6,
     *' (This point was obtained using a ',A3,' iteration)')
 985  format(1X,'Current point (first ',I6, ' components): ',
     */,6(1X,1PD11.4))
 986  format(1X,'Current gradient (first ',I6, ' components): ',
     */,6(1X,1PD11.4))
 987  format(1X,'Current continuous projected gradient (first ',I6, 
     *' components): ',/,6(1X,1PD11.4))
 988  format(1X,'Current free variables (first ',I6,
     *', total number ',I6,'): ',/,10(1X,I6))
 990  format(/1X,'Flag of GENCAN= ',I2,
     *' (convergence with Euclidian-norm of the projected gradient',
     */,1X,'smaller than ',1PD11.4,')')
 991  format(/1X,'Flag of GENCAN= ',I2,
     *' (convergence with sup-norm of the projected gradient',
     */,1X,'smaller than ',1PD11.4,')')
 992  format(/1X,'Flag of GENCAN= ',I2,
     *' (The algorithm stopped by lack of enough progress. This means',
     */,1X,'that  f(x_k) - f(x_{k+1}) .le. ',1PD11.4,
     *' * max [ f(x_j)-f(x_{j+1}, j < k ]',/,1X,'during ',I7,
     *' consecutive iterations')
 993  format(/1X,'Flag of GENCAN= ',I2,
     *' (The algorithm stopped because the order of the',
     */,1X,'Euclidian-norm of the continuous projected gradient did',
     *' not change during ',/,1X,I7,' consecutive iterations.',
     *' Probably, an exagerately small norm of the',/,1X,'continuous',
     *' projected gradient is required for declaring convergence')
 994  format(/1X,'Flag of GENCAN= ',I2,
     *' (The algorithm stopped because the functional value is',
     */,1X,'smaller than ',1PD11.4)
 996  format(/1X,'Flag of GENCAN= ',I2,
     *' (Too small step in a line search. After having made at least ',
     */,1X,I7,' interpolations, the steplength becames small. Small',
     *' means that we were',/,1X,'at point x with direction d and took',
     *' a step  alpha such that',/,1X,'alpha * |d_i| .lt. max [',
     *1PD11.4,' * |x_i|,',1PD11.4,' ] for all i')
 997  format(/1X,'Flag of GENCAN= ',I2,
     *' (It was exceeded the maximum allowed number of iterations',
     */,1X,'(maxit=',I7,')')
 998  format(/1X,'Flag of GENCAN= ',I2,
     *' (It was exceeded the maximum allowed number of functional',
     */,1X,'evaluations (maxfc=',I7,')')
 999  format(/1X,'Flag of GENCAN= ',I2,
     *' (Error in evalf subroutine)')
 1000 format(/1X,'Flag of GENCAN= ',I2,
     *' (Error in evalg subroutine)')
 1001 format(/1X,'Flag of GENCAN= ',I2,
     *' (Error in evalhd subroutine)')
 1002 format(1X,'Functional value: ', 1PD11.4,
     */,1X,'Euclidian-norm of the continuous projected gradient: ',
     *1PD11.4,
     */,1X,'Euclidian-norm of the internal projection of gp: ',1PD11.4,
     */,1X,'Sup-norm of the continuous projected gradient: ',1PD11.4,
     */,1X,'Free variables at this point: ',I7,
     *' (over a total of ',I7,')',
     */,1X,'SPG iterations: ',I7,
     */,1X,'TN iterations: ',I7,
     */,1X,'Functional evaluations: ',I7,
     */,1X,'Gradient evaluations: ',I7,
     */,1X,'Conjugate gradient iterations: ',I7)

      end


c     ******************************************************
c     ******************************************************
      subroutine spgls(n,x,f,g,l,u,lambda,fmin,maxfc,iprint,fcnt,flag,
     *xtrial,d,gamma,sigma1,sigma2,nint,mininterp,sterel,steabs,epsrel,
     *epsabs,infty) 
 
      integer n,maxfc,iprint,fcnt,flag,mininterp
      double precision x(n),f,g(n),l(n),u(n),lambda,fmin,xtrial(n),d(n),
     *gamma,sigma1,sigma2,nint,sterel,steabs,epsrel,epsabs,infty

c     Safeguarded quadratic interpolation, used in SPG.
c
c     On Entry
c
c     n    integer
c          the order of the x
c
c     x    double precision x(n)
c          current point
c
c     f    double precision
c          function value at the current point
c
c     g    double precision g(n)
c          gradient vector at the current point
c
c     l    double precision l(n)
c          lower bounds
c
c     u    double precision u(n)
c          upper bounds
c
c     lambda double precision
c          spectral steplength
c
c     fmin double precision
c          functional value for the stopping criteria f <= fmin
c
c     maxfc integer
c          maximum number of funtion evaluations
c
c     output logical
c          TRUE: print some information at each iteration,
c          FALSE: no print.
c
c     xtrial, d double precision xtrial(n), d(n)
c          working vectors
c
c     gamma double precision
c          constant for the Armijo crtierion
c          f(x + alpha d) <= f(x) + gamma * alpha * <\nabla f(x),d>
c
c          RECOMMENDED: gamma = 10^{-4}
c
c     sigma1, sigma2 double precision
c          constant for the safeguarded interpolation
c          if alpha_new \notin [sigma1, sigma*alpha] then we take
c          alpha_new = alpha / nint
c
c          RECOMMENDED: sigma1 = 0.1 and sigma2 = 0.9
c
c     nint double precision
c          constant for the interpolation. See the description of
c          sigma1 and sigma2 above. Sometimes we take as a new trial
c          step the previous one divided by nint
c
c          RECOMMENDED: nint = 2.0
c
c     mininterp integer
c          constant for testing if, after having made at least mininterp
c          interpolations, the steplength is soo small. In that case 
c          failure of the line search is declared (may be the direction
c          is not a descent direction due to an error in the gradient 
c          calculations) 
c
c          RECOMMENDED: mininterp = 4
c
c     sterel, steabs double precision
c          this constants mean a ``relative small number'' and ``an 
c          absolute small number'' for the increments in finite
c          difference approximations of derivatives
c
c          RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} 
c
c     epsrel, epsabs, infty  double precision
c          this constants mean a ``relative small number'', ``an 
c          absolute small number'', and ``infinite or a very big
c          number''. Basically, a quantity A is considered negligeble 
c          with respect to another quantity B if 
c          |A| < max ( epsrel * |B|, epsabs ) 
c
c          RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20} and 
c          infty = 10^{+20}
c
c     On Return
c
c     x    double precision
c          final estimation of the solution
c
c     f    double precision
c          functional value at the final estimation 
c
c     fcnt integer
c          number of function evaluations used in the line search   
c
c     flag integer
c          This output parameter tells what happened in this 
c          subroutine, according to the following conventions:
c
c          0= convergence with an Armijo-like criterion
c             (f(xnew) <= f(x) + 1.0d-4 * alpha * <g,d>);
c
c          4= the algorithm stopped because the functional value
c             is very small (f <= fmin);
c
c          6= too small step in the line search. After having made at 
c             least mininterp interpolations, the steplength becames 
c             small. ``small steplength'' means that we are at point
c             x with direction d and step alpha, and, for all i, 
c
c             |alpha * d(i)| .le. max ( epsrel * |x(i)|, epsabs ). 
c 
c             In that case failure of the line search is declared 
c             (maybe the direction is not a descent direction 
c             due to an error in the gradient calculations). Use
c             mininterp > maxfc for inhibit this criterion;
c
c          8= it was achieved the maximum allowed number of
c             function evaluations (maxfc);
c
c          9= error in evalf subroutine.

c     We use the Armijo parameter gamma=1.0d-4. You can replace it by 
c     some other number in (0,1) if you find it convenient.

      logical samep
      integer i,interp,inform
      double precision ftrial,gtd,alpha,atemp

c     Initialization

      interp= 0

c     Compute first trial point, spectral projected gradient 
c     direction, and directional derivative <g,d>.

      alpha= 1.0d0

      gtd= 0.0d0
      do i= 1, n
          xtrial(i)= dmin1(u(i),dmax1(l(i),x(i)-lambda*g(i)))
          d(i)= xtrial(i)-x(i)
          gtd= gtd + g(i)*d(i)
      end do

      call evalq(n,xtrial,ftrial,inform)
      fcnt= fcnt + 1

      if (inform.ne.0) then

          flag = 9

          if (iprint.ge.1) then
              write(*,1000) flag
              write(10,1000) flag
          end if

          go to 500

      end if

c     Print initial information

      if (iprint.ge.1) then
          write(*,980) lambda
          write(*,999) alpha,ftrial,fcnt

          write(10,980) lambda
          write(10,999) alpha, ftrial, fcnt
      end if

c     Main loop

 100  continue

c     Test Armijo stopping criterion

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

          f= ftrial

          do i= 1, n
              x(i)= xtrial(i)
          end do

          flag= 0

          if (iprint.ge.1) then
              write(*,990) flag
              write(10,990) flag
          end if

          go to 500

      end if

c     Test whether f is very small

      if (ftrial.le.fmin) then

          f= ftrial

          do i= 1, n
              x(i)= xtrial(i)
          end do

          flag= 4

          if (iprint.ge.1) then
              write(*,994) flag
              write(10,994) flag
          end if

          go to 500

      end if

c     Test whether the number of functional evaluations is exhausted

      if (fcnt.ge.maxfc) then

          if (ftrial.lt.f) then

              f= ftrial

              do i= 1, n
                  x(i)= xtrial(i)
              end do

          end if

          flag= 8

          if (iprint.ge.1) then
              write(*,998) flag
              write(10,998) flag
          end if

          go to 500

      end if

c     Compute new step (safeguarded quadratic interpolation)

      interp= interp + 1

      if (alpha.lt.sigma1) then
          alpha= alpha/nint      

      else
          atemp= (-gtd*alpha**2)/(2.0d0*(ftrial-f-alpha*gtd))

          if (atemp.lt.sigma1 .or. atemp.gt.sigma2*alpha) then
              alpha= alpha/nint

          else
              alpha= atemp
          end if
      end if

c     Compute new trial point

      do i= 1, n
          xtrial(i)= x(i) + alpha*d(i)
      end do

      call evalq(n,xtrial,ftrial,inform)
      fcnt= fcnt + 1

      if (inform.ne.0) then

          flag = 9

          if (iprint.ge.1) then
              write(*,1000) flag
              write(10,1000) flag
          end if

          go to 500

      end if

c     Print information of this iteration

      if (iprint.ge.1) then
          write(*,999)  alpha,ftrial,fcnt
          write(10,999) alpha,ftrial,fcnt
      end if

c     Test whether at least mininterp interpolations were made and 
c     two consecutive iterates are close enough

      samep= .true.
      do i= 1, n
         if (dabs(alpha*d(i)).gt.dmax1(epsrel*dabs(x(i)),epsabs)) then
             samep= .false.
         end if
      end do

      if (interp.gt.mininterp.and.samep) then

          if (ftrial.lt.f) then

              f= ftrial

              do i= 1, n
                  x(i)= xtrial(i)
              end do

          end if

          flag= 6

          if (iprint.ge.1) then
              write(*,996)  flag
              write(10,996) flag
          end if
  
          go to 500

      end if

c     Iterate

      go to 100

c     Return

 500  continue

      return

c     Non-executable statements

 980  format(/,6x,'SPG (spectral steplength ',1PD11.4,')',
     */,/,6x,'SPG Line search')
 999  format(6x,'Alpha= ',1PD11.4,' F= ',1PD11.4,' FE= ',I5)
 990  format(6x,'Flag of SPG Line search= ',I2,
     *' (Convergence with an Armijo-like criterion)')
 994  format(6x,'Flag of SPG Line search= ',I2,
     *' (Small functional value, smaller than ',/,6X,'parameter fmin)')
 996  format(6x,'Flag of SPG Line search= ',I2,
     *' (Too small step in the interpolation process)')
 998  format(6x,'Flag of SPG Line search= ',I2,
     *' (Too many functional evaluations)')
 1000 format(6x,'Flag of SPG Line search= ',I2,
     *' (Error in evalf subroutine)')

      end


c     ******************************************************
c     ******************************************************
      subroutine cg(nind,ind,n,x,g,delta,l,u,eps,epsnqmp,maxitnqmp,
     *maxit,nearlyq,gtype,htvtype,trtype,iprint,ncomp,s,iter,rbdtype,
     *rbdind,flag,w,y,r,d,sprev,theta,sterel,steabs,epsrel,epsabs,infty)

      logical nearlyq
      integer nind,ind(nind),n,maxitnqmp,maxit,gtype,htvtype,trtype,
     *iprint,ncomp,iter,rbdtype,rbdind,flag
      double precision x(n),g(n),delta,l(n),u(n),eps,epsnqmp,s(n),w(n),
     *y(n),r(n),d(n),sprev(n),theta,sterel,steabs,epsrel,epsabs,infty

c     This subrotuine implements the conjugate-gradient method for 
c     minimizing the quadratic approximation q(s) of f(x) at x, where
c
c     q(s) = 1/2 s^T H s + g^T s,
c
c        H = \nabla^2 f(x),
c
c        g = \nabla f(x),
c
c     subject to ||s||_2 <= delta and l <= x + s <= u.
c
c     The method returns an approximation s to the solution such 
c     that ||H s + g||_2 <= eps * ||g||_2; or converges to the 
c     boundary of ||s||_2 <= delta and l <= x + s <= u; or finds 
c     a point s and a direction d such that q(s + alpha d) = q(s) 
c     for any alpha, i.e., d^T H d = g^T d = 0.
c
c     On Entry
c
c     nind integer
c          number of free variables (this is thee dimension in 
c          which this subroutine will work)
c
c     ind  integer ind(n)
c          array which contains, in the first nind positions, the
c          identifiers of the free variables
c
c     n    integer
c          dimension of the full space
c
c     x    double precision x(n)
c          point at which f function is being approximated by the
c          quadratic model
c
c          The first nind positions of x contains the free variables 
c          x_ind(1), x_ind(2), ..., x_ind(nind).
c
c     g    double precision g(n)
c          linear coeficient of the quadratic function
c
c          This is \nabla f(x) and it also contains in the first nind
c          positions the components g_ind(1), g_ind(2), ..., g_ind(nind).
c    
c          IMPORTANT: the linear algebra of this subroutine lies in 
c          a space of dimension nind. The value of the full dimension n,
c          the non-free variables (which are at the end of array x) and
c          its gradient components (which are at the and of array g)
c          are, at this moment, being used to approximate the hessian
c          times vector products by incremental quotients.
c
c     delta double precision
c          trust region radius (||s||_2 <= delta) 
c    
c     l    double precision l(n)
c          lower bounds on x + s. It components are ordered in the
c          same way as x and g.
c
c     u    double precision u(n)
c          upper bounds on x + s. It components are ordered in the
c          same way as x, g and l.
c
c     eps  double precision
c          tolerance for the stopping criterion 
c          ||H s + g||_2 < eps * ||g||_2
c
c     epsnqmp double precision
c          see below
c
c     maxitnqmp integer
c          This and the previous one parameter are used for a stopping
c          criterion of the conjugate gradient subalgorithm. If the
c          progress in the quadratic model is less or equal than a 
c          fraction of the best progress ( epsnqmp * bestprog ) during
c          maxitnqmp consecutive iterations then CG is stopped by not
c          enough progress of the quadratic model.
c
c          RECOMMENDED: epsnqmp = 1.0d-4, maxitnqmp = 5
c
c     maxit integer
c          maximum number of iterations allowed
c
c     nearlyq logical
c          if function f is (nearly) quadratic, use the option 
c          nerlyq = TRUE. Otherwise, keep the default option.
c
c          if, in an iteration of CG we find a direction d such
c          that d^T H d <= 0 then we take the following decision:
c
c          (i) if nearlyq = TRUE then take direction d and try to 
c          go to the boundary chosing the best point among the two
c          point at the boundary and the current point. 
c
c          (ii) if nearlyq = FALSE then we stop at the current point.
c
c          RECOMMENDED: nearlyq = FALSE
c
c     gtype integer
c          type of gradient calculation
c          gtype = 0 means user suplied evalg subroutine,
c          gtype = 1 means central diference approximation.
c
c          RECOMMENDED: gtype = 0
c
c          (provided you have the evalg subroutine)
c
c     htvtype integer
c          type of gradient calculation
c          htvtype = 0 means user suplied evalhd subroutine,
c          htvtype = 1 means incremental quotients approximation.
c
c          RECOMMENDED: htvtype = 1
c
c          (you take some risk using this option but, a menos que
c          voce tenha uma boa evalhd subroutine, incremental
c          quotients is a very cheap option)
c
c     trtype integer
c          type of trust-region radius
c          trtype = 0 means 2-norm trust-region
c          trtype = 1 means infinite-norm trust-region
c
c          RECOMMENDED: trtype = 0
c
c     output logical
c          TRUE: print some information at each iteration,
c          FALSE: no print.
c
c     w,y,r,d,sprev double precision w(n),y(n),r(n),d(n),sprev(n)
c          working vectors
c
c     theta double precision
c          constant for the angle condition, i.e., at iteration k
c          we need a direction d_k such that
c          <g_k,d_k> <= -theta ||g||_2 ||d_k||_2,
c          where g_k is \nabla f(x_k)
c
c          RECOMMENDED: theta = 10^{-6}
c
c     sterel, steabs double precision
c          this constants mean a ``relative small number'' and ``an 
c          absolute small number'' for the increments in finite
c          difference approximations of derivatives
c
c          RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} 
c
c     epsrel, epsabs, infty  double precision
c          this constants mean a ``relative small number'', ``an 
c          absolute small number'', and ``infinite or a very big
c          number''. Basically, a quantity A is considered negligeble 
c          with respect to another quantity B if 
c          |A| < max ( epsrel * |B|, epsabs ) 
c
c          RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20} and 
c          infty = 10^{+20}
c
c     On Return
c
c     s    double precision s(n)
c          final estimation of the solution
c
c     iter integer
c          number of conjugate-gradient iterations performed
c
c     flag integer
c          termination parameter:
c
c          0= convergence with ||H s + g||_2 <= eps * ||g||_2;
c
c          1= convergence to the boundary of ||s||_2 <= delta;
c
c          2= convergence to the boundary of l - x <= s <= u - x;
c
c          3= stopping with s= s_k  such that 
c             <g,s_k> <= -theta ||g||_2 ||s_k||_2 and 
c             <g,s_{k+1}> > -theta ||g||_2 ||s_{k+1}||_2;
c
c          4= not enough progress of the quadratic model during
c             maxitnqmp iterations, i.e., during maxitnqmp iterations
c             |q-qprev| .le. max ( epsrel * |q|, epsabs );
c
c          6= very similar consecutive iterates, 
c             for two consecutive iterates x and y, for all i
c             |x(i)-y(i)| .le. max ( epsrel * |x(i)|, epsabs );
c
c          7= stopping with d such that d^T H d = 0 and g^T d = 0;
c
c          8= too many iterations;
c
c         11= error in evalhd subroutine.

      character*5 rbdtypea
      logical goth,samep
      integer i,itnqmp,rbdposatype,rbdnegatype,rbdposaind,rbdnegaind,
     *inform
      double precision rnorm2,rnorm2prev,beta,alpha,amax,amax1,amax2, 
     *dtw,dnorm2,dts,aa,bb,cc,dd,gts,dtr,gnorm2,snorm2,q,snorm2prev,
     *qprev,amax1n,amax2n,amaxn,amax2x,amax2nx,qamax,qamaxn,norm2s,
     *bestprog,currprog

c     =======================================================
c     Initialization
c     =======================================================

      goth= .false.

      gnorm2= norm2s(nind,g)

      iter= 0

      itnqmp= 0
      qprev= infty
      bestprog= 0.0d0

      do i= 1, nind
          s(i)= 0.0d0
          r(i)= g(i)
      end do

      q= 0.0d0

      gts= 0.0d0
      snorm2= 0.0d0
      rnorm2= gnorm2

c     =======================================================
c     Print initial information
c     =======================================================

      if (iprint.ge.1) then
          write(*,980) maxit,eps
          if (trtype.eq.0) then
              write(*,981) delta
          else if (trtype.eq.1) then
              write(*,982) delta
          else
              write(*,983)
          end if
          write(*,984) iter,rnorm2,dsqrt(snorm2),q

          write(10,980) maxit,eps
          if (trtype.eq.0) then
              write(10,981) delta
          else if (trtype.eq.1) then
              write(10,982) delta
          else
              write(10,983)
          end if
          write(10,984) iter,rnorm2,dsqrt(snorm2),q

      end if

c     =======================================================
c     Main loop
c     =======================================================

 100  continue

c     =======================================================
c     Test stopping criteria
c     =======================================================

c     if ||r||_2 = ||H s + g||_2 <= eps * ||g||_2 then stop

      if (rnorm2.le.eps**2*gnorm2.or.
     *(rnorm2.le.1.0d-10.and.iter.ne.0)) then

          flag= 0

          if (iprint.ge.1) then
              write(*,990)  flag
              write(10,990) flag
          end if
  
          go to 500

      end if

c     if the maximum number of iterations was achieved then stop

      if (iter.ge.maxit) then

          flag= 8

          if (iprint.ge.1) then
              write(*,998)  flag
              write(10,998) flag
          end if
  
          go to 500

      end if

c     =======================================================
c     Compute direction
c     =======================================================

      if (iter.eq.0) then

          do i= 1, nind
              d(i)= -r(i)
          end do

          dnorm2= rnorm2
          dtr= -rnorm2

      else

          beta= rnorm2/rnorm2prev

          do i= 1, nind
              d(i)= -r(i) + beta*d(i)
          end do

          dnorm2= rnorm2 - 2.0d0*beta*(dtr+alpha*dtw) + beta**2*dnorm2
          dtr= -rnorm2 + beta*(dtr+alpha*dtw)

      end if

c     Force d to be a descent direction of q(s), i.e.,
c     <\nabla q(s), d> = <H s + g, d> = <r, d> \le 0.

      if (dtr.gt.0.0d0) then

          do i= 1, nind
              d(i)= -d(i)
          end do

          dtr= -dtr

      end if

c     =======================================================
c     Compute d^T H d
c     =======================================================

c     w= A d

      call calchd(nind,ind,x,d,g,n,x,gtype,htvtype,w,y,sterel,steabs,
     *goth,inform)

      if (inform.ne.0) then
          flag = 11

          if (iprint.ge.1) then
              write(*,1011) flag
              write(10,1011) flag
          end if

          go to 500
      end if

c     Compute d^T w and ||w||^2

      dtw= 0.0d0
      do i= 1, nind
          dtw= dtw + d(i)*w(i)
      end do 

c     =======================================================
c     Compute maximum step
c     =======================================================

c     amax1 > 0 and amax1n < 0 are the values of alpha such 
c     that ||s + alpha * d||_2 or ||s + alpha * d||_\infty = delta 

      dts= 0.0d0
      do i= 1, nind
          dts= dts + d(i)*s(i)
      end do

c     2-norm trust-region radius

      if (trtype.eq.0) then

          aa= dnorm2
          bb= 2.0d0*dts
          cc= snorm2-delta**2
          dd= dsqrt(bb**2-4.0d0*aa*cc)

          amax1= (-bb+dd)/(2.0d0*aa)
          amax1n= (-bb-dd)/(2.0d0*aa)

c     infinite-norm trust-region radius

      else if (trtype.eq.1) then

          amax1= infty
          amax1n= -infty

          do i= 1, nind
              if (d(i).gt.0.0d0) then
                  amax1= dmin1(amax1, (delta-s(i))/d(i))
                  amax1n= dmax1(amax1n, (-delta-s(i))/d(i))
              else if (d(i).lt.0.0d0) then
                  amax1= dmin1(amax1, (-delta-s(i))/d(i))
                  amax1n= dmax1(amax1n, (delta-s(i))/d(i))
              end if
          end do

      end if

c     amax2 > 0 and amax2n < 0 are the maximum and the minimum
c     values of alpha such that l - x <= s + alpha * d <= u - x,
c     respectively

      amax2= infty
      amax2n= -infty

      do i= 1, nind
          if (d(i).gt.0.0d0) then
              if (u(i).lt.infty) then
                  amax2x= (u(i)-x(i)-s(i))/d(i)
                  if (amax2x.lt.amax2) then
                      amax2= amax2x
                      rbdposaind= i
                      rbdposatype= 2
                  end if
              end if
              if (l(i).gt.-infty) then    
                  amax2nx= (l(i)-x(i)-s(i))/d(i)
                  if (amax2nx.gt.amax2n) then
                      amax2n= amax2nx
                      rbdnegaind= i
                      rbdnegatype= 1
                  end if
              end if
          else if (d(i).lt.0.0d0) then
              if (l(i).gt.-infty) then    
                  amax2x= (l(i)-x(i)-s(i))/d(i)
                  if (amax2x.lt.amax2) then
                      amax2= amax2x
                      rbdposaind= i
                      rbdposatype= 1
                  end if
              end if
              if (u(i).lt.infty) then
                  amax2nx= (u(i)-x(i)-s(i))/d(i)
                  if (amax2nx.gt.amax2n) then
                      amax2n= amax2nx
                      rbdnegaind= i
                      rbdnegatype= 2
                  end if
              end if
          end if
      end do

c     Compute amax as the minimum among amax1 and amax2, and
c     amaxn as the minimum among amax1n and amax2n. Moreover
c     change amaxn by -amaxn to have amax and amaxn as maximum
c     steps along d direction (and not -d in the case of amaxn)

      amax=   dmin1(amax1,  amax2)
      amaxn=  dmax1(amax1n, amax2n)

c     =======================================================
c     Compute the step (and the quadratic functional value at 
c     the new point)
c     =======================================================

      qprev= q

c     If d^T H d > 0 then take the conjugate gradients step

      if (dtw.gt.0.0d0) then

          alpha= dmin1(amax, rnorm2/dtw)

          q= q + 0.5d0*alpha**2*dtw + alpha*dtr

c     If d^T H d <= 0 and function f is nearly quadratic then 
c     take the point with the minimum functional value (q) among 
c     the current one and the ones which are at the boundary, i.e., 
c     the best one between q(s), q(s + amax*d) and q(s + amaxn*d).

      else

          qamax= q + 0.5d0*amax**2*dtw + amax*dtr

c         If we are at iteration zero then take the maximum positive
c         step in the minus gradient direction

          if (iter.eq.0) then

              alpha= amax
              q= qamax

c         If we are not in the first iteration then if function f is 
c         nearly quadratic and q(s + amax * d) or q(s + amaxn * d) is 
c         smaller than q(s), go to the best point in the boundary

          else if (nearlyq.and.(qamax.lt.q.or.qamaxn.lt.q)) then

              qamaxn= q + 0.5d0*amaxn**2*dtw + amaxn*dtr

              if (qamax.lt.qamaxn) then
                  alpha= amax
                  q= qamax
              else
                  alpha= amaxn
                  q= qamaxn
              end if

c         Else, stop at the current point

          else

              flag= 7

              if (iprint.ge.1) then
                  write(*,997)  flag
                  write(10,997) flag
              end if
  
              go to 500

          end if

      end if

c     =======================================================
c     Compute new s
c     =======================================================

      do i= 1, nind
          sprev(i)= s(i)
          s(i)= s(i) + alpha*d(i)
      end do

      snorm2prev= snorm2
      snorm2= snorm2 + alpha**2*dnorm2 + 2.0d0*alpha*dts

c     =======================================================
c     Compute the residual r = H s + g
c     =======================================================

      rnorm2prev= rnorm2

      do i= 1, nind
          r(i)= r(i) + alpha*w(i)
      end do

      rnorm2= norm2s(nind,r)

c     =======================================================
c     Increment number of iterations
c     =======================================================

      iter= iter + 1

c     =======================================================
c     Print information of this iteration
c     =======================================================

      if (iprint.ge.1) then
          write(*,984)  iter,dsqrt(rnorm2),dsqrt(snorm2),q
          write(10,984) iter,dsqrt(rnorm2),dsqrt(snorm2),q
      end if

c     =======================================================
c     Test other stopping criteria
c     =======================================================

c     Test angle condition

      gts= 0.0d0
      do i= 1, nind
          gts= gts + g(i)*s(i)
      end do

      if (gts.gt.0.0d0 .or. gts**2.lt.theta**2*gnorm2*snorm2) then

          do i= 1, nind
              s(i)= sprev(i)
          end do

          snorm2= snorm2prev

          q= qprev

          flag= 3

          if (iprint.ge.1) then
              write(*,993)  flag
              write(10,993) flag
          end if

          go to 500

      end if

c     If we are in the boundary of the box also stop

      if (alpha.eq.amax2.or.alpha.eq.amax2n) then

          if (alpha.eq.amax2) then
              rbdind= rbdposaind
              rbdtype= rbdposatype
          else ! if (alpha.eq.amax2n) then
              rbdind= rbdnegaind
              rbdtype= rbdnegatype
          end if

          if (rbdtype.eq.1) then
              rbdtypea= 'lower'
          else ! if (rbdtype.eq.2) then
              rbdtypea= 'upper'
          end if
 
          flag= 2

          if (iprint.ge.1) then
              write(*,992)  flag,ind(rbdind),rbdtypea
              write(10,992) flag,ind(rbdind),rbdtypea
          end if
  
          go to 500

      end if

c     If we are in the boundary of the trust region then stop

      if (alpha.eq.amax1.or.alpha.eq.amax1n) then

          flag= 1

          if (iprint.ge.1) then
              write(*,991)  flag
              write(10,991) flag
          end if
  
          go to 500

      end if

c     If two consecutive iterates are much close then stop

      samep= .true.
      do i= 1, nind
         if (dabs(alpha*d(i)).gt.dmax1(epsrel*dabs(s(i)),epsabs)) then
              samep= .false.
          end if
      end do

      if (samep) then

          flag= 6

          if (iprint.ge.1) then
              write(*,996)  flag
              write(10,996) flag
          end if
  
          go to 500

      end if

c     Test whether we performed many iterations without good progress
c     of the quadratic model

c     if (dabs(q-qprev).le.dmax1(epsrel*dabs(qprev),epsabs)) then

c         itnqmp= itnqmp + 1

c         if (itnqmp.ge.maxitnqmp) then

c             flag= 4

c             if (iprint.ge.1) then
c                 write(*,994)  flag,itnqmp
c                 write(10,994) flag,itnqmp
c             end if
  
c             go to 500

c         endif

c     else
c         itnqmp= 0
c     endif

c     Test whether we performed many iterations without good progress
c     of the quadratic model 

      currprog= qprev - q
      bestprog= dmax1(currprog, bestprog)

      if (currprog.le.epsnqmp*bestprog) then

          itnqmp= itnqmp + 1

          if (itnqmp.ge.maxitnqmp) then
              flag= 4

              if (iprint.ge.1) then
                  write(*,994)  flag,itnqmp,epsnqmp,bestprog
                  write(10,994) flag,itnqmp,epsnqmp,bestprog
              end if

              go to 500
          endif

      else
          itnqmp= 0
      endif

c     =======================================================
c     Iterate
c     =======================================================

      go to 100

c     =======================================================
c     End of main loop
c     =======================================================

c     =======================================================
c     Return
c     =======================================================

 500  continue

c     Print final information

      if (iprint.ge.1) then
          write(*,985) min0(nind,ncomp),(s(i),i=1,min0(nind,ncomp))
          write(10,985) min0(nind,ncomp),(s(i),i=1,min0(nind,ncomp))
      end if

      return

c     Non-executable statements

 980  format(/,6x,'Conjugate gradients (maxit= ',I7,' acc= ',1PD11.4,
     *')')
 981  format(6x,'Using Euclidian trust region (delta= ',1PD11.4,
     *')')
 982  format(6x,'Using sup-norm trust region (delta= ',1PD11.4,')')
 983  format(6x,'Unknown trust-region type')
 984  format(6x,'CG iter= ',I5,' rnorm: ',1PD11.4,' snorm= ',1PD11.4,
     *' q= ',1PD11.4)
 985  format(/,6x,'Truncated Newton direction (first ',I6, 
     *' components): ',/,1(6x,6(1PD11.4,1x)))
 990  format(6x,'Flag of CG= ',I2,' (Convergence with small residual)')
 991  format(6x,'Flag of CG= ',I2,
     *' (Convergence to the trust region boundary)')
 992  format(6x,'Flag of CG= ',I2,
     *' (Convergence to the boundary of the box constraints,',/,6x,
     *'taking step >= 1, variable ',I7,' will reaches its ',A5,
     *' bound)')
 993  format(6x,'Flag of CG= ',I2,
     *' (The next CG iterate will not satisfy the angle condition)')
 994  format(6x,'Flag of CG= ',I2,
     *' (Not enough progress in the quadratic model. This means',/,6x,
     *'that the progress of the last ',I7,' iterations was smaller ', 
     *'than ',/,6x,1PD11.4,' times the best progress (',1PD11.4,')')
 996  format(6x,'Flag of CG= ',I2,
     *' (Very near consecutive iterates)')
 997  format(6x,'Flag of CG= ',I2,
     *' (d such that d^T H d = 0 and g^T d = 0 was found)')
 998  format(6x,'Flag of CG= ',I2,' (Too many GC iterations)')
 1011 format(6x,'Flag of CG= ',I2,' (Error in evalhd subroutine)')

      end


c     ******************************************************
c     ******************************************************
      subroutine tnls(nind,ind,n,x,f,g,d,l,u,amax,rbdtype,rbdind,fmin,
     *maxfc,gtype,iprint,fcnt,gcnt,intcnt,exgcnt,exbcnt,flag,xplus,
     *xtemp,gamma,beta,sigma1,sigma2,nint,next,mininterp,sterel,steabs,
     *epsrel,epsabs,infty)

      integer nind,ind(nind),n,rbdtype,rbdind,maxfc,gtype,iprint,fcnt,
     *gcnt,intcnt,exgcnt,exbcnt,flag,mininterp
      double precision x(n),f,g(n),d(n),l(n),u(n),amax,fmin,xplus(n),
     *xtemp(n),gamma,beta,sigma1,sigma2,nint,next,sterel,steabs,epsrel,
     *epsabs,infty

c     This subrotuine implements the line search described in a working 
c     paper by E. G. Birgin and J. M. Martinez.
c
c     On Entry
c
c     nind integer
c          number of free variables (this is thee dimension in 
c          which this subroutine will work)
c
c     ind  integer ind(n)
c          array which contains, in the first nind positions, the
c          identifiers of the free variables
c
c     n    integer
c          dimension of the full space
c
c     x    double precision x(n)
c          current point
c
c          The first nind positions of x contains the free variables 
c          x_ind(1), x_ind(2), ..., x_ind(nind).
c
c     f    double precision
c          functional value at x
c
c     g    double precision g(n)
c          gradient vector at x
c
c          It also contains in the first nind positions the components 
c          g_ind(1), g_ind(2), ..., g_ind(nind).
c    
c          IMPORTANT: the linear algebra of this subroutine lies in 
c          a space of dimension nind. The value of the full dimension n,
c          the non-free variables (which are at the end of array x) and
c          its gradient components (which are at the and of array g)
c          are also used and actualized any time the gradient is being
c          evaluated.
c
c     d    double precision d(nind)
c          descent direction 
c    
c     l    double precision l(nind)
c          lower bounds on x. It components are ordered in the
c          same way as x and g.
c
c     u    double precision u(nind)
c          upper bounds on x. It components are ordered in the
c          same way as x, g and l.
c
c     fmin double precision
c          functional value for the stopping criteria f <= fmin
c
c     maxfc integer
c          maximum number of funtion evaluations
c
c     gtype integer
c          type of gradient calculation
c          gtype = 0 means user suplied evalg subroutine,
c          gtype = 1 means central diference approximation.
c
c          RECOMMENDED: gtype = 0
c
c          (provided you have the evalg subroutine)
c
c     output logical
c          TRUE: print some information at each iteration,
c          FALSE: no print.
c
c     xplus, xtemp double precision xplus(nind),xtemp(nind)
c          working vectors
c
c     gamma double precision
c          constant for the Armijo crtierion
c          f(x + alpha d) <= f(x) + gamma * alpha * <\nabla f(x),d>
c
c          RECOMMENDED: gamma = 10^{-4}
c
c     beta double precision
c          constant for the beta condition
c          <d_k, g(x_k + d_k)>  <  beta * <d_k,g_k>
c          if (x_k + d_k) satisfies the Armijo condition but does not
c          satisfy the beta condition then the point is accepted, but
c          if it satisfied the Armijo condition and also satisfies the
c          beta condition then we know that there is the possibility
c          for a succesful extrapolation
c
c          RECOMMENDED: beta = 0.5
c
c     sigma1, sigma2 double precision
c          constant for the safeguarded interpolation
c          if alpha_new \notin [sigma1, sigma*alpha] then we take
c          alpha_new = alpha / nint
c
c          RECOMMENDED: sigma1 = 0.1 and sigma2 = 0.9
c
c     nint double precision
c          constant for the interpolation. See the description of
c          sigma1 and sigma2 above. Sometimes we take as a new trial
c          step the previous one divided by nint
c
c          RECOMMENDED: nint = 2.0
c
c     next double precision
c          constant for the extrapolation
c          when extrapolating we try alpha_new = alpha * next
c
c          RECOMMENDED: next = 2.0
c
c     mininterp integer
c          constant for testing if, after having made at least mininterp
c          interpolations, the steplength is soo small. In that case 
c          failure of the line search is declared (may be the direction
c          is not a descent direction due to an error in the gradient 
c          calculations) 
c
c          RECOMMENDED: mininterp = 4
c
c     sterel, steabs double precision
c          this constants mean a ``relative small number'' and ``an 
c          absolute small number'' for the increments in finite
c          difference approximations of derivatives
c
c          RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} 
c
c     epsrel, epsabs, infty  double precision
c          this constants mean a ``relative small number'', ``an 
c          absolute small number'', and ``infinite or a very big
c          number''. Basically, a quantity A is considered negligeble 
c          with respect to another quantity B if 
c          |A| < max ( epsrel * |B|, epsabs ) 
c
c          RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20} and 
c          infty = 10^{+20}
c
c     On Return
c
c     x    double precision x(n)
c          new current point
c
c     f    double precision
c          functional value at x
c
c     g    double precision g(n)
c          gradient vector at x
c
c     fcnt integer
c          number of funtional evaluations used in this line search
c
c     gcnt integer
c          number of gradient evaluations used in this line search
c
c     intcnt integer
c          number of interpolations
c
c     exgcnt integer
c          number of good extrapolations
c
c     exbcnt integer
c          number of bad extrapolations
c
c     flag integer
c          This output parameter tells what happened in this 
c          subroutine, according to the following conventions:
c
c          0= convergence with an Armijo-like criterion
c             (f(xnew) <= f(x) + 1.0d-4 * alpha * <g,d>);
c
c          4= the algorithm stopped because the functional value
c             is very small (f <= fmin);
c
c          6= soo small step in the line search. After having made at 
c             least mininterp interpolations, the steplength becames 
c             small. ``small steplength'' means that we are at point
c             x with direction d and step alpha, and, for all i, 
c
c             |alpha * d(i)| .le. max ( epsrel * |x(i)|, epsabs ). 
c 
c             In that case failure of the line search is declared 
c             (may be the direction is not a descent direction 
c             due to an error in the gradient calculations). Use
c             mininterp > maxfc for inhibit this criterion;
c
c          8= it was achieved the maximum allowed number of
c             function evaluations (maxfc);
c
c          9= error in evalf subroutine;
c
c         10= error in evalg subroutine.

c     Armijo parameter and beta parameter. You can modify them if 
c     necessary

      logical samep
      integer i,interp,inform
      double precision fplus,atemp,ftemp,gtd,gptd,alpha,fbext

c     =======================================================
c     Compute directional derivative
c     =======================================================

      gtd= 0.0d0
      do i= 1, nind
          gtd= gtd + g(i)*d(i)
      end do

c     =======================================================
c     Compute first trial
c     =======================================================

      alpha= dmin1(1.0d0, amax)

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

      if (alpha.eq.amax) then
          if (rbdtype.eq.1) then
              xplus(rbdind)= l(rbdind)
          else ! if (rbdtype.eq.2) then
              xplus(rbdind)= u(rbdind)
          end if
      end if

      call calcf(nind,ind,xplus,n,x,fplus,inform)
      fcnt= fcnt + 1

      if (inform.ne.0) then

          flag = 9

          if (iprint.ge.1) then
              write(*,1000) flag
              write(10,1000) flag
          end if

          go to 500

      end if

c     Print initial information

      if (iprint.ge.1) then
          write(*,980) amax
          write(*,999) alpha,fplus,fcnt

          write(10,980) amax
          write(10,999) alpha,fplus,fcnt
      end if

c     =======================================================
c     Test Armijo and beta-condition and decide for accepting 
c     the trial point, interpolate or extrapolate.
c     =======================================================

      if (amax.gt.1.0d0) then

c         x + d belongs to the interior of the feasible set (amax > 1)

          if (iprint.ge.1) then
              write(*,*)  '     x+d belongs to int of the feasible set'
              write(10,*) '     x+d belongs to int of the feasible set'
          end if

c         Verify Armijo

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

c             Armijo condition holds  
              if (iprint.ge.1) then
                  write(*,*)  '     Armijo condition holds' 
                  write(10,*) '     Armijo condition holds' 
              end if

              call calcg(nind,ind,xplus,n,x,gtype,g,sterel,steabs,
     *        inform)
              gcnt= gcnt + 1

              if (inform.ne.0) then

                  flag = 10

                  if (iprint.ge.1) then
                      write(*,1001) flag
                      write(10,1001) flag
                  end if

                  go to 500

              end if

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

c             Verify directional derivative (beta condition)

              if (gptd.ge.beta*gtd) then

c                 Step = 1 was ok, finish the line search
                  if (iprint.ge.1) then
                      write(*,*)  '     The beta condition is also true'
                      write(*,*)  '     Line search is over'
                      write(10,*) '     The beta condition is also true'
                      write(10,*) '     Line search is over'
                  end if

                  f= fplus

                  do i= 1, nind
                      x(i)= xplus(i)
                  end do

                  flag= 0

                  if (iprint.ge.1) then
                      write(*,990)  flag
                      write(10,990) flag
                  end if

                  go to 500

              else

c                 Extrapolate
                  if (iprint.ge.1) then
                      write(*,*) '     The beta-condition does not hold'
                      write(*,*) '     We will extrapolate'
                      write(10,*)'     The beta-condition does not hold'
                      write(10,*)'     We will extrapolate'
                  end if

c                 f before extrapolation
                  fbext= fplus

                  go to 100

              end if

          else

c             Interpolate
              if (iprint.ge.1) then
                  write(*,*)  '     Armijo does not hold'
                  write(*,*)  '     We will interpolate'
                  write(10,*) '     Armijo does not hold'
                  write(10,*) '     We will interpolate'
              end if

              go to 200

          end if

      else

c         x + d does not belong to the feasible set (amax <= 1)
          if (iprint.ge.1) then
              write(*,*)  '     x+d does not belong to box-interior'
              write(10,*) '     x+d does not belong to box-interior'
          end if

          if (fplus.lt.f) then

c             Extrapolate
              if (iprint.ge.1) then
                  write(*,*)  '     f(x+d) < f(x)'
                  write(*,*)  '     We will extrapolate'
                  write(10,*) '     f(x+d) < f(x)'
                  write(10,*) '     We will extrapolate'
              end if

c             f before extrapolation
              fbext= fplus

              go to 100

          else

c             Interpolate
              if (iprint.ge.1) then
                  write(*,*)  '     f(x+d) >= f(x)'
                  write(*,*)  '     We will interpolate'
                  write(10,*) '     f(x+d) >= f(x)'
                  write(10,*) '     We will interpolate'
              end if

              go to 200

          end if

      end if


c     =======================================================
c     Extrapolation
c     =======================================================
 100  continue

c     Test f going to -inf

 120  if (fplus.le.fmin) then

c         Finish the extrapolation with the current point

          f= fplus

          do i= 1, nind
              x(i)= xplus(i)
          end do

          call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
          gcnt= gcnt + 1

          if (inform.ne.0) then

              flag = 10

              if (iprint.ge.1) then
                  write(*,1001) flag
                  write(10,1001) flag
              end if

              go to 500

          end if

          if (f.lt.fbext) then
              exgcnt= exgcnt + 1
          else
              exbcnt= exbcnt + 1
          end if

          flag= 4

          if (iprint.ge.1) then
              write(*,994) flag
              write(10,994) flag
          end if

          go to 500

      end if

c     Test maximum number of functional evaluations

      if (fcnt.ge.maxfc) then

c         Finish the extrapolation with the current point

          f= fplus

          do i= 1, nind
              x(i)= xplus(i)
          end do

          call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
          gcnt= gcnt + 1

          if (inform.ne.0) then

              flag = 10

              if (iprint.ge.1) then
                  write(*,1001) flag
                  write(10,1001) flag
              end if

              go to 500

          end if

          if (f.lt.fbext) then
              exgcnt= exgcnt + 1
          else
              exbcnt= exbcnt + 1
          end if

          flag= 8

          if (iprint.ge.1) then
              write(*,998) flag
              write(10,998) flag
          end if

          go to 500

      end if

c     Chose new step 

      if (alpha.lt.amax.and.next*alpha.gt.amax) then
          atemp= amax
      else
          atemp= next*alpha
      end if

c     Compute new trial point

      do i= 1, nind
          xtemp(i)= x(i) + atemp*d(i)
      end do

      if (atemp.eq.amax) then
          if (rbdtype.eq.1) then
              xtemp(rbdind)= l(rbdind)
          else ! if (rbdtype.eq.2) then
              xtemp(rbdind)= u(rbdind)
          end if
      end if

c     Project

      if (atemp.gt.amax) then
          do i= 1, nind
              xtemp(i)= dmax1(l(i),dmin1(xtemp(i), u(i)))
          end do
      end if

c     Test if this is not the same point as the previous one.
c     This test is performed only when alpha >= amax.

      if(alpha.ge.amax) then

          samep= .true.
          do i= 1, nind
              if (dabs(xtemp(i)-xplus(i)).gt.
     *        dmax1(epsrel*dabs(xplus(i)),epsabs)) then
                  samep= .false.
              end if
          end do

          if (samep) then

c             Finish the extrapolation with the current point

              f= fplus

              do i= 1, nind
                  x(i)= xplus(i)
              end do

              call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
              gcnt= gcnt + 1

              if (inform.ne.0) then

                  flag = 10

                  if (iprint.ge.1) then
                      write(*,1001) flag
                      write(10,1001) flag
                  end if

                  go to 500

              end if

              if (f.lt.fbext) then
                  exgcnt= exgcnt + 1
              else
                  exbcnt= exbcnt + 1
              end if

              flag= 0

              if (iprint.ge.1) then
                  write(*,990) flag
                  write(10,990) flag
              end if

              go to 500

          end if

      end if

c     Evaluate function

      call calcf(nind,ind,xtemp,n,x,ftemp,inform)
      fcnt= fcnt + 1

      if (inform.ne.0) then

          flag = 9

          if (iprint.ge.1) then
              write(*,1000) flag
              write(10,1000) flag
          end if

          go to 500

      end if

c     Print information of this iteration

      if (iprint.ge.1) then
          write(*,999)  atemp,ftemp,fcnt
          write(10,999) atemp,ftemp,fcnt
      end if

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

      if (ftemp.lt.fplus) then

          alpha= atemp

          fplus= ftemp

          do i= 1, nind
              xplus(i)= xtemp(i)
          end do

          go to 120

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

      else

          f= fplus

          do i= 1, nind
              x(i)= xplus(i)
          end do

          call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
          gcnt= gcnt + 1

          if (inform.ne.0) then

              flag = 10

              if (iprint.ge.1) then
                  write(*,1001) flag
                  write(10,1001) flag
              end if

              go to 500

          end if

          if (f.lt.fbext) then
              exgcnt= exgcnt + 1
          else
              exbcnt= exbcnt + 1
          end if

          flag= 0

          if (iprint.ge.1) then
              write(*,990) flag
              write(10,990) flag
          end if

          go to 500

      end if
c     =======================================================
c     End of extrapolation
c     =======================================================

c     =======================================================
c     Interpolation
c     =======================================================
 200  continue
      intcnt= intcnt + 1

c     Initialization

      interp= 0

c     Test f going to -inf

 210  if (fplus.le.fmin) then

c         Finish the interpolation with the current point

          f= fplus

          do i= 1, nind
              x(i)= xplus(i)
          end do

          call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
          gcnt= gcnt + 1

          if (inform.ne.0) then

              flag = 10

              if (iprint.ge.1) then
                  write(*,1001) flag
                  write(10,1001) flag
              end if

              go to 500

          end if

          flag= 4

          if (iprint.ge.1) then
              write(*,994) flag
              write(10,994) flag
          end if

          go to 500

      end if

c     Test maximum number of functional evaluations

      if (fcnt.ge.maxfc) then

c         As this is an abrupt termination then the current point of 
c         the interpolation may be worst than the intial one

c         If the current point is better than the initial one then
c         finish the interpolation with the current point else
c         discard all we did inside this line search and finish with
c         the initial point

          if (fplus.lt.f) then

              f= fplus

              do i= 1, nind
                  x(i)= xplus(i)
              end do

              call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
              gcnt= gcnt + 1

              if (inform.ne.0) then

                  flag = 10

                  if (iprint.ge.1) then
                      write(*,1001) flag
                      write(10,1001) flag
                  end if

                  go to 500

              end if

          end if

          flag= 8

          if (iprint.ge.1) then
              write(*,998) flag
              write(10,998) flag
          end if

          go to 500

      end if

c     Compute new step

      interp= interp + 1

      if (alpha.lt.sigma1) then
          alpha= alpha/nint      

      else
          atemp= (-gtd*alpha**2)/(2.0d0*(fplus-f-alpha*gtd))

          if (atemp.lt.sigma1 .or. atemp.gt.sigma2*alpha) then
              alpha= alpha/nint

          else
              alpha= atemp
          end if
      end if

c     Compute new trial point

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

      call calcf(nind,ind,xplus,n,x,fplus,inform)
      fcnt= fcnt + 1

      if (inform.ne.0) then

          flag = 9

          if (iprint.ge.1) then
              write(*,1000) flag
              write(10,1000) flag
          end if

          go to 500

      end if

c     Print information of this iteration

      if (iprint.ge.1) then
          write(*,999)  alpha,fplus,fcnt
          write(10,999) alpha,fplus,fcnt
      end if

c     Test Armijo condition

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

c         Finish the line search
      
          f= fplus

          do i= 1, nind
              x(i)= xplus(i)
          end do

          call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
          gcnt= gcnt + 1

          if (inform.ne.0) then

              flag = 10

              if (iprint.ge.1) then
                  write(*,1001) flag
                  write(10,1001) flag
              end if

              go to 500

          end if

          flag= 0

          if (iprint.ge.1) then
              write(*,990) flag
              write(10,990) flag
          end if

          go to 500

      end if

c     Test whether at least mininterp interpolations were made and 
c     two consecutive iterates are much close

      samep= .true.
      do i= 1, nind
         if (dabs(alpha*d(i)).gt.dmax1(epsrel*dabs(x(i)),epsabs)) then
             samep= .false.
         end if
      end do

      if (interp.gt.mininterp.and.samep) then

c         As this is an abrupt termination then the current point of 
c         the interpolation may be worst than the intial one

c         If the current point is better than the initial one then
c         finish the interpolation with the current point else
c         discard all we did inside this line search and finish with
c         the initial point

c         if (fplus.lt.f) then

c             f= fplus

c             do i= 1, nind
c                 x(i)= xplus(i)
c             end do

c             call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
c             gcnt= gcnt + 1

c             if (inform.ne.0) then

c                 flag = 10

c                 if (iprint.ge.1) then
c                     write(*,1001) flag
c                     write(10,1001) flag
c                 end if

c                 go to 500

c             end if

c         end if

c         The previous lines were commented because, as it is been 
c         used, this subroutine must return with the intial point 
c         in case of finding a very small interpolation step

          flag= 6

          if (iprint.ge.1) then
              write(*,996)  flag
              write(10,996) flag
          end if
  
          go to 500

      end if

c     Else, iterate

      go to 210
c     =======================================================
c     End of interpolation
c     =======================================================

 500  continue

c     =======================================================
c     Return
c     =======================================================

      return

c     Non-executable statements

 980  format(/,6x,'TN Line search (alphamax= ',1PD11.4,')')
 999  format(6x,'Alpha= ',1PD11.4,' F= ',1PD11.4,' FE= ',I5)
 990  format(6x,'Flag of TN Line search= ',I2,
     *' (Convergence with an Armijo-like criterion)')
 994  format(6x,'Flag of TN Line search= ',I2,
     *' (Small functional value, smaller than ',/,6X,'parameter fmin)')
 996  format(6x,'Flag of TN Line search= ',I2,
     *' (Too small step in the interpolation process)')
 998  format(6x,'Flag of TN Line search= ',I2,
     *' (Too many functional evaluations)')
 1000 format(6x,'Flag of TN Line search= ',I2,
     *' (Error in evalf subroutine)')
 1001 format(6x,'Flag of TN Line search= ',I2,
     *' (Error in evalg subroutine)')

      end


c     ******************************************************
c     ******************************************************
      subroutine Gshrink(nind,ind,n,x,g,l,u)

      integer n, nind, ind(nind)
      double precision x(n), g(n), l(n), u(n)

c     Shrink all vectors x, g, l and u from the full dimension
c     space (dimension n) to the reduced space (dimension nind).
 
      integer i, indi
      double precision temp

c     Shrink x, g, l and u for the reduced space

      do i= 1, nind
           indi= ind(i)
           if (i.ne.indi) then
               temp= x(indi)
               x(indi)= x(i)
               x(i)= temp

               temp= g(indi)
               g(indi)= g(i)
               g(i)= temp

               temp= l(indi)
               l(indi)= l(i)
               l(i)= temp

               temp= u(indi)
               u(indi)= u(i)
               u(i)= temp
          end if
      end do

      return

      end


c     ******************************************************
c     ******************************************************
      subroutine Gexpand(nind,ind,n,x,g,l,u)

      integer n, nind, ind(nind)

      double precision x(n), g(n), l(n), u(n)

c     Expands vectors x, g, l and u from the reduced space
c     (dimension nind) to the full space (dimension n).

      integer i, indi
      double precision temp

c     Expand x, g, l and u to the full space

      do i= nind, 1, -1
          indi= ind(i)
          if (i.ne.indi) then
              temp= x(indi)
              x(indi)= x(i)
              x(i)= temp

              temp= g(indi)
              g(indi)= g(i)
              g(i)= temp

              temp= l(indi)
              l(indi)= l(i)
              l(i)= temp

              temp= u(indi)
              u(indi)= u(i)
              u(i)= temp
          end if
      end do
     
      return

      end
 

c     ******************************************************
c     ******************************************************
      subroutine calcf(nind,ind,x,n,xc,al,inform)

      integer nind,ind(nind),n,inform
      double precision x(n),xc(n),al

c     This subroutines computes the augmented lagrangean. It is
c     called from the reduced space (dimension nind), expands
c     the point x where the function will be evaluated and call
c     the subroutine evalf to compute the augmented lagrangean
c     function. Finally, returns vector x to the reduced space 
c     (shrinking). 
c
c     About subroutines named calc[something]. The subroutines
c     whos names start with ``calc'' work in (are called from)
c     the reduced space. Their task is (i) expand the arguments 
c     to the full space, (ii) call the corresponding ``eval'' 
c     subroutine (which works in the full space), and (iii) 
c     shrink the parameters again and also shrink an eventualy
c     output of the ``eval'' subroutines. Subroutines of this 
c     type are: calcf, calcg, calchd and calchddiff (the last 
c     one called from calchd) and the corresponding subroutines
c     in the full space are evalf, evalg and evalhd.

      integer i,indi
      double precision temp

c     Complete x

      do i= nind+1, n
          x(i)= xc(i)
      end do

c     Expand x to the full space

      do i= nind, 1, -1
          indi= ind(i)
          if (i.ne.indi) then
              temp= x(indi)
              x(indi)= x(i)
              x(i)= temp
          end if
      end do
      
c     Compute f

      call evalq(n,x,al,inform)

c     Shrink x to the reduced space

      do i= 1, nind
          indi= ind(i)
          if (i.ne.indi) then
              temp= x(indi)
              x(indi)= x(i)
              x(i)= temp
          end if
      end do

      return

      end


c     ******************************************************
c     ******************************************************
      subroutine calcg(nind,ind,x,n,xc,gtype,nal,sterel,steabs,inform)

      integer nind,ind(nind),n,gtype,inform
      double precision x(n),xc(n),nal(n),sterel,steabs

c     This subroutine computes the augmented lagrangean gradient 
c     vector nal(x). The way it is computed depends on parameter gtype:
c
c     gtype = 0: a user defined function, called evalg, is used,
c     gtype = 1: central finite differences are used.
c
c     It is called from the reduced space (dimension nind), expands
c     the point x where the gradient will be evaluated and call
c     the subroutine evalg to compute the augmented lagrangean 
c     gradient vector. Finally, shrinks vectors x and nal to the
c     reduced space. 
c
c     About subroutines named calc[something]. The subroutines
c     whos names start with ``calc'' work in (are called from)
c     the reduced space. Their task is (i) expand the arguments 
c     to the full space, (ii) call the corresponding ``eval'' 
c     subroutine (which works in the full space), and (iii) 
c     shrink the parameters again and also shrink an eventualy
c     output of the ``eval'' subroutines. Subroutines of this 
c     type are: calcf, calcg, calchd and calchddiff (the last 
c     one called from calchd) and the corresponding subroutines
c     in the full space are evalf, evalg and evalhd.

      integer i,indi
      double precision temp

c     Complete x

      do i= nind+1, n
          x(i)= xc(i)
      end do

c     Expand x to the full space

      do i= nind, 1, -1
          indi= ind(i)
          if (i.ne.indi) then
              temp= x(indi)
              x(indi)= x(i)
              x(i)= temp
          end if
      end do
      
c     Compute gradient vector

c     True gradient

      if (gtype.eq.0) then
          call evalnq(n,x,nal,inform)

c     Central finite diferences approximation

      else if (gtype.eq.1) then
          call evalgdiff(n,x,nal,sterel,steabs,inform)

      end if

c     shrink x and g to the reduced space

      do i= 1, nind
          indi= ind(i)
          if (i.ne.indi) then
              temp= x(indi)
              x(indi)= x(i)
              x(i)= temp

              temp= nal(indi)
              nal(indi)= nal(i)
              nal(i)= temp
          end if    
      end do

      return

      end


c     ******************************************************
c     ******************************************************
      subroutine calchd(nind,ind,x,d,g,n,xc,gtype,htvtype,hd,xtemp,
     *sterel,steabs,goth,inform)

      logical goth
      integer nind,ind(nind),n,gtype,htvtype,inform
      double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n),sterel,steabs

c     This subroutine computes the product Hessian times vector d.
c     The way it is computed depends on parameter htvtype:
c
c     mvptype = 0: a user subroutine, called evalhd, for computing 
c     the hessian times vector d product is called,
c
c     mvptype = 1: a incremental quotients aproximation is used.
c
c     If a user subroutine will be used then this subroutine (which
c     is called from the reduced space) expands vectors x and d, 
c     calls the user supplied subroutine to compute the hessian times
c     vector d product, and shrinks vectors x, d and hd. 
c
c     About subroutines named calc[something]. The subroutines
c     whos names start with ``calc'' work in (are called from)
c     the reduced space. Their task is (i) expand the arguments 
c     to the full space, (ii) call the corresponding ``eval'' 
c     subroutine (which works in the full space), and (iii) 
c     shrink the parameters again and also shrink an eventualy
c     output of the ``eval'' subroutines. Subroutines of this 
c     type are: calcf, calcg, calchd and calchddiff (the last 
c     one called from calchd) and the corresponding subroutines
c     in the full space are evalf, evalg and evalhd.

      integer i,indi
      double precision temp

c     =======================================================
c     User defined subroutine to compute the Hessian times
c     vector d product (it works in the full dimension space)
c     =======================================================

      if (htvtype.eq.0) then

c         Complete d with zeroes

          do i= nind+1, n
              d(i)= 0.0d0
          end do

c         Complete x

          do i= nind+1, n
              x(i)= xc(i)
          end do

c         Expand x, d and g to the full space

          do i= nind, 1, -1
              indi= ind(i)
              if (i.ne.indi) then
                  temp= x(indi)
                  x(indi)= x(i)
                  x(i)= temp

                  temp= d(indi)
                  d(indi)= d(i)
                  d(i)= temp

                  temp= g(indi)
                  g(indi)= g(i)
                  g(i)= temp
              end if
          end do

c         Call the user defined subroutine

          call evalhd(nind,ind,x,d,g,n,xc,gtype,htvtype,hd,xtemp,
     *    sterel,steabs,goth,inform)

c         Shrink x, d, g and hd to the reduced space

          do i= 1, nind
              indi= ind(i)
              if (i.ne.indi) then
                  temp= x(indi)
                  x(indi)= x(i)
                  x(i)= temp

                  temp= d(indi)
                  d(indi)= d(i)
                  d(i)= temp

                  temp= g(indi)
                  g(indi)= g(i)
                  g(i)= temp

                  temp= hd(indi)
                  hd(indi)= hd(i)
                  hd(i)= temp
              end if    
          end do
      

c     =======================================================
c     Incremental quotients approximation
c     (it works in the reduced space)
c     =======================================================

      else if (htvtype.eq.1) then
          call calchddiff(nind,ind,x,d,g,n,xc,gtype,hd,xtemp,sterel,
     *    steabs,inform)

      end if

      return

      end


c ****************************************************
c ****************************************************
      subroutine calchddiff(nind,ind,x,d,g,n,xc,gtype,hd,xtemp,sterel,
     *steabs,inform)

      integer n,nind,ind(n),gtype,inform
      double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n),sterel,steabs

c     This subroutine computes the Hessian times vector d product
c     by means of a ``directional finite difference''. The idea is 
c     that, at the current point x, the product H d is the limit of 
c 
c     [ Gradient(x + t d) - Gradient(x) ] / t
c
c     In this implementation we use
c
c     t = max(steabs, sterel ||x||_\infty) / ||d||_\infty
c
c     provided that d not equal 0, of course. 
c
c     So, we evaluate the Gradient at the auxiliary point x + t d 
c     and use the quotient above to approximate H d. 
c
c     About subroutines named calc[something]. The subroutines
c     whos names start with ``calc'' work in (are called from)
c     the reduced space. Their task is (i) expand the arguments 
c     to the full space, (ii) call the corresponding ``eval'' 
c     subroutine (which works in the full space), and (iii) 
c     shrink the parameters again and also shrink an eventualy
c     output of the ``eval'' subroutines. Subroutines of this 
c     type are: calcf, calcg, calchd and calchddiff (the last 
c     one called from calchd) and the corresponding subroutines
c     in the full space are evalf, evalg and evalhd.
c
c     On Entry
c
c     n     integer
c           order of the x
c
c     x     double precision x(n)
c           point for which Hessian(x) times d will be approximated
c
c     d     double precision d(n)
c           vector for which the Hessian times vetor product will
c           be approximated
c
c     g     double precision g(n)
c           gradient at x
c
c     xtemp double precision xtemp(n)
c           working vector
c
c     sterel, steabs double precision
c           this constants mean a ``relative small number'' and 
c           ``an absolute small number''
c
c     On Return
c
c     hd    double precision hd(n)
c           approximation of H d

      integer i
      double precision xinfn,dinfn,step

c     Compute incremental quotients step

      xinfn= 0.0d0
      dinfn= 0.0d0
      do i= 1, nind
          xinfn= dmax1(xinfn, dabs(x(i)))
          dinfn= dmax1(dinfn, dabs(d(i)))
      end do

      step= dmax1(sterel*xinfn,steabs) / dinfn 

c     Set the point in which the gradient will be 
c     evaluated: xtemp = x + step * d

      do i= 1, nind
          xtemp(i)= x(i) + step*d(i)
      end do

c     Evaluate the gradient at (x + step * d)

      call calcg(nind,ind,xtemp,n,xc,gtype,hd,sterel,steabs,inform)

c     Compute incremental quotients

      do i= 1, nind
          hd(i)= (hd(i)-g(i)) / step
      end do 

      return

      end


c ****************************************************
c ****************************************************
      subroutine evalgdiff(n,x,nal,sterel,steabs,inform)

      integer n,inform
      double precision x(n),nal(n),sterel,steabs

c     Computes the augmented lagrangean gradient vector by central 
c     finite differences. This subroutine, which works in the full 
c     space, is prepared to replace the subroutine evalg (to 
c     evaluate the gradient vector) in the case of the lastest have 
c     not being provided by the user.

      integer j
      double precision tmp,step,fplus,fminus

      do j= 1, n
          tmp= x(j)
          step= dmax1(steabs, sterel*dabs(tmp))

          x(j)= tmp + step
          call evalq(n,x,fplus,inform)

          x(j)= tmp - step
          call evalq(n,x,fminus,inform)

          nal(j)= (fplus-fminus) / (2.0d0*step)
          x(j)= tmp
      end do

      return

      end


c ****************************************************
c ****************************************************

      double precision function norm2s(n,x)

      integer n
      double precision x(n)
      double precision hsldnrm2

      norm2s= hsldnrm2(n,x,1) ** 2

      return 

      end

c ****************************************************
c ****************************************************

      DOUBLE PRECISION FUNCTION HSLDNRM2(N,DX,INCX)
      DOUBLE PRECISION ZERO,ONE
      PARAMETER (ZERO=0.0D0,ONE=1.0D0)
      DOUBLE PRECISION CUTLO,CUTHI
      PARAMETER (CUTLO=8.232D-11,CUTHI=1.304D19)
      INTEGER INCX,N
      DOUBLE PRECISION DX(*)
      DOUBLE PRECISION HITEST,SUM,XMAX
      INTEGER I,J,NEXT,NN
      INTRINSIC DABS,DSQRT,FLOAT
      IF (N.GT.0) GO TO 10
      HSLDNRM2 = ZERO
      GO TO 300
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
      NN = N*INCX
      I = 1
   20 GO TO NEXT
   30 IF (DABS(DX(I)).GT.CUTLO) GO TO 85
      ASSIGN 50 TO NEXT
      XMAX = ZERO
   50 IF (DX(I).EQ.ZERO) GO TO 200
      IF (DABS(DX(I)).GT.CUTLO) GO TO 85
      ASSIGN 70 TO NEXT
      GO TO 105
  100 I = J
      ASSIGN 110 TO NEXT
      SUM = (SUM/DX(I))/DX(I)
  105 XMAX = DABS(DX(I))
      GO TO 115
   70 IF (DABS(DX(I)).GT.CUTLO) GO TO 75
  110 IF (DABS(DX(I)).LE.XMAX) GO TO 115
      SUM = ONE + SUM* (XMAX/DX(I))**2
      XMAX = DABS(DX(I))
      GO TO 200
  115 SUM = SUM + (DX(I)/XMAX)**2
      GO TO 200
   75 SUM = (SUM*XMAX)*XMAX
   85 HITEST = CUTHI/DFLOAT(N)
      DO 95 J = I,NN,INCX
        IF (DABS(DX(J)).GE.HITEST) GO TO 100
   95 SUM = SUM + DX(J)**2
      HSLDNRM2 = DSQRT(SUM)
      GO TO 300
  200 CONTINUE
      I = I + INCX
      IF (I.LE.NN) GO TO 20
      HSLDNRM2 = XMAX*DSQRT(SUM)
  300 CONTINUE
      RETURN
      END

c ****************************************************
c ****************************************************

c     Modifications introduced from March 1st to March 21th of 2002
c     in ocassion of the ISPG development:
c
c     1) Comments of some new parameters introduced in the previous
c     modification
c
c     2) As it was, in the first iteration of GENCAN (when kappa
c     takes value equal 1) and for one-dimensional faces, cgmaxit
c     (the maximum number of Conjugate Gradient iterations to compute
c     the internal to the face truncated-Newton direction) was being 0.
c     As it is obviously wrong, we add a max between what was being
c     computed and one to allow at least one CG iteration.
c
c     3) Parameter inform in subroutines evalf, evalg and evalhd 
c     supplied by the user was added
c
c     Modifications introduced from May 31th to November 2nd of 2001
c     in ocassion of the ALGENCAN development:
c
c     Fixed bugs:
c
c     1) The first spectral steplength was not been projected in the
c     [lammin,lammax] interval.
c
c     2) The conjugate gradients accuracy (cgeps) which is linearly
c     dependent of the euclidian-norm of the projected gradient, was
c     also not been projected in the interval [cgepsi,cgepsf].
c
c     3) Conjugate gradients said that it was being used an euclidian
c     norm trust region when it has really being used an infinite norm
c     trust region and viceversa.
c
c     4) Sometimes, the analytic gradient has been used although the
c     user choose the finite differences option.
c
c     Modifications:
c
c     1) To avoid roundoff errors, an explicit detection of at least one
c     variable reaching its bound when a maximum step is being made was
c     added.
c
c     2) The way in which two points were considered very similar in, for
c     example, the interpolations and the extrapolations (which was 
c     dependent of the infinity norm of the points) showed to be very
c     scale dependent. A new version which test the difference coordinate
c     to coordinate was done. In this was the calculus of the current
c     point x and the descent direction sup-norm is not done any
c     more.
c
c     3) The same constants epsrel and epsabs were used as small
c     relative and absolute values for, for example, detecting similar
c     points and for finite diferences. Now, 
c     epsrel and epsabs are used for detecting similar points 
c     (and the recommended values are 10^{-10} and 10^{-20}, respectively)
c     and new constants sterel ans steabs were introduced for finite
c     differences (and the recommended values are 10^{-7} and 10^{-10},
c     respectively).
c
c     4) Two new stopping criteria for CG were added: (i) we stop if
c     two consecutive iterates are too  close; and (ii) we also
c     stop if there is no enough quadratic model progress during
c     maxitnqmp iterations.
c
c     5) The linear relation between the conjugate gradient accuracy
c     and the norm of the projected gradient can be computed using
c     the euclidian- and the sup-norm of the projected gradient (only
c     euclidian-norm version was present in the previous version. The
c     linear relation is such that the CG accuracy is cgepsi when the
c     projected gradient norm value is equal to the value corresponding
c     to the initial guess and the CG accuracy is cgepsf when the
c     projected gradient norm value is cgrelf).
c
c     6) Inside Conjugate Gradients, the euclidian-norm is been
c     computed using an algorithm developed by C.L.LAWSON, 1978 JAN 08.
c     Numerical experiments showed that the performance of GENCAN depends
c     basically on the conjugate gradients performance and stopping
c     criteria and that the conjugate gradients depends on the way the
c     euclidian-norm is been computed. These things deserve further 
c     research.
c
c     7) In the Augmented Lagrangean algorithm ALGENCAN, which uses
c     GENCAN to solve the bounded constrained subproblems, the maximum
c     number of Conjugate Gradients iterations (cgmaxit), which in
c     this version is linearly dependent of the projected gradient norm, 
c     was set to 2 * (# of free variables). As CG is not using restarts
c     we do not know very well what this means. On the other hand,
c     the accuracy (given by cgeps) continues being more strict when 
c     we are near to the solution and less strict when we ar far from 
c     the solution. 
c
c     8) Many things in the output were changed.


c ****************************************************
c ****************************************************
      subroutine evalhd(nind,ind,y,dy,nl,m,yc,gtype,htvtype,hd,temp,
     *sterel,steabs,goth,inform)

C     Computing H d, where H = lambda A A^t

C     PARAMETERS
      integer mmax
      parameter (mmax    =   100 000)
      integer nmax
      parameter (nmax    =     4 000)
      integer anzemax
      parameter (anzemax =  3 * mmax)

C     SCALAR ARGUMENTS
      logical goth
      integer nind,m,gtype,htvtype,inform
      double precision sterel,steabs

C     ARRAY ARGUMENTS
      integer ind(nind)
      double precision y(m),dy(m),nl(m),yc(m),hd(m),temp(m)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS
      double precision Atd(nmax)

C     COMMON SCALARS
      integer anze,n
      double precision lambda

C     COMMON ARRAYS
      integer arow(anzemax),acol(anzemax)
      double precision g(nmax),bbar(mmax),aval(anzemax)

C     COMMON BLOCKS
      common /ldata1/ arow,acol,aval,anze,n
      common /ldata2/ g,bbar,lambda

      inform = 0

C     lambda A A^t d

C     A^t d

      do i = 1,n
          Atd(i) = 0.0d0
      end do

      do i = 1,anze
          Atd(acol(i)) = Atd(acol(i)) + aval(i) * dy(arow(i))
      end do

C     A A^t d

      do i = 1,m
          hd(i) = 0.0d0
      end do

      do i = 1,anze
          hd(arow(i)) = hd(arow(i)) + aval(i) * Atd(acol(i))
      end do

C     lambda A A^t d

      do i = 1,m
          hd(i) = lambda * hd(i)
      end do

      end

      subroutine evalq(m,y,l,inform)

C     Computing -L(y) = 0.5 lambda || (gk + A^t y) ||_2^2 + bbar^t y

C     PARAMETERS
      integer mmax
      parameter (mmax    =   100 000)
      integer nmax
      parameter (nmax    =     4 000)
      integer anzemax
      parameter (anzemax =  3 * mmax)

C     SCALAR ARGUMENTS
      integer m,inform
      double precision l

C     ARRAY ARGUMENTS
      double precision y(m)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS
      double precision d(nmax)

C     COMMON SCALARS
      integer anze,n
      double precision lambda

C     COMMON ARRAYS
      integer arow(anzemax),acol(anzemax)
      double precision g(nmax),bbar(mmax),aval(anzemax)

C     COMMON BLOCKS
      common /ldata1/ arow,acol,aval,anze,n
      common /ldata2/ g,bbar,lambda

      inform = 0

C     d = g + A^t y

      do i = 1,n
          d(i) = g(i)
      end do

      do i = 1,anze
          d(acol(i)) = d(acol(i)) + aval(i) * y(arow(i))
      end do

C     l = 0.5 lambda d^t d + bbar^t y

      l = 0.0d0
      do i = 1,n
          l = l + d(i) ** 2
      end do
      l = 0.5d0 * lambda * l

      do i = 1,m
          l = l + bbar(i) * y(i)
      end do

      return

      end
    
      subroutine evalnq(m,y,nl,inform)

C     Computing - \nabla L(y) = lambda A (gk + A^t y) + bbar

C     PARAMETERS
      integer mmax
      parameter (mmax    =   100 000)
      integer nmax
      parameter (nmax    =     4 000)
      integer anzemax
      parameter (anzemax =  3 * mmax)

C     SCALAR ARGUMENTS
      integer m,inform

C     ARRAY ARGUMENTS
      double precision y(m),nl(m)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS
      double precision d(nmax)

C     COMMON SCALARS
      integer anze,n
      double precision lambda

C     COMMON ARRAYS
      integer arow(anzemax),acol(anzemax)
      double precision g(nmax),bbar(mmax),aval(anzemax)

C     COMMON BLOCKS
      common /ldata1/ arow,acol,aval,anze,n
      common /ldata2/ g,bbar,lambda

      inform = 0

C     d = g + A^t y

      do i = 1,n
          d(i) = g(i)
      end do

      do i = 1,anze
          d(acol(i)) = d(acol(i)) + aval(i) * y(arow(i))
      end do

C     nl = bbar + lambda A d

      do i = 1,m
          nl(i) = bbar(i)
      end do

      do i = 1,anze
          nl(arow(i)) = nl(arow(i)) + lambda * aval(i) * d(acol(i))
      end do

      return

      end
    
