C     Last update: May 14th, 2003.

C     ALGENCAN main subroutine. This subroutine implements the 67 Augmented
C     Lagrangian methods for inequaltity constraints and bound-constrained
C     problems described in:
C
C     E. G. Birgin, R. Castillo and J. M. Martínez, "Numerical comparison 
C     of Augmented Lagrangian algorithms for nonconvex problems", 
C     Computational Optimization and Applications 31, pp. 31-56, 2005.

      subroutine algencan(pname,n,x,l,u,m,c,mu,mubar,rho,zeta,epsfc,
     +maxit,maxfc,epsnqmp,maxitnqmp,al,nal,gpeucn2,gpsupn,constrv,
     +outiter,totiter,totfcnt,totgcnt,totcgit,totspgit,tottnit,
     +inform,nearlyq,epsgpen,epsgpsn,epsnfp,fmin,delta0,cgmaxit,cgscre,
     +cggpnf,cgepsi,cgepsf,maxitnfp,maxitngp,gtype,htvtype,trtype,
     +iprint,w1,w2,cprev,eta,delmin,lammax,lammin,theta,gamma,beta,
     +sigma1,sigma2,nint,next,mininterp,maxextrap,ncomp,sterel,steabs,
     +epsrel,epsabs,infrel,infabs)

C     SCALAR ARGUMENTS
      logical nearlyq
      character * 10 pname
      integer n,m,maxit,maxfc,maxitnqmp,outiter,totiter,totfcnt,totgcnt,
     +        totcgit,totspgit,tottnit,inform,cgmaxit,maxitnfp,maxitngp,
     +        gtype,htvtype,trtype,iprint,mininterp,maxextrap,ncomp,
     +        cgscre
      double precision epsfc,epsnqmp,al,gpeucn2,gpsupn,
     +        constrv,epsgpen,epsgpsn,epsnfp,fmin,delta0,cggpnf,
     +        cgepsi,cgepsf,eta,delmin,lammax,lammin,theta,gamma,beta,
     +        sigma1,sigma2,nint,next,sterel,steabs,epsrel,epsabs,
     +        infrel,infabs

C     ARRAY ARGUMENTS
      integer w2(n)
      double precision x(n),l(n),u(n),c(m),mu(m),mubar(m),zeta(m),
     +        rho(m),nal(n),w1(10*n+maxitngp),cprev(m)

C     Solves the general-constrained minimization problem
C
C                  Minimize f(x)
C
C                  subject to 
C
C                  c(x) <= 0
C                  l <= x <= u
C     
C     using different Augmented Lagrangean methods (as part of a working 
C     paper by E. G. Birgin, R. Castillo and J. M. Martinez).

C     inform integer
C          This output parameter tells what happened in this 
C          subroutine, according to the following conventions:
C 
C          0 = convergence with feasibility and complementarity;
C
C          1 = it was achieved the maximum allowed number of
C              outer iterations (maxoutit);
C
C          2 = it was achieved the maximum allowed total number of
C              inner iterations (maxtotit);
C
C          3 = it was achieved the maximum allowed total number of
C              functional evaluations (maxtotfc);
C
C          4 = the algorithm stopped by ``lack of feasibility progress'',
C              that means that c(x_k) > epsfc and c(x_k) >= c(x_{k-1}) 
C              during maxitncp consecutive iterations;
C
C        -90 = error in evalfc subroutine;
C
C        -91 = error in evalgjc subroutine;
C
C        -92 = error in evalp subroutine;
C
C        -93 = error in evaldpdy subroutine;
C
C        -94 = error in evalhd subroutine.

C     COMMON SCALAR
      integer ptype,ttype
      double precision tconst,umuini,urhoini,urhomul,urhofrac,umumin,
     +        umumax

C     LOCAL SCALARS
      integer i,iter,fcnt,gcnt,cgcnt,spgiter,spgfcnt,maxtotit,maxtotfc,
     +        tniter,tnfcnt,tnstpcnt,tnintcnt,tnexgcnt,tnexbcnt,tnintfe,
     +        tnexgfe,tnexbfe,maxoutit,nprint,mprint,icrit,maxitncp,flag
      double precision constrvprev,complv,f,v,vprev,uw,uwprev,muprev
C     logical interior

C     COMMON BLOCKS
      common /p/ ptype,ttype,tconst
      common /user/umuini,urhoini,urhomul,urhofrac,umumin,umumax

C     CONSTANTS

C     just for printing
      nprint   = min0(n,ncomp)
      mprint   = min0(m,ncomp)

      maxoutit =       12
      maxitncp = maxoutit
      maxtotit =   100000
      maxtotfc   = 100000

C     INTIALIZATION

      icrit    = 0

      outiter  = 0
      totiter  = 0
      totspgit = 0
      tottnit  = 0
      totfcnt  = 0
      totgcnt  = 0
      totcgit  = 0

      if (iprint.ge.1) then
          write(*, 900) ptype,ttype,pname,n,m,epsfc,epsfc,epsgpsn
          write(10,900) ptype,ttype,pname,n,m,epsfc,epsfc,epsgpsn
      end if

C     COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS

      call evalfc(n,x,f,m,c,flag)

      if ( flag .ne. 0 ) then
        ! ERROR IN THE OBJECTIVE FUNCTION OR THE CONSTRAINTS 
        ! EVALUATION, STOP
          inform = -90

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

          go to 500
      end if

C     SET INITIAL LAGRANGEAN MULTIPLIERS APPROXIMATION, INITIAL PENALTY 
C     PARAMETER AND ZETA

      do i = 1,m
          mu(i)    = umuini
          mubar(i) = max( umumin, min( umumax, mu(i) ) )
          rho(i)   = urhoini
          if ( ptype .eq. 9 .and. ttype .ge. 16 .and. ttype .le. 18 ) 
     +    then
              rho(i) = max( rho(i), 2.0d0 * mubar(i) )
          end if 
          zeta(i)  = rho(i) 
      end do

C     COMPUTE MAXIMUM FEASIBILITY VIOLATION AND COMPLEMENTARITY VIOLATION

      constrv = 0.0d0
      complv  = 0.0d0
      do i = 1,m
          constrv = max( constrv, c(i) )
          if ( c(i) .lt. -epsfc .and. mu(i) .gt. epsfc ) then
              complv = max( complv, - mu(i) * c(i) )
          end if
      end do

C     PRINT INITIAL INFORMATION

      if (iprint.ge.1) then
          write(*,910) outiter,nprint,(x(i),i=1,nprint)
          write(*,920) f,constrv,complv,mprint,(mu(i),i=1,mprint)

          write(10,910) outiter,nprint,(x(i),i=1,nprint)
          write(10,920) f,constrv,complv,mprint,(mu(i),i=1,mprint)
      end if

C     MAIN LOOP

 100  continue

C     TEST STOPPING CRITERIA

      if (outiter .ge. maxoutit) then
        ! MAXIMUM NUMBER OF OUTER ITERATIONS EXCEEDED, STOP
          inform = 1

          if (iprint.ge.1) then
              write(*, 961) inform,maxoutit
              write(10,961) inform,maxoutit
          end if

          go to 500
      end if
 
      if (totiter .gt. maxtotit) then
        ! MAXIMUM NUMBER OF INNER ITERATIONS EXCEEDED, STOP
          inform = 2

          if (iprint.ge.1) then
              write(*, 962) inform,maxtotit
              write(10,962) inform,maxtotit
          end if

          go to 500
      end if
 
      if (totfcnt .gt. maxtotfc) then
        ! MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED, STOP
          inform = 3

          if (iprint.ge.1) then
              write(*, 963) inform,maxtotfc
              write(10,963) inform,maxtotfc
          end if

          go to 500
      end if

C     DO AN OUTER ITERATION

      outiter= outiter + 1      

      if (iprint.ge.1) then
          write(*, 940) outiter
          write(10,940) outiter
      end if

C     CALL THE INNER-SOLVER

      if (iprint.ge.1) then
          write(*, 950) epsgpsn
          write(10,950) epsgpsn
      end if

      call gencan(n,x,l,u,epsgpen,epsgpsn,epsnfp,maxitnfp,maxitngp,fmin,
     +maxit,maxfc,delta0,cgmaxit,cgscre,cggpnf,cgepsi,cgepsf,
     +epsnqmp,maxitnqmp,nearlyq,gtype,htvtype,trtype,iprint,al,nal,
     +gpeucn2,gpsupn,iter,fcnt,gcnt,cgcnt,spgiter,spgfcnt,tniter,tnfcnt,
     +tnstpcnt,tnintcnt,tnexgcnt,tnexbcnt,tnintfe,tnexgfe,tnexbfe,
     +inform,w1(1),w1(n+1),w1(2*n+1),w2(1),w1(3*n+1),w1(3*n+maxitngp+1),
     +eta,delmin,lammax,lammin,theta,gamma,beta,sigma1,sigma2,nint,next,
     +mininterp,maxextrap,ncomp,sterel,steabs,epsrel,epsabs,infrel,
     +infabs)

      totiter  = totiter  + iter
      totspgit = totspgit + spgiter
      tottnit  = tottnit  + tniter
      totfcnt  = totfcnt  + fcnt
      totgcnt  = totgcnt  + gcnt
      totcgit  = totcgit  + cgcnt

      if ( inform .lt. 0 ) then

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

          go to 500
      end if

C     COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS

      do i = 1,m
          cprev(i) = c(i)
      end do

      call evalfc(n,x,f,m,c,flag)

      if ( flag .ne. 0 ) then
        ! ERROR IN THE OBJECTIVE FUNCTION OR THE CONSTRAINTS 
        ! EVALUATION, STOP
          inform = -90

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

          go to 500
      end if

C     COMPUTE NEW LAGRANGEAN MULTIPLIERS APPROXIMATION, UPDATE PENALTY 
C     PARAMETERS AND ZETA

      do i = 1,m
          muprev = mu(i)

          call evaldpdy(c(i),rho(i),mubar(i),zeta(i),mu(i),inform)
          mubar(i) = max( umumin, min( umumax, mu(i) ) )

          if ( inform .lt. 0 ) then
              if (iprint.ge.1) then
                  write(*, 1000) inform
                  write(10,1000) inform
              end if
              go to 500
          end if

          v = max( 0.0d0, c(i) )
          vprev = max( 0.0d0, cprev(i) )

          if ( c(i) .lt. -epsfc .and. mu(i) .gt. epsfc ) then
              uw = - mu(i) * c(i)
          else
              uw = 0.0d0
          end if

          if ( cprev(i) .lt. -epsfc .and. muprev .gt. epsfc ) then
              uwprev = - muprev * cprev(i)
          else
              uwprev = 0.0d0
          end if

          if ( v .gt. urhofrac * vprev .or. uw .gt. urhofrac * uwprev ) 
     +    then
              rho(i) = urhomul * rho(i)
          end if

          if ( ptype .eq. 9 .and. ttype .ge. 16 .and. ttype .le. 18 ) 
     +    then
              rho(i) = max( rho(i), 2.0d0 * mubar(i) )
          end if 

          zeta(i) = rho(i)
      end do

C     COMPUTE MAXIMUM FEASIBILITY VIOLATION AND COMPLEMENTARITY VIOLATION

      constrvprev= constrv

      constrv = 0.0d0
      complv  = 0.0d0
      do i = 1,m
          constrv = max( constrv, c(i) )
          if ( c(i) .lt. -epsfc .and. mu(i) .gt. epsfc ) then
              complv = max( complv, - mu(i) * c(i) )
          end if
      end do

C     PRINT ITERATION INFORMATION

      if (iprint.ge.1) then
          write(*, 910) outiter,nprint,(x(i),i=1,nprint)
          write(*, 920) f,constrv,complv,mprint,(mu(i),i=1,mprint)
          write(*, 925) mprint,(rho(i),i=1,mprint)
          write(*, 930) al,gpsupn,totiter,totfcnt,totgcnt,totcgit

          write(10,910) outiter,nprint,(x(i),i=1,nprint)
          write(10,920) f,constrv,complv,mprint,(mu(i),i=1,mprint)
          write(10,925) mprint,(rho(i),i=1,mprint)
          write(10,930) al,gpsupn,totiter,totfcnt,totgcnt,totcgit
      end if

C     TEST FEASIBILITY AND COMPLEMENTARITY

      if ( constrv .le. epsfc .and. complv .le. epsfc ) then
        ! THE POINT IS FEASIBLE AND COMPLEMENTARITY IS SATISFIED

          flag= 0

          if (iprint.ge.1) then
              write(*, 960) flag,epsfc,epsfc,gpsupn,epsgpsn
              write(10,960) flag,epsfc,epsfc,gpsupn,epsgpsn
          end if

          go to 500
      end if

C     TEST IF THE INFEASIBILITY DECREASED SUFFICIENTLY

      if ( constrv .gt. epsfc .and. constrv .ge. constrvprev ) then

          icrit = icrit + 1

          if ( icrit .ge. maxitncp) then
            ! LACK OF FEASIBILITY PROGRESS, STOP
              flag = 4

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

              go to 500
          end if

      else
          icrit = 0
      end if

C     ITERATE

      go to 100

C     STOP

 500  continue

      return

 900  format( /,' ****** Output of the Augmented Lagrangian',
     +          ' defined by P = ',I1,' Theta = ',I2,' ******'/,
     +        /,' Problem: ',A10,/,
     +        /,' Main characteristics of the problem: ',/,
     +        /,' Number of variables                  = ',4X,I7,
     +        /,' Number of inequality constraints     = ',4X,I7,
     +        /,' Required constraints infeasibility   = ',1PD11.4,
     +        /,' Required complementarity violation   = ',1PD11.4,
     +        /,' Required projected gradient sup-norm = ',1PD11.4)

 910  format( /,' OUTER ITERATION ',I11,/,
     +        /,' Current point (first ',I6, ' components): ',
     +        /,6(1X,1PD11.4))
 920  format( /,' Objective function value      = ',1PD11.4,
     +        /,' Constraints infeasibility     = ',1PD11.4,
     +        /,' Complementarity violation     = ',1PD11.4,/,
     +        /,' Current Lagrange multipliers estimates',
     +          ' (first ',I6, ' components): ',/,6(1X,1PD11.4))
 925  format( /,' Current penalty parameters (first ',I6,
     +          ' components): ',/,6(1X,1PD11.4))
 930  format( /,' Augmented lagrangian value    = ',1PD11.4,
     +        /,' Projected gradient sup-norm   = ',1PD11.4,/,
     +        /,' No. of inner iterations       = ',I11,
     +        /,' No. of functional evaluations = ',I11,
     +        /,' No. of gradient evaluations   = ',I11,
     +        /,' No. of CG iterations          = ',I11)

 940  format( /,' Computing Outer Iteration ',I2,' ...')

 950  format( /,' MINIMIZATION OF THE AUGMENTED LAGRANGIAN ',/,
     +        /,' Required gradient  sup-norm',1PD11.4)

 960  format( /,' Flag of ALGENCAN = ',I3,/,
     +        /,' Solution was found !!! ',/,
     +          ' The infeasibility of the constraints is smaller',
     +          ' than ',1PD11.4,' and the',/,' complementarity',
     +          ' violation is smaller than ',1PD11.4,'. The',
     +          ' sup-norm of',/,' the projected continuous gradient',
     +          ' is ',1PD11.4,' (and something smaller',/,
     +          ' than ',1PD11.4,' was requested).')

 961  format( /,' Flag of ALGENCAN = ',I3,/,
     +        /,' It was exceeded the maximum allowed number of outer',
     +          ' iterations',/,' (maxoutit = ',I7,')') 

 962  format( /,' Flag of ALGENCAN = ',I3,/,
     +        /,' It was exceeded the maximum allowed total number of',
     +          ' inner iterations',/,' (maxtotit = ',I7,')') 

 963  format( /,' Flag of ALGENCAN = ',I3,/,
     +        /,' It was exceeded the maximum allowed total number of',
     +          ' functional evaluations',/,' (maxtotfc = ',I7,')') 

 964  format( /,' Flag of ALGENCAN = ',I3,/,
     +        /,' Constraint violations have not decreased',
     +          ' substantially over',I7,' major',/,' iterations.')

 1000 format( /,' Flag of ALGENCAN = ',I3,' Fatal Error',/,
     +        /,' The following codes means: ',/,
     +        /,' -90 : error in evalfc   subroutine', 
     +        /,' -91 : error in evalgjc  subroutine', 
     +        /,' -92 : error in evalp    subroutine', 
     +        /,' -93 : error in evaldpdy subroutine', 
     +        /,' -94 : error in evalhd   subroutine',/) 

      end

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

c     Last update: March 22th, 2002

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,inform,s,y,d,ind,lastgpns,w,eta,delmin,lammax,lammin,
     *theta,gamma,beta,sigma1,sigma2,nint,next,mininterp,maxextrap,
     *ncomp,sterel,steabs,epsrel,epsabs,infrel,infabs)

C     SCALAR ARGUMENTS
      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,inform,mininterp,
     +        maxextrap,ncomp
      double precision epsgpen,epsgpsn,epsnfp,fmin,udelta0,cggpnf,
     +        cgepsi,cgepsf,epsnqmp,f,gpeucn2,gpsupn,eta,delmin,lammax,
     +        lammin,theta,gamma,beta,sigma1,sigma2,nint,next,sterel,
     +        steabs,epsrel,epsabs,infrel,infabs

C     ARRAY ARGUMENTS
      integer ind(n)
      double precision x(n),l(n),u(n),g(n),s(n),y(n),d(n),
     +        lastgpns(0:maxitngp-1),w(5*n)

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", to appear in Computational Optimization and
c     Applications, 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          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     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     maxextrap integer
c          constant to limit the number of extrapolations in the truncated 
c          Newton direction
c
c          RECOMMENDED: mininterp = maxfc 
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, infrel, infabs  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}, 
c          infrel = 10^{+20} and infabs = 10^{+99}
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     inform 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        < 0 = error in evalal or evalnal subroutines.

C     LOCAL SCALARS
      character * 3 ittype
      integer i,itnfp,nind,nprint,cgmaxit,cgiter,fcntprev,
     +        tnintprev,tnexgprev,tnexbprev,rbdtype,rbdind,infotmp
      double precision lambda,gieucn2,gpi,fprev,sts,sty,epsgpen2,delta,
     +        currprog,bestprog,cgeps,xnorm,ometa2,amax,amaxx,acgeps,
     +        bcgeps,kappa,gpeucn20,gpsupn0,gpnmax

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, and
      fprev= infabs
      bestprog= 0.0d0
      itnfp= 0

c     for testing progress in the projected gradient norm.
      do i = 0,maxitngp-1
          lastgpns(i) = infabs
      end do

c     Print problem information

      if(iprint.ge.2) 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 evalal(n,x,f,inform)
      fcnt= fcnt + 1

      if ( inform .lt. 0 ) then

          if (iprint.ge.2) then
              write(*, 1000) inform
              write(10,1000) inform
          end if

          return
      end if

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

      if ( inform .lt. 0 ) then

          if (iprint.ge.2) then
              write(*,1000)  inform
              write(10,1000) inform
          end if

          return
      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.2) 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
          inform= 4

          if (iprint.ge.2) then
              write(*,994)  inform,fmin
              write(10,994) inform,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
          inform= 8

          if (iprint.ge.2) then
              write(*,998)  inform,maxfc
              write(10,998) inform,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
          inform= 0

          if (iprint.ge.2) then
              write(*,990)  inform,epsgpen
              write(10,990) inform,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
          inform= 1

          if (iprint.ge.2) then
              write(*,991)  inform,epsgpsn
              write(10,991) inform,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
              inform= 2

              if (iprint.ge.2) then
                  write(*,992)  inform,epsnfp,maxitnfp
                  write(10,992) inform,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

      gpnmax = 0.0d0
      do i = 0,maxitngp-1
          gpnmax = max( gpnmax, lastgpns(i) )
      end do

      lastgpns(mod( iter,maxitngp )) = gpeucn2

      if ( gpeucn2 .ge. gpnmax ) then

          inform = 3

          if ( iprint .ge. 2 ) then
              write(*, 993) inform,maxitngp
              write(10,993) inform,maxitngp
          end if

          go to 500

      endif

c     Test whether the number of iterations is exhausted

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

          if (iprint.ge.2) then
              write(*,997)  inform,maxit
              write(10,997) inform,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,inform,
     *    w(1),w(n+1),gamma,sigma1,sigma2,nint,mininterp,sterel,steabs,
     *    epsrel,epsabs,infrel,infabs) 

          spgfcnt= spgfcnt + (fcnt-fcntprev) 

          if ( inform .lt. 0 ) then

              if (iprint.ge.2) then
                  write(*, 1000) inform
                  write(10,1000) inform
              end if

              return
          end if

c         Compute the gradient of the new iterate

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

              if (iprint.ge.2) then
                  write(*, 1000) inform
                  write(10,1000) inform
              end if

              return
          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
c                 delta= dmax1(delmin,100.0d0*dmax1(1.0d0,xnorm))
                  delta= dmax1(delmin,10.0d0*dmax1(1.0d0,xnorm))
              else
                  delta= udelta0
              end if
          else
c             delta= dmax1(delmin,1.0d0*dsqrt(sts))
              delta= dmax1(delmin,10.0d0*dsqrt(sts))
          end if

c         Shrink the point, its gradient and the bounds

          call shrink(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,inform,w(1),w(n+1),w(2*n+1),w(3*n+1),w(4*n+1),
     *    theta,sterel,steabs,epsrel,epsabs,infrel,infabs)

	  cgcnt= cgcnt + cgiter

          if ( inform .lt. 0 ) then

              if (iprint.ge.2) then
                  write(*, 1000) inform
                  write(10,1000) inform
              end if

              return

          end if

c         Compute maximum step

          if ( inform .eq. 2 ) then
              amax = 1.0d0
          else
              amax = infabs
              do i = 1,nind
c                 if (d(i).gt.0.0d0.and.u(i).lt.infrel) then
                  if ( d(i) .gt. 0.0d0 ) then
                      amaxx = ( u(i) - x(i) ) / d(i)
                      if ( amaxx .lt. amax ) then
                          amax    = amaxx
                          rbdind  = i
                          rbdtype = 2
                      end if
c                 else if (d(i).lt.0.0d0.and.l(i).gt.-infrel) then
                  else if ( d(i) .lt. 0.0d0 ) 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,
     +    inform,w(1),w(n+1),w(2*n+1),gamma,beta,sigma1,sigma2,nint,
     +    next,mininterp,maxextrap,sterel,steabs,epsrel,epsabs,infrel,
     +    infabs)

          if ( inform .lt. 0 ) then

              if (iprint.ge.2) then
                  write(*, 1000) inform
                  write(10,1000) inform
              end if

              return

          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 expand(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 ( inform .eq. 6 ) then

              if (iprint.ge.2) 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,
     *        inform,w(1),w(n+1),gamma,sigma1,sigma2,nint,mininterp,
     *        sterel,steabs,epsrel,epsabs,infrel,infabs) 

              spgfcnt= spgfcnt + (fcnt-fcntprev)

              if ( inform .lt. 0 ) then

                  if (iprint.ge.2) then
                      write(*, 1000) inform
                      write(10,1000) inform
                  end if

                  return
              end if

c             Compute the gradient in the new iterate

              infotmp = inform

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

                  if (iprint.ge.2) then
                      write(*,1000)  inform
                      write(10,1000) inform
                  end if

                  return
              end if

              inform = infotmp

          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.2) 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 ( inform .eq. 4 ) then

          if (iprint.ge.2) then
              write(*, 994) inform,fmin
              write(10,994) inform,fmin
          end if

          go to 500

      else if ( inform .eq. 6 ) then

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

          go to 500

      else if ( inform .eq. 8 ) then

          if (iprint.ge.2) then
              write(*, 998) inform,maxfc
              write(10,998) inform,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.2) 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 = ',I3,
     *' (convergence with Euclidian-norm of the projected gradient',
     */,1X,'smaller than ',1PD11.4,')')
 991  format(/1X,'Flag of GENCAN = ',I3,
     *' (convergence with sup-norm of the projected gradient',
     */,1X,'smaller than ',1PD11.4,')')
 992  format(/1X,'Flag of GENCAN= ',I3,
     *' (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 = ',I3,
     *' (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 = ',I3,
     *' (The algorithm stopped because the functional value is',
     */,1X,'smaller than ',1PD11.4)
 996  format(/1X,'Flag of GENCAN = ',I3,
     *' (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 = ',I3,
     *' (It was exceeded the maximum allowed number of iterations',
     */,1X,'(maxit=',I7,')')
 998  format(/1X,'Flag of GENCAN = ',I3,
     *' (It was exceeded the maximum allowed number of functional',
     */,1X,'evaluations (maxfc=',I7,')')
 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)

 1000 format(/1X,'Flag of GENCAN = ',I3,' Fatal Error')

      end


c     ******************************************************
c     ******************************************************
      subroutine spgls(n,x,f,g,l,u,lambda,fmin,maxfc,iprint,fcnt,inform,
     *xtrial,d,gamma,sigma1,sigma2,nint,mininterp,sterel,steabs,epsrel,
     *epsabs,infrel,infabs) 

C     SCALAR ARGUMENTS
      integer n,maxfc,iprint,fcnt,inform,mininterp
      double precision f,lambda,fmin,gamma,sigma1,sigma2,nint,sterel,
     +        steabs,epsrel,epsabs,infrel,infabs

C     ARRAY ARGUMENTS
      double precision x(n),g(n),l(n),u(n),xtrial(n),d(n)
 
c     Safeguarded quadratic interpolation, used in Spectral Projected 
c     Gradient directions.
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, infrel, infabs 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}, 
c          infrel = 10^{+20} and infabs = 10^{+99}
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     inform 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        < 0 = error in evalal subroutine.

C     LOCAL SCALARS
      logical samep
      integer i,interp
      double precision ftrial,gtd,alpha,atemp

c     Print presentation information

      if (iprint.ge.3) then
          write(*, 980) lambda
          write(10,980) lambda
      end if

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 evalal(n,xtrial,ftrial,inform)
      fcnt= fcnt + 1

      if ( inform .lt. 0 ) then

          if (iprint.ge.3) then
              write(*, 1000) inform
              write(10,1000) inform
          end if

          return

      end if

c     Print information of the first trial

      if (iprint.ge.3) then
          write(*, 999) alpha,ftrial,fcnt
          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

          inform = 0

          if (iprint.ge.3) then
              write(*, 990) inform
              write(10,990) inform
          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

          inform = 4

          if (iprint.ge.3) then
              write(*, 994) inform
              write(10,994) inform
          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

          inform = 8

          if (iprint.ge.3) then
              write(*, 998) inform
              write(10,998) inform
          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 evalal(n,xtrial,ftrial,inform)
      fcnt= fcnt + 1

      if ( inform .lt. 0 ) then

          if (iprint.ge.3) then
              write(*, 1000) inform
              write(10,1000) inform
          end if

          return

      end if

c     Print information of the current trial

      if (iprint.ge.3) 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

          inform = 6

          if (iprint.ge.3) then
              write(*, 996) inform
              write(10,996) inform
          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 = ',I3,
     *          ' (Convergence with an Armijo-like criterion)')
 994  format(6x,'Flag of SPG Line search = ',I3,
     *          ' (Small functional value, smaller than ',/,
     *       6X,'parameter fmin)')
 996  format(6x,'Flag of SPG Line search = ',I3,
     *          ' (Too small step in the interpolation process)')
 998  format(6x,'Flag of SPG Line search = ',I3,
     *          ' (Too many functional evaluations)')
 1000 format(6x,'Flag of SPG Line search = ',I3,' Fatal Error')

      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,inform,w,y,r,d,sprev,theta,sterel,steabs,epsrel,epsabs,
     *infrel,infabs)

C     SCALAR ARGUMENTS
      logical nearlyq
      integer nind,n,maxitnqmp,maxit,gtype,htvtype,trtype,iprint,ncomp,
     +        iter,rbdtype,rbdind,inform
      double precision delta,eps,epsnqmp,theta,sterel,steabs,epsrel,
     +        epsabs,infrel,infabs

C     ARRAY ARGUMENTS
      integer ind(nind)
      double precision x(n),g(n),l(n),u(n),s(n),w(n),y(n),r(n),d(n),
     +        sprev(n)

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, infrel, infabs 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}, 
c          infrel = 10^{+20} and infabs = 10^{+99}
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     inform 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        < 0 = error in evalhd subroutine.

C     LOCAL SCALARS
      character * 5 rbdtypea
      logical goth,samep
      integer i,itnqmp,rbdposatype,rbdnegatype,rbdposaind,rbdnegaind
      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    = infabs
      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. 3 ) 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

          inform = 0

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

      end if

c     if the maximum number of iterations was achieved then stop

      if ( iter .ge. maxit ) then

          inform = 8

          if ( iprint .ge. 3 ) then
              write(*, 998) inform
              write(10,998) inform
          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 .lt. 0 ) then

          if ( iprint .ge. 3 ) then
              write(*, 1000) inform
              write(10,1000) inform
          end if

          return
      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 = sqrt( 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  =  infabs
          amax1n = -infabs

          do i = 1,nind
              if ( d(i) .gt. 0.0d0 ) then
                  amax1  = min( amax1,  (   delta - s(i) ) / d(i) )
                  amax1n = max( amax1n, ( - delta - s(i) ) / d(i) )
              else if ( d(i) .lt. 0.0d0 ) then
                  amax1  = min( amax1,  ( - delta - s(i) ) / d(i) )
                  amax1n = max( 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  =   infabs
      amax2n = - infabs

      do i = 1,nind
          if ( d(i) .gt. 0.0d0 ) then
c             if (u(i).lt.infrel) then
                  amax2x = ( u(i) - x(i) - s(i) ) / d(i)
                  if ( amax2x .lt. amax2 ) then
                      amax2       = amax2x
                      rbdposaind  = i
                      rbdposatype = 2
                  end if
c             end if
c             if (l(i).gt.-infrel) then    
                  amax2nx = ( l(i) - x(i) - s(i) ) / d(i)
                  if ( amax2nx .gt. amax2n ) then
                      amax2n      = amax2nx
                      rbdnegaind  = i
                      rbdnegatype = 1
                  end if
c             end if
          else if ( d(i) .lt. 0.0d0 ) then
c             if (l(i).gt.-infrel) then    
                  amax2x = ( l(i) - x(i) - s(i) ) / d(i)
                  if ( amax2x .lt. amax2 ) then
                      amax2       = amax2x
                      rbdposaind  = i
                      rbdposatype = 1
                  end if
c             end if
c             if (u(i).lt.infrel) then
                  amax2nx = ( u(i) - x(i) - s(i) ) / d(i)
                  if ( amax2nx .gt. amax2n ) then
                      amax2n      = amax2nx
                      rbdnegaind  = i
                      rbdnegatype = 2
                  end if
c             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  = min( amax1 , amax2  )
      amaxn = max( 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 = min( 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 

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

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

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

c         Else, stop at the current point

              else

                  inform = 7

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

              end if

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

          inform = 3

          if ( iprint .ge. 3 ) then
              write(*, 993) inform
              write(10,993) inform
          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

C         HABIA!

c         ACA HAY UN ERROR SI ALPHA = 1.0D+20 Y AMAX2 TAMBIEN ENTONCES
c         NO ES VERDAD QUE ESTAMOS EN EL BORDE DE LA CAJA Y RBDIND, QUE
c         SERA USADO COMO INDICE DE VETOR EN BREVE, NO ESTA DEFINIDO !!!
c         PARA RESOLVERLO TUVE QUE USAR 1.0D+99 Y -1.0D+99 PARA INICIALIZAR
c         EL CALCULO DE LOS AMAX, EN LUGAR DE 1.0D+20 Y -1.0D+20.
c         OTRA SOLUCION ES QUE EL RADIO DE LA REGION DE CONFIANZA NO SEA
c         NUNCA MAYOR QUE 1.0D+20.

          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
 
          inform = 2

          if ( iprint .ge. 3 ) then
              write(*, 992) inform,ind(rbdind),rbdtypea
              write(10,992) inform,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

          inform = 1

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

      end if

c     If two consecutive iterates are much close then stop

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

      if ( samep ) then

          inform = 6

          if ( iprint .ge. 3 ) then
              write(*, 996) inform
              write(10,996) inform
          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             inform = 4

c             if (iprint.ge.3) then
c                 write(*,994)  inform,itnqmp
c                 write(10,994) inform,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 = max( currprog, bestprog )

      if ( currprog .le. epsnqmp * bestprog ) then

          itnqmp = itnqmp + 1

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

              if ( iprint .ge. 3 ) then
                  write(*, 994) inform,itnqmp,epsnqmp,bestprog
                  write(10,994) inform,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. 3 ) 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 = ',I3,' (Convergence with small residual)')
 991  format(6x,'Flag of CG = ',I3,
     *' (Convergence to the trust region boundary)')
 992  format(6x,'Flag of CG = ',I3,
     *' (Convergence to the boundary of the box constraints,',/,6x,
     *'taking step >= 1, variable ',I6,' will reaches its ',A5,
     *' bound)')
 993  format(6x,'Flag of CG = ',I3,
     *' (The next CG iterate will not satisfy the angle condition)')
 994  format(6x,'Flag of CG = ',I3,
     *' (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 = ',I3,
     *' (Very near consecutive iterates)')
 997  format(6x,'Flag of CG= ',I3,
     *' (d such that d^T H d = 0 and g^T d = 0 was found)')
 998  format(6x,'Flag of CG = ',I3,' (Too many GC iterations)')
 1000 format(6x,'Flag of CG = ',I3,' Fatal Error')

      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,inform,xplus,
     +xtemp,xbext,gamma,beta,sigma1,sigma2,nint,next,mininterp,
     +maxextrap,sterel,steabs,epsrel,epsabs,infrel,infabs)

C     SCALAR ARGUMENTS
      integer nind,n,rbdtype,rbdind,maxfc,gtype,iprint,fcnt,gcnt,intcnt,
     +        exgcnt,exbcnt,inform,mininterp,maxextrap
      double precision f,amax,fmin,gamma,beta,sigma1,sigma2,nint,next,
     +        sterel,steabs,epsrel,epsabs,infrel,infabs

C     ARRAY ARGUMENTS
      integer ind(nind)
      double precision x(n),g(n),d(n),l(n),u(n),xplus(n),xtemp(n),
     +        xbext(n)

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, xbext double precision xplus(nind),xtemp(nind),xbext(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     maxextrap integer
c          constant to limit the number of extrapolations
c
c          RECOMMENDED: maxextrap = 1000 (a big number)
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, infrel, infabs 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}, 
c          infrel = 10^{+20} and infabs = 10^{+99}
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     inform 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        < 0 = error in evalal or evalnal subroutines.

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

C     LOCAL SCALARS
      logical samep
      integer i,interp,extrap
      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 = min( 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 calcal(nind,ind,xplus,n,x,fplus,inform)
      fcnt = fcnt + 1

      if ( inform .lt. 0 ) then
 
          if ( iprint .ge. 3 ) then
              write(*, 1000) inform
              write(10,1000) inform
          end if

          return

      end if

c     Print initial information

      if ( iprint .ge. 3 ) 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
          if ( iprint .ge. 3 ) 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. 3 ) then
                  write(*, *) '     Armijo condition holds' 
                  write(10,*) '     Armijo condition holds' 
              end if

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

              if ( inform .lt. 0 ) then

                  if ( iprint .ge. 3 ) then
                      write(*, 1000) inform
                      write(10,1000) inform
                  end if

                  return

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

                  inform = 0

                  if ( iprint .ge. 3 ) then
                      write(*, 990) inform
                      write(10,990) inform
                  end if

                  go to 500

              else

c                 Extrapolate
                  if ( iprint .ge. 3 ) 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 and x before extrapolation
                  fbext = fplus

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

                  go to 100

              end if

          else

c             Interpolate
              if ( iprint .ge. 3 ) 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. 3 ) 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. 3 ) 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 and x before extrapolation
              fbext = fplus

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

              go to 100

          else

c             Interpolate
              if ( iprint .ge. 3 ) 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

      extrap = 0

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

          if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then

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

              if ( inform .lt. 0 ) then
 
                  if ( iprint .ge. 3 ) then
                      write(*, 1000) inform
                      write(10,1000) inform
                  end if

                  return

              end if

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

          end if

          inform = 4

          if ( iprint .ge.3 ) then
              write(*, 994) inform
              write(10,994) inform
          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

c         If extrap=0 and amax>1 the gradient was computed for testing 
c         the beta condition and it is not necessary to compute it again
          if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then

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

              if ( inform .lt. 0 ) then

                  if ( iprint .ge. 3 ) then
                      write(*, 1000) inform
                      write(10,1000) inform
                  end if

                  return

              end if

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

          end if

          inform = 8

          if ( iprint .ge. 3 ) then
              write(*, 998) inform
              write(10,998) inform
          end if

          go to 500

      end if

c     Test if the maximum number of extrapolations was exceeded

      if ( extrap .ge. maxextrap ) then

c         Finish the extrapolation with the current point

          f = fplus

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

c         If extrap=0 and amax>1 the gradient was computed for testing 
c         the beta condition and it is not necessary to compute it again
          if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then

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

              if ( inform .lt. 0 ) then

                  if ( iprint .ge. 3 ) then
                      write(*, 1000) inform
                      write(10,1000) inform
                  end if

                  return

              end if

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

          end if

          inform = 7

          if ( iprint .ge. 3 ) then
              write(*, 997) inform
              write(10,997) inform
          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) = max( l(i), min( 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 .gt. amax ) then

          samep = .true.
          do i = 1,nind
              if ( abs( xtemp(i) - xplus(i) ) .gt.
     +             max( epsrel * abs( 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

c             If extrap=0 and amax>1 the gradient was computed for testing 
c             the beta condition and it is not necessary to compute it again
              if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then

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

                  if ( inform .lt. 0 ) then

                      if ( iprint .ge. 3 ) then
                          write(*, 1000) inform
                          write(10,1000) inform
                      end if

                      return

                  end if

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

              end if

              inform = 0

              if ( iprint .ge. 3 ) then
                  write(*, 990) inform
                  write(10,990) inform
              end if

              go to 500

          end if

      end if

c     Evaluate function

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

      if ( inform .lt. 0 ) then

c         if ( iprint .ge. 3 ) then
c             write(*, 1000) inform
c             write(10,1000) inform
c         end if

c         return

c         If the objective function is not well defined in an 
c         extrapolated point, we discard all the extrapolated points
c         and return to a safe region (where the point before
c         starting the extrapolations is)

          f = fbext

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

c         If extrap=0 and amax>1 the gradient was computed for testing 
c         the beta condition and it is not necessary to compute it again
          if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then

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

              if ( inform .lt. 0 ) then

                  if ( iprint .ge. 3 ) then
                      write(*, 1000) inform
                      write(10,1000) inform
                  end if

                  return

              end if

              exbcnt = exbcnt + 1

          end if

          inform = 0

          if ( iprint .ge. 3 ) then
              write(*, 1010) inform
              write(10,1010) inform
          end if

          go to 500

      end if

c     Print information of this iteration

      if ( iprint .ge. 3 ) 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

          extrap = extrap + 1

          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

c         If extrap=0 and amax>1 the gradient was computed for testing 
c         the beta condition and it is not necessary to compute it again
          if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then

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

              if ( inform .lt. 0 ) then

                  if ( iprint .ge. 3 ) then
                      write(*, 1000) inform
                      write(10,1000) inform
                  end if

                  return

              end if

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

          end if

          inform = 0

          if ( iprint .ge.3 ) then
              write(*, 990) inform
              write(10,990) inform
          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 calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
          gcnt = gcnt + 1

          if ( inform .lt. 0 ) then

              if ( iprint .ge. 3 ) then
                  write(*, 1000) inform
                  write(10,1000) inform
              end if

              return

          end if

          inform = 4

          if ( iprint .ge. 3 ) then
              write(*, 994) inform
              write(10,994) inform
          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 calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
              gcnt = gcnt + 1

              if ( inform .lt. 0 ) then

                  if ( iprint .ge. 3 ) then
                      write(*, 1000) inform
                      write(10,1000) inform
                  end if

                  return

              end if

          end if

          inform = 8

          if ( iprint .ge. 3 ) then
              write(*, 998) inform
              write(10,998) inform
          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 calcal(nind,ind,xplus,n,x,fplus,inform)
      fcnt = fcnt + 1

      if ( inform .lt. 0 ) then

          if ( iprint .ge. 3 ) then
              write(*, 1000) inform
              write(10,1000) inform
          end if

          return

      end if

c     Print information of this iteration

      if ( iprint .ge. 3 ) 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 calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
          gcnt = gcnt + 1

          if ( inform .lt. 0 ) then

              if ( iprint .ge. 3 ) then
                  write(*, 1000) inform
                  write(10,1000) inform
              end if

              return

          end if

          inform = 0

          if ( iprint .ge. 3 ) then
              write(*, 990) inform
              write(10,990) inform
          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 ( abs( alpha * d(i) ) .gt. 
     +        max( epsrel * abs( 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 calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform)
c             gcnt = gcnt + 1

c             if ( inform .lt. 0 ) then 

c                 if ( iprint .ge. 3 ) then
c                     write(*, 1000) inform
c                     write(10,1000) inform
c                 end if

c                 return

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. From
c         that intial point, something different will be tried.

          inform = 6

          if ( iprint .ge. 3 ) then
              write(*, 996) inform
              write(10,996) inform
          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= ',I3,
     +          ' (Convergence with an Armijo-like criterion)')
 994  format(6X,'Flag of TN Line search= ',I3,
     +          ' (Small functional value, smaller than ',/,
     +       6X,'parameter fmin)')
 996  format(6X,'Flag of TN Line search= ',I3,
     +          ' (Too small step in the interpolation process)')
 997  format(6X,'Flag of TN Line search= ',I3,
     +          ' (Too many extrapolations)')
 998  format(6X,'Flag of TN Line search= ',I3,
     +          ' (Too many functional evaluations)')
 1000 format(6X,'Flag of TN Line search = ',I3,' Fatal Error')
 1010 format(6X,'Flag of TN Line search= ',I3,
     +          ' (Fatal Error in an extrapolated point)')

      end


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

C     SCALAR ARGUMENTS
      integer n,nind

C     ARRAY ARGUMENTS
      integer 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).

C     LOCAL SCALARS
      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 expand(nind,ind,n,x,g,l,u)

C     SCALAR ARGUMENTS
      integer n, nind

C     ARRAY ARGUMENTS
      integer 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).

C     LOCAL SCALARS
      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 calcal(nind,ind,x,n,xc,al,inform)

C     SCALAR ARGUMENTS
      integer nind,n,inform
      double precision al

C     ARRAY ARGUMENTS
      integer ind(nind)
      double precision x(n),xc(n)

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 evalal 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: calcal, calcnal, 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     LOCAL SCALARS
      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 evalal(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 calcnal(nind,ind,x,n,xc,gtype,nal,sterel,steabs,inform)

C     SCALAR ARGUMENTS
      integer nind,n,gtype,inform
      double precision sterel,steabs

C     ARRAY ARGUMENTS
      integer ind(nind)
      double precision x(n),xc(n),nal(n)

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 evalnal, 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 evalnal 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: calcal, calcnal, 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     LOCAL SCALARS
      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 evalnal(n,x,nal,inform)

c     Central finite diferences approximation

      else if (gtype.eq.1) then
          call evalnaldiff(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)

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

C     ARRAY ARGUMENTS
      integer ind(nind)
      double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n)

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: calcal, calcnal, 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     LOCAL SCALARS
      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)

C     SCALAR ARGUMENTS
      integer n,nind,gtype,inform
      double precision sterel,steabs

C     ARRAY ARGUMENTS
      integer ind(n)
      double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n)

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: calcal, calcnal, 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

C     LOCAL SCALARS
      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 calcnal(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 evalal(n,x,al,inform)

C     SCALAR ARGUMENTS
      double precision al
      integer n,inform

C     ARRAY ARGUMENTS
      double precision x(n)

C     This subroutine computes the augmented lagrangian function
C
C     On Entry:
C
C     n     integer,
C           size of the problem,
C
C     x     double precision x(n),
C           point at which the function will be evaluated.
C
C     On Return
C
C     al    double precision,
C           augmented lagrangean value at x,
C
C     inform integer,
C           termination parameter:
C           0= the function was succesfully evaluated,
C           1= some error occurs in the function evaluation.

C     PARAMETERS
      integer mmax
      parameter (mmax=100000)

C     COMMON SCALARS
      integer m
 
C     COMMON ARRAYS
      double precision mubar(mmax),zeta(mmax),rho(mmax)

C     LOCAL SCALARS
      integer i,flag
      double precision f,p

C     LOCAL ARRAYS
      double precision c(mmax)

C     COMMON BLOCKS
      common /auglag/mubar,zeta,rho,m

C     INITIALIZE

      inform = 0

C     COMPUTE OBJECTIVE FUNTION AND CONSTRAINTS

      call evalfc(n,x,f,m,c,flag)
      if ( flag .ne. 0 ) then
          inform = -90
          return
      end if

C     SAVE THE CONSTRAINTS FOR THE GRADIENT EVALUATION
      call savec(n,x,m,c)

C     COMPUTES AL(x) = f(x) + \sum_i=1^m P(c(i),rho(i),mu(i))

      al = f

      do i = 1,m

          call evalp(c(i),rho(i),mubar(i),zeta(i),p,inform)     
          if ( inform .lt. 0 ) then
              return
          end if

          al = al + p

      end do

      end


c ****************************************************
c ****************************************************
      subroutine evalnal(n,x,nal,inform)

C     SCALAR ARGUMENTS
      integer n,inform

C     ARRAY ARGUMENTS
      double precision nal(n),x(n)

C     This subroutine computes the gradient of the 
C     augmented lagrangian function.
C
C     On Entry:
C
C     n     integer,
C           size of the problem,
C
C     x     double precision x(n),
C           point at which the gradient will be evaluated.
C
C     On Return
C
C     nal   double precision nal(n),
C           gradient of the augmented lagrangean at x,
C
C     inform integer,
C           termination parameter:
C           0= the gradient was succesfully evaluated,
C           1= some error occurs in the gradient evaluation.

C     PARAMETERS
      integer mmax,jcnzemax
      parameter (mmax=100000,jcnzemax=100000)

C     COMMON SCALARS
      integer m
 
C     COMMON ARRAYS
      double precision mubar(mmax),zeta(mmax),rho(mmax)

C     LOCAL SCALARS
      logical got
      integer i,jcnze,flag
      double precision f

C     LOCAL ARRAYS
      integer jcvar(jcnzemax),jcfun(jcnzemax)
      double precision c(mmax),dpdc(mmax),jcval(jcnzemax)
 
C     COMMON BLOCKS
      common /auglag/mubar,zeta,rho,m

C     INITIALIZE

      inform = 0

C     COMPUTE CONSTRAINTS

      call gotc(n,x,m,c,got)

      if ( .not. got ) then
          call evalfc(n,x,f,m,c,flag)
          if ( flag .ne. 0 ) then
              inform = -90
              return
          end if
      end if

C     COMPUTE dP/dc

      do i = 1,m

          call evaldpdy(c(i),rho(i),mubar(i),zeta(i),dpdc(i),inform)
          if ( inform .lt. 0 ) then
              return
          end if

      end do

C     COMPUTE \nabla L(x) = \nabla f(x) + \sum_{i=1}^m dPdc * dcdx

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

      call evalgjc(n,x,m,jcvar,jcfun,jcval,jcnze,flag)
      if ( flag .ne. 0 ) then
          inform = -91
          return
      end if

      do i = 1,jcnze
          if ( jcfun(i).eq.0 ) then
              nal(jcvar(i)) = nal(jcvar(i)) + jcval(i)
          else
              nal(jcvar(i)) = 
     +        nal(jcvar(i)) + dpdc(jcfun(i)) * jcval(i)
          end if
      end do

      end


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

C     SCALAR ARGUMENTS
      integer n,inform
      double precision sterel,steabs

C     ARRAY ARGUMENTS
      double precision x(n),nal(n)

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 evalnal (to 
c     evaluate the gradient vector) in the case of the lastest have 
c     not being provided by the user.

C     LOCAL SCALARS
      integer j
      double precision tmp,step,fplus,fminus

      inform = 0

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

          x(j)= tmp + step
          call evalal(n,x,fplus,inform)
          if ( inform .lt. 0 ) then
              return
          end if

          x(j)= tmp - step
          call evalal(n,x,fminus,inform)
          if ( inform .lt. 0 ) then
              return
          end if

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

      return

      end


c ****************************************************
c ****************************************************
      subroutine evalp(y,rho,mu,zeta,p,inform)

C     SCALAR ARGUMENTS
      integer inform
      double precision y,rho,mu,zeta,p

C     COMMON SCALARS
      integer ptype,ttype
      double precision tconst

C     LOCAL SCALARS
      double precision t,v,tmp,ytildv,thetav,thetav2,c

C     COMMON BLOCKS
      common /p/ ptype,ttype,tconst

      inform = 0

      if ( ptype .eq. 0 ) then

c         This PHR but replacing max(0,a) = 0.5 ( a + |a| )
c         and approximating |a| by sqrt( a^2 + 1/rho ).

          if ( ttype .eq. 1 ) then
              c = 1.0d0
          else if ( ttype .eq. 2 ) then
              c = 1.0d-1
          else if ( ttype .eq. 3 ) then
              c = 1.0d-2
          else if ( ttype .eq. 4 ) then
              c = 1.0d-3
          else if ( ttype .eq. 5 ) then
              c = 1.0d-4
          else if ( ttype .eq. 6 ) then
              c = 1.0d-5
          else if ( ttype .eq. 7 ) then
              c = 1.0d-6
          else if ( ttype .eq. 8 ) then
              c = 0.0d0
          end if

          tmp = mu + rho * y
          p = ( 0.5d0 / rho ) * 
     +        ( 0.25d0 * ( tmp + sqrt( tmp ** 2 + c / rho ) ) ** 2 - 
     +          mu ** 2 )

      else if ( ptype .eq. 1 ) then

          if ( y .ge. 0.0d0 ) then
              p = mu * y + 0.5d0 * rho * y ** 2 + rho ** 2 * y ** 3
          else if ( y .ge. -mu / rho ) then
              p = mu * y + 0.5d0 * rho * y ** 2
          else
              p = -0.5d0 * mu ** 2 / rho
          end if

      else if ( ptype .eq. 2 ) then

          if ( y .ge. 0.0d0 ) then
              p = mu * y + mu * rho * y ** 2 + rho ** 2 * y ** 3 / 6.0d0
          else
              p = mu * y / ( 1.0d0 - rho * y )
          end if

      else if ( ptype .eq. 3 ) then

          if ( y .ge. 0.0d0 ) then
              p = mu * y + mu * rho * y ** 2 
          else
              p = ( mu * y ) / ( 1.0d0 - rho * y )
          end if

      else if ( ptype .eq. 4 .and. ttype .eq. 1 ) then

          if ( mu + rho * y .ge. 0.0d0 ) then
              p = mu * y + 0.5d0 * rho * y ** 2
          else
              p = - 0.5d0 * mu ** 2 / rho
          end if

      else if ( ptype .eq. 4 .and. ttype .eq. 2 ) then

          tmp = sign( 1.0d0, y ) * abs( y ) ** ( 1.0d0 / 3.0d0 ) *
     +          rho ** ( 1.0d0 / 3.0d0 ) + mu
          if ( tmp .ge. 0.0d0 ) then
              p = (3.0d0 * tmp + mu) * y / 4.0d0
          else
              p = - 0.25d0 * mu ** 4 / rho
          end if

      else if ( ptype .eq. 4 .and. ttype .eq. 3 ) then

          if ( mu + sinh( rho * y ) .ge. 0.0d0 ) then
              p = mu * y + ( cosh( rho * y ) - 1.0d0 ) / rho
          else
              tmp = log ( -mu + sqrt( mu ** 2 + 1.0d0 ) )
              p = mu * tmp / rho + ( cosh( tmp ) - 1.0d0 ) / rho
          end if

      else if ((ptype .eq. 5 .or. ptype .eq. 8) .and. ttype .eq. 4) then

          t = rho * y
          v = exp( t ) - 1.0d0
          p = mu * tconst * v / rho

      else if ( ptype .eq. 5 .and. ttype .eq. 5 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              v = - log( 1.0d0 - t )
              p = mu * tconst * v / rho
          else
              v = exp( 2.0d0 * t - 1.0d0 ) + log(2.0d0) - 1.0d0
              p = mu * tconst * v / rho
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 6 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              v = - log( 1.0d0 - t )
              p = mu * tconst * v / rho
          else
              v = 2.0d0 * rho * y ** 2 + ( log( 2.0d0 ) - 0.5d0 ) / rho
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 7 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              v = y / ( 1.0d0 - t )
              p = mu * tconst * v
          else 
              v = exp ( 4.0d0 * t - 2.0d0 )
              p = mu * tconst * v / rho
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 8 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              v = y / ( 1.0d0 - t )
              p = mu * tconst * v
          else
              v = 8.0d0 * rho * y ** 2 - 4.0d0 * y + 1.0d0 / rho
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 9 ) then

          t = rho * y
          if ( t .le. -0.5d0 ) then
              v = - 0.25d0 * log ( -2.0d0 * t ) - 0.375d0
              p = mu * tconst * v / rho
          else
              v = 0.5d0 * rho * y ** 2 + y
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 10 ) then

          t = rho * y
          tmp = ( 1.0d0 + t ) + sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 )
          v = 0.0625d0 * tmp ** 2 + log ( 0.25d0 * tmp ) - 1.0d0
          p = mu * tconst * v / rho

      else if ( ptype .eq. 5 .and. ttype .eq. 11 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              v = max( 0.0d0, t + 0.5d0 ) ** 3 / 6.0d0 - 
     +            1.0d0 / 24.0d0
              p = mu * tconst * v / rho
          else
              v = 0.5d0 * rho * y ** 2
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 12 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              v = exp ( t )
              p = mu * tconst * v / rho
          else
              v = 0.5d0 * exp( 0.5d0 ) * 
     +            ( rho * y ** 2 + y + 1.25d0 / rho )
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 13 ) then

          t = rho * y
          if ( t .le. -0.5d0 ) then
              v = - log ( - t ) - 1.0d0 
              p = mu * tconst * v / rho
          else
              v = 2.0d0 * rho * y ** 2 + 4.0d0 * y + 
     +            ( 0.5d0 + log ( 2.0d0 ) ) / rho
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 14 ) then

          t = rho * y
          if ( t .le. -0.5d0 ) then
              v = - 1.0d0 / t
              p = mu * tconst * v / rho
          else
              v = 8.0d0 * rho * y ** 2 + 12.0d0 * y + 6.0d0 / rho
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 15 ) then

          t = rho * y
          if ( t .le. -1.0d0 ) then
              v = 4.0d0 / ( 1.0d0 - t ) - 2.0d0
              p = mu * tconst * v / rho
          else if ( t .le. -0.25d0 ) then
              v = - log ( - t )
              p = mu * tconst * v / rho
          else
              v = 8.0d0 * rho * y ** 2 + 8.0d0 * y + 
     *            ( 1.5d0 + log ( 4.0d0 ) ) / rho
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 16 ) then

          t = rho * y
          v = 0.5d0 * ( y + sqrt( y ** 2 + 4.0d0 / rho ** 2 ) )
          p = mu * tconst * v

      else if ( ptype .eq. 5 .and. ttype .eq. 17 ) then

          t = rho * y
          if ( t .le. 0.0d0 ) then
              v = log ( 1.0d0 + exp ( t ) )
              p = mu * tconst * v / rho
          else
              v = y + log ( 1.0d0 + exp ( - t ) ) / rho
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 5 .and. ttype .eq. 18 ) then

          t = rho * y
          if ( t .le. 0.0d0 ) then
              v = 0.5d0 * exp ( t )
              p = mu * tconst * v / rho
          else
              v = y + 0.5d0 * exp ( - t ) / rho
              p = mu * tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 4 ) then

          t = rho * mu * y
          v = exp( t ) - 1.0d0
          p = tconst * v / rho

      else if ( ptype .eq. 6 .and. ttype .eq. 5 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              v = -log( 1.0d0 - t )
              p = tconst * v / rho
          else
              v = exp( 2.0d0 * t - 1.0d0 ) + log( 2.0d0 ) - 1.0d0
              p = tconst * v / rho
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 6 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              v = -log( 1.0d0 - t )
              p = tconst * v / rho
          else
              v = 2.0d0 * rho * ( mu * y ) ** 2 + 
     +            ( log( 2.0d0 ) - 0.5d0 ) / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 7 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              v = mu * y / ( 1.0d0 - t )
              p = tconst * v
          else
              v = exp ( 4.0d0 * t - 2.0d0 )
              p = tconst * v / rho
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 8 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              v = mu * y / ( 1.0d0 - t )
              p = tconst * v
          else
              v = 8.0d0 * rho * ( mu * y ) ** 2 - 4.0d0 * mu * y + 
     +            1.0d0 / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 9 ) then

          t = rho * mu * y
          if ( t .le. -0.5d0 ) then
              v = -0.25d0 * log ( -2.0d0 * t ) - 0.375d0
              p = tconst * v / rho
          else
              v = 0.5d0 * rho * ( mu * y ) ** 2 + mu * y
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 10 ) then

          t = rho * mu * y
          tmp = ( 1.0d0 + t ) + sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 )
          v = 0.0625d0 * tmp ** 2 + log ( 0.25d0 * tmp ) - 1.0d0
          p = tconst * v / rho

      else if ( ptype .eq. 6 .and. ttype .eq. 11 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              v = max( 0.0d0, t + 0.5d0 ) ** 3 / 6.0d0 - 1.0d0 / 24.0d0
              p = tconst * v / rho
          else
              v = 0.5d0 * rho * ( mu * y ) ** 2
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 12 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              v = exp( t )
              p = tconst * v / rho
          else
              v = exp( 0.5d0 ) * ( 0.5d0 * rho * ( mu * y ) ** 2 + 
     +            0.5d0 * mu * y + 0.625d0 / rho )
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 13 ) then

          t = rho * mu * y
          if ( t .le. -0.5d0 ) then
              v = - log ( - t ) - 1.0d0 
              p = tconst * v / rho
          else
              v = 2.0d0 * rho * ( mu * y ) ** 2 + 4.0d0 * mu * y + 
     +            ( 0.5d0 + log ( 2.0d0 ) ) / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 14 ) then

          t = rho * mu * y
          if ( t .le. -0.5d0 ) then
              v = - 1.0d0 / t
              p = tconst * v / rho
          else
              v = 8.0d0 * rho * ( mu * y ) ** 2 + 12.0d0 * mu * y + 
     +            6.0d0 / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 15 ) then

          t = rho * mu * y
          if ( t .le. -1.0d0 ) then
              v = 4.0d0 / ( 1.0d0 - t ) - 2.0d0
              p = tconst * v / rho
          else if ( t .le. -0.25d0 ) then
              v = - log ( -t )
              p = tconst * v / rho
          else
              v = 8.0d0 * mu * y * ( rho * mu * y + 1.0d0 ) + 
     +            ( 1.5d0 + log ( 4.0d0 ) ) / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 16 ) then

          t = rho * mu * y
          v = 0.5d0 * ( mu * y + sqrt( ( mu * y ) ** 2 + 
     +        4.0d0 / rho ** 2 ) )
          p = tconst * v

      else if ( ptype .eq. 6 .and. ttype .eq. 17 ) then

          t = rho * mu * y
          if ( t .le. 0.0d0 ) then
              v = log ( 1.0d0 + exp ( t ) )
              p = tconst * v / rho
          else
              v = mu * y + log ( 1.0d0 + exp ( -t ) ) / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 18 ) then

          t = rho * mu * y
          if ( t .le. 0.0d0 ) then
              v = 0.5d0 * exp ( t )
              p = tconst * v / rho
          else
              v = mu * y + 0.5d0 * exp ( -t ) / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 4 ) then

          t = rho * y / mu
          v = exp( t ) - 1.0d0
          p = mu ** 2 * tconst * v / rho

      else if ( ptype .eq. 7 .and. ttype .eq. 5 ) then

          if ( rho * y  .le. 0.5d0 * mu ) then
              t = rho * y / mu
              v = - log( 1.0d0 - t )
              p = mu ** 2 * tconst * v / rho
          else
              t = rho * y / mu
              v = exp( 2.0d0 * t - 1.0d0 ) + log(2.0d0) - 1.0d0
              p = mu ** 2 * tconst * v / rho
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 6 ) then

          if ( rho * y .lt. 0.5d0 * mu ) then
              t = rho * y / mu
              v = - log( 1.0d0 - t )
              p = mu ** 2 * tconst * v / rho
          else
              v = 2.0d0 * rho * y ** 2 + 
     +            ( log( 2.0d0 ) - 0.5d0 ) * mu ** 2 / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 7 ) then

          if ( rho * y .le. 0.5d0 * mu ) then
              v = mu ** 2 * y / ( mu - rho * y )
              p = tconst * v
          else
              t = rho * y / mu
              v = exp ( 4.0d0 * t - 2.0d0 )
              p = mu ** 2 * tconst * v / rho
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 8 ) then

          if ( rho * y .lt. 0.5d0 * mu ) then
              v = mu ** 2 * y / ( mu - rho * y )
              p = tconst * v
          else
              v = 8.0d0 * rho * y ** 2 - 4.0d0 * y * mu + mu ** 2 / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 9 ) then

          if ( rho * y .lt. - 0.5d0 * mu ) then
              t = rho * y / mu
              v = - 0.25d0 * log ( -2.0d0 * t ) - 0.375d0
              p = mu ** 2 * tconst * v / rho
          else
              v = 0.5d0 * rho * y ** 2 + y * mu
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 10 ) then

          t = rho * y / mu
          tmp = ( 1.0d0 + t ) + sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 )
          v = 0.0625d0 * tmp ** 2 + log ( 0.25d0 * tmp ) - 1.0d0
          p = mu ** 2 * tconst * v / rho

      else if ( ptype .eq. 7 .and. ttype .eq. 11 ) then

          if ( rho * y .lt. 0.5d0 * mu ) then
              t = rho * y / mu
              v = max( 0.0d0, t + 0.5d0 ) ** 3 / 6.0d0 - 1.0d0 / 24.0d0
              p = mu ** 2 * tconst * v / rho
          else
              v = 0.5d0 * rho * y ** 2
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 12 ) then

          if ( rho * y .lt. 0.5d0 * mu ) then
              t = rho * y / mu
              v = exp( t )
              p = mu ** 2 * tconst * v / rho
          else
              v = exp( 0.5d0 ) * ( 0.5d0 * rho * y ** 2 + 
     +            0.5d0 * y * mu + 0.625d0 * mu ** 2 / rho )
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 13 ) then

          if ( rho * y .lt. - 0.5d0 * mu ) then
              t = rho * y / mu
              v = - log ( - t ) - 1.0d0 
              p = mu ** 2 * tconst * v / rho
          else
              v = 2.0d0 * rho * y ** 2 + 4.0d0 * y * mu + 
     +            ( 0.5d0 + log ( 2.0d0 ) ) * mu ** 2 / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 14 ) then

          if ( rho * y .lt. - 0.5d0 * mu ) then
              v = - mu ** 3 / ( rho ** 2 * y )
              p = tconst * v
          else
              v = 8.0d0 * rho * y ** 2 + 12.0d0 * y * mu + 
     +            6.0d0 * mu ** 2 / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 15 ) then

          if ( rho * y .lt. - mu ) then
              t = rho * y / mu
              v = 4.0d0 / ( 1.0d0 - t ) - 2.0d0
              p = mu ** 2 * tconst * v / rho
          else if ( rho * y .lt. - 0.25d0 * mu ) then
              t = rho * y / mu
              v = - log ( - t )
              p = mu ** 2 * tconst * v / rho
          else
              v = 8.0d0 * y * ( rho * y + mu ) + 
     +            ( 1.5d0 + log ( 4.0d0 ) ) * mu ** 2 / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 16 ) then

          if ( y .eq. 0.0d0 ) then
              v = 1.0d0
              p = mu ** 2 * tconst * v / rho
          else 
              v = 0.5d0 * mu * ( y + abs( y ) * 
     +            sqrt( 1.0d0 + 4.0d0 * mu ** 2 / ( rho * y ) ** 2 ) )
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 17 ) then

          if ( y .le. 0.0d0 ) then
              t = rho * y / mu
              v = log ( 1.0d0 + exp ( t ) )
              p = mu ** 2 * tconst * v / rho
          else
              t = rho * y / mu
              v = mu * y + mu ** 2 * log ( 1.0d0 + exp ( -t ) ) / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 18 ) then

          if ( y .le. 0.0d0 ) then
              t = rho * y / mu
              v = 0.5d0 * exp ( t )
              p = mu ** 2 * tconst * v / rho
          else
              t = rho * y / mu
              v = mu * y + 0.5d0 * mu ** 2 * exp ( -t ) / rho
              p = tconst * v
          end if

      else if ( ptype .eq. 8 .and. ttype .eq. 5 ) then

          if ( mu .le. 2.0d0 ) then
              t = rho * y + 1.0d0 - 1.0d0 / mu
              if ( t .le. 0.5d0 ) then
                  p = -( log( 1.0d0 / mu - rho * y ) + log( mu ) ) / rho
              else
                  p = ( exp( 2.0d0 * ( rho * y - 1.0d0 / mu + 0.5d0 ) ) 
     +                  + log( 2.0d0 ) - 1.0d0 - log( mu ) ) / rho
              end if
          else 
              t = rho * y + 0.5d0 * ( log ( 0.5d0 * mu ) + 1.0d0 )
              if ( t .le. 0.5d0 ) then
                  p = ( 1.0d0 - log( 1.0d0 - t ) - 0.5d0 * mu - 
     +                  log( 2.0d0 ) ) / rho
              else
                  p = ( 0.5d0 * mu * ( exp( 2.0d0 * rho * y ) - 1.0d0 ) 
     +                  + log( 2.0d0 ) + 1.0d0 ) / rho
              end if
          end if

      else if ( ptype .eq. 8 .and. ttype .eq. 7 ) then

          if ( mu .le. 4.0d0 ) then
              t = rho * y + 1.0d0 - 1.0d0 / sqrt ( mu )
              if ( t .le. 0.5d0 ) then
                  p = mu * y / ( 1.0d0 - sqrt( mu ) * rho * y )
              else 
                  p = ( exp ( 4.0d0 * t - 2.0d0 ) - 
     +                  sqrt ( mu ) + 1.0d0 ) / rho
              end if
          else 
              t = rho * y + 0.25d0 * log ( 0.25d0 * mu ) + 0.5d0
              if ( t .le. 0.5d0 ) then
                  p = ( 1.0d0 / ( 1.0d0 - t ) - 1.0d0 - 
     +                  0.25d0 * mu ) / rho
              else 
                  p = ( exp ( 4.0d0 * rho * y ) - 1.0d0 ) * 
     +                0.25d0 * mu / rho
              end if
          end if

      else if ( ptype .eq. 8 .and. ttype .eq. 9 ) then

          if ( mu .le. 0.5d0 ) then
              t = rho * y - 0.25d0 / mu
              if ( t .le. - 0.5d0 ) then
                  p = log( - 2.0d0 * t ) + log( 2.0d0 * mu ) 
                  p = - 0.25d0 * p / rho
              else 
                  p = ( 0.5d0 * t ** 2 + t - 0.25d0 * 
     +                  ( log( mu ) - 1.5d0 - log ( 0.5d0 ) ) ) / rho
              end if

          else 
              t = rho * y + mu - 1.0d0
              if ( t .le. - 0.5d0 ) then
                  p = log ( - 2.0d0 * t ) + 2.0d0 * mu ** 2 - 0.5d0
                  p = - 0.25d0 * p / rho
              else 
                  p = 0.5d0 * rho * y ** 2 + y * mu
              end if
          end if

      else if ( ( ptype .eq. 8 .and. ttype .eq.  6 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq.  8 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 10 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 11 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 12 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 13 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 14 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 15 ) ) then

          call ytild(mu,zeta,ytildv,inform)
          call theta(rho * y + ytildv,thetav,inform)
          call theta(ytildv,thetav2,inform)
          p = ( thetav - thetav2 ) / rho

      else if ( ( ptype .eq. 9 .and. ttype .eq. 16 ) .or.
     +          ( ptype .eq. 9 .and. ttype .eq. 17 ) .or.
     +          ( ptype .eq. 9 .and. ttype .eq. 18 ) ) then

          call ytild(mu,zeta,ytildv,inform)
          call theta(rho * y + ytildv,thetav,inform)
          call theta(ytildv,thetav2,inform)
          p = zeta * ( thetav - thetav2 ) / rho

      else
          inform = -92
          return
      end if

c     This just true if p = NaN or |p| > 1.0d+20. In this case, 
c     project onto the box [-1.0d+20, 1.0d+20].
      if ( .not. p .gt. -1.0d+20 .or. .not. p .lt. 1.0d+20 ) then
          p = max( -1.0d+20, min( p, 1.0d+20 ) )
      end if

      end

c ****************************************************
c ****************************************************
      subroutine evaldpdy(y,rho,mu,zeta,dpdy,inform)

C     SCALAR ARGUMENTS
      integer inform
      double precision y,rho,mu,zeta,dpdy

C     COMMON SCALARS
      integer ptype,ttype
      double precision tconst

C     LOCAL SCALARS
      double precision t,dtdy,tmp,dtmpdy,aux,ytildv,c

C     COMMON BLOCKS
      common /p/ ptype,ttype,tconst

      inform = 0

      if ( ptype .eq. 0 ) then

          if ( ttype .eq. 1 ) then
              c = 1.0d0
          else if ( ttype .eq. 2 ) then
              c = 1.0d-1
          else if ( ttype .eq. 3 ) then
              c = 1.0d-2
          else if ( ttype .eq. 4 ) then
              c = 1.0d-3
          else if ( ttype .eq. 5 ) then
              c = 1.0d-4
          else if ( ttype .eq. 6 ) then
              c = 1.0d-5
          else if ( ttype .eq. 7 ) then
              c = 1.0d-6
          else if ( ttype .eq. 8 ) then
              c = 0.0d0
          end if

          tmp = mu + rho * y
          aux = sqrt( tmp ** 2 + c / rho )
          dpdy = 0.25d0 * aux * ( 1.0d0 + tmp / aux ) ** 2

      else if ( ptype .eq. 1 ) then

          if ( y .ge. 0.0d0 ) then
              dpdy = mu + rho * y + 3.0d0 * rho ** 2 * y ** 2
          else if ( y .ge. - mu / rho ) then
              dpdy = mu + rho * y
          else
              dpdy = 0.0d0
          end if

      else if ( ptype .eq. 2 ) then

          if ( y .ge. 0.0d0 ) then
              dpdy = mu + 2.0d0 * mu * rho * y + 
     +               ( rho * y ) ** 2 / 3.0d0
          else
              tmp = 1.0d0 - rho * y
              dpdy = ( mu / tmp ) * ( 1.0d0 + ( rho * y ) / tmp )
          end if

      else if ( ptype .eq. 3 ) then

          if ( y .ge. 0.0d0 ) then
              dpdy = mu + 2.0d0 * mu * rho * y 
          else
              tmp = 1.0d0 - rho * y
              dpdy = ( mu / tmp ) * ( 1.0d0 + ( rho * y ) / tmp )
          end if

      else if ( ptype .eq. 4 .and. ttype .eq. 1 ) then

          dpdy = max( 0.0d0, mu + rho * y )

      else if ( ptype .eq. 4 .and. ttype .eq. 2 ) then

          dpdy = max( 0.0d0, sign( 1.0d0, y ) * 
     +           abs( y ) ** ( 1.0d0 / 3.0d0) * 
     +           rho ** ( 1.0d0 / 3.0d0 ) + mu)

      else if ( ptype .eq. 4 .and. ttype .eq. 3 ) then

          dpdy = max ( 0.0d0, mu + sinh(rho * y) )

      else if ((ptype .eq. 5 .or. ptype .eq. 8) .and. ttype .eq. 4) then

          t = rho * y
          dtdy = exp ( t )
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 5 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t )
          else
              dtdy = 2.0d0 * exp( 2.0d0 * t - 1.0d0 )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 6 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t )
          else
              dtdy = 4.0d0 * t
           end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 7 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2
          else
              dtdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 8 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2
          else
              dtdy = 16.0d0 * t - 4.0d0
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 9 ) then

          t = rho * y
          if ( t .le. -0.5d0 ) then
              dtdy = - 0.25d0 / t
          else
              dtdy = 1.0d0 + t
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 10 ) then

          t = rho * y
          aux = sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 )
          tmp = ( 1.0d0 + t ) + aux
          dtmpdy = 1.0d0 + ( 1.0d0 + t ) / aux
          dtdy = 0.125d0 * tmp * dtmpdy + dtmpdy / tmp
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 11 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              dtdy = 0.5d0 * max(0.0d0, t + 0.5d0) ** 2 
          else
              dtdy = t
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 12 ) then

          t = rho * y
          if ( t .le. 0.5d0 ) then
              dtdy = exp ( t )
          else
              dtdy = exp( 0.5d0 ) * ( t + 0.5d0 )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 13 ) then

          t = rho * y
          if ( t .le. -0.5d0 ) then
              dtdy = - 1.0d0 / t
          else
              dtdy = 4.0d0 * t + 4.0d0
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 14 ) then

          t = rho * y
          if ( t .le. -0.5d0 ) then
              dtdy = 1.0d0 / t ** 2
          else
              dtdy = 16.0d0 * t + 12.0d0
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 15 ) then

          t = rho * y
          if ( t .le. -1.0d0 ) then
              dtdy = 4.0d0 / ( 1.0d0 - t ) ** 2
          else if ( t .le. -0.25d0 ) then
              dtdy = - 1.0d0 / t
          else
              dtdy = 16.0d0 * t + 8.0d0
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 16 ) then

          t = rho * y
          if ( abs( t ) .le. 1.0d0 ) then
              tmp = sqrt( t ** 2 + 4.0d0 )
          else
              tmp = abs( t ) * sqrt( 1.0d0 + 4.0d0 / t ** 2 )
          end if
          dtdy = 0.5d0 * ( 1.0d0 + t / tmp )
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 17 ) then

          t = rho * y
          if ( t .le. 0.0d0 ) then
              dtdy = exp ( t ) / ( 1.0d0 + exp ( t ) )
          else
              dtdy = 1.0d0 / ( 1.0d0 + exp ( -t ) )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 5 .and. ttype .eq. 18 ) then

          t = rho * y
          if ( t .le. 0.0d0 ) then
              dtdy = 0.5d0 * exp ( t )
          else
              dtdy = 1.0d0 - 0.5d0 * exp ( -t )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 4 ) then

          t = rho * mu * y
          dtdy = exp ( t )
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 5 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t )
          else
              dtdy = 2.0d0 * exp( 2.0d0 * t - 1.0d0 )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 6 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t )
          else
              dtdy = 4.0d0 * t
           end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 7 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2
          else
              dtdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 8 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2
          else
              dtdy = 16.0d0 * t - 4.0d0
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 9 ) then

          t = rho * mu * y
          if ( t .le. -0.5d0 ) then
              dtdy = - 0.25d0 / ( rho * y )
              dpdy = tconst * dtdy
          else
              dtdy = 1.0d0 + t
              dpdy = mu * tconst * dtdy
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 10 ) then

          t = rho * mu * y
          aux = sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 )
          tmp = ( 1.0d0 + t ) + aux
          dtmpdy = 1.0d0 + ( 1.0d0 + t ) / aux
          dtdy = 0.125d0 * tmp * dtmpdy + dtmpdy / tmp
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 11 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              dtdy = 0.5d0 * max(0.0d0, t + 0.5d0) ** 2 
          else
              dtdy = t
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 12 ) then

          t = rho * mu * y
          if ( t .le. 0.5d0 ) then
              dtdy = exp ( t )
          else
              dtdy = exp( 0.5d0 ) * ( t + 0.5d0 )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 13 ) then

          t = rho * mu * y
          if ( t .le. -0.5d0 ) then
              dtdy = - 1.0d0 / ( rho * y )
              dpdy = tconst * dtdy
          else
              dtdy = 4.0d0 * t + 4.0d0
              dpdy = mu * tconst * dtdy
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 14 ) then

          t = rho * mu * y
          if ( t .le. -0.5d0 ) then
              dtdy = 1.0d0 / ( mu * ( rho * y ) ** 2 )
              dpdy = tconst * dtdy
          else
              dtdy = 16.0d0 * t + 12.0d0
              dpdy = mu * tconst * dtdy
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 15 ) then

          t = rho * mu * y
          if ( t .le. -1.0d0 ) then
              dtdy = 4.0d0 / ( 1.0d0 - t ) ** 2
              dpdy = mu * tconst * dtdy
          else if ( t .le. -0.25d0 ) then
              dtdy = - 1.0d0 / ( rho * y )
              dpdy = tconst * dtdy
         else
              dtdy = 16.0d0 * t + 8.0d0
              dpdy = mu * tconst * dtdy
          end if

      else if ( ptype .eq. 6 .and. ttype .eq. 16 ) then

          t = rho * mu * y
          if ( abs( t ) .le. 1.0d0 ) then
              tmp = sqrt( t ** 2 + 4.0d0 )
          else
              tmp = abs( t ) * sqrt( 1.0d0 + 4.0d0 / t ** 2 )
          end if
          dtdy = 0.5d0 * ( 1.0d0 + t / tmp )
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 17 ) then

          t = rho * mu * y
          if ( t .le. 0.0d0 ) then
              dtdy = exp ( t ) / ( 1.0d0 + exp ( t ) )
          else
              dtdy = 1.0d0 / ( 1.0d0 + exp ( -t ) )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 6 .and. ttype .eq. 18 ) then

          t = rho * mu * y
          if ( t .le. 0.0d0 ) then
              dtdy = 0.5d0 * exp ( t )
          else
              dtdy = 1.0d0 - 0.5d0 * exp ( -t )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 7 .and. ttype .eq. 4 ) then

          t = rho * y / mu
          dtdy = exp ( t )
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 7 .and. ttype .eq. 5 ) then

          t = rho * y / mu
          if ( rho * y .le. 0.5d0 * mu ) then
              dtdy = 1.0d0 / ( 1.0d0 - t )
          else
              dtdy = 2.0d0 * exp( 2.0d0 * t - 1.0d0 )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 7 .and. ttype .eq. 6 ) then

          if ( rho * y .lt. 0.5d0 * mu ) then
              t = rho * y / mu
              dtdy = 1.0d0 / ( 1.0d0 - t )
              dpdy = mu * tconst * dtdy
          else
              dtdy = 4.0d0 * rho * y
              dpdy = tconst * dtdy
           end if

      else if ( ptype .eq. 7 .and. ttype .eq. 7 ) then

          t = rho * y / mu
          if ( rho * y .le. 0.5d0 * mu ) then
              dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2
          else
              dtdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 7 .and. ttype .eq. 8 ) then

          if ( rho * y .lt. 0.5d0 * mu ) then
              t = rho * y / mu
              dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2
              dpdy = mu * tconst * dtdy
          else
              dtdy = 16.0d0 * rho * y - 4.0d0 * mu
              dpdy = tconst * dtdy
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 9 ) then

          if ( rho * y .lt. - 0.5d0 * mu ) then
              dtdy = - 0.25d0 * mu ** 2 / ( rho * y )
              dpdy = tconst * dtdy
          else
              dtdy = mu + rho * y
              dpdy = tconst * dtdy
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 10 ) then

          t = rho * y / mu
          aux = sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 )
          tmp = ( 1.0d0 + t ) + aux
          dtmpdy = 1.0d0 + ( 1.0d0 + t ) / aux
          dtdy = 0.125d0 * tmp * dtmpdy + dtmpdy / tmp
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 7 .and. ttype .eq. 11 ) then

          if ( rho * y .lt. 0.5d0 * mu ) then
              t = rho * y / mu
              dtdy = 0.5d0 * max(0.0d0, t + 0.5d0) ** 2 
              dpdy = mu * tconst * dtdy
          else
              dtdy = rho * y
              dpdy = tconst * dtdy
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 12 ) then

          if ( rho * y .lt. 0.5d0 * mu ) then
              t = rho * y / mu
              dtdy = exp ( t )
              dpdy = mu * tconst * dtdy
          else
              dtdy = exp( 0.5d0 ) * ( rho * y + 0.5d0 * mu )
              dpdy = tconst * dtdy
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 13 ) then

          if ( rho * y .lt. - 0.5d0 * mu ) then
              dtdy = - mu ** 2 / ( rho * y )
              dpdy = tconst * dtdy
          else
              dtdy = 4.0d0 * rho * y + 4.0d0 * mu
              dpdy = tconst * dtdy
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 14 ) then

          if ( rho * y .lt. - 0.5d0 * mu ) then
              dtdy = mu ** 3 / ( rho * y ) ** 2
              dpdy = tconst * dtdy
          else
              dtdy = 16.0d0 * rho * y + 12.0d0 * mu
              dpdy = tconst * dtdy
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 15 ) then

          if ( rho * y .lt. - mu ) then
              t = rho * y / mu
              dtdy = 4.0d0 / ( 1.0d0 - t ) ** 2
              dpdy = mu * tconst * dtdy
          else if ( rho * y .lt. - 0.25d0 * mu ) then
              dtdy = - mu ** 2 / ( rho * y )
              dpdy = tconst * dtdy
          else
              dtdy = 16.0d0 * rho * y + 8.0d0 * mu
              dpdy = tconst * dtdy
          end if

      else if ( ptype .eq. 7 .and. ttype .eq. 16 ) then

          dtdy = 0.5d0 * ( 1.0d0 + sign( 1.0d0, y ) / 
     +           sqrt( 1.0d0 + 4.0d0 * mu ** 2 / ( rho * y ) ** 2 ) )
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 7 .and. ttype .eq. 17 ) then

          t = rho * y / mu
          if ( y .le. 0.0d0 ) then
              dtdy = exp ( t ) / ( 1.0d0 + exp ( t ) )
          else
              dtdy = 1.0d0 / ( 1.0d0 + exp ( -t ) )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 7 .and. ttype .eq. 18 ) then

          t = rho * y / mu
          if ( y .le. 0.0d0 ) then
              dtdy = 0.5d0 * exp ( t )
          else
              dtdy = 1.0d0 - 0.5d0 * exp ( -t )
          end if
          dpdy = mu * tconst * dtdy

      else if ( ptype .eq. 8 .and. ttype .eq. 5 ) then

          if ( mu .le. 2.0d0 ) then
              t = 1.0d0 + rho * y - 1.0d0 / mu
              if ( t .le. 0.5d0 ) then
                  dpdy = 1.0d0 / ( 1.0d0 / mu - rho * y )
              else
                  dpdy = 2.0d0 * 
     +                   exp( 2.0d0 * ( rho * y - 1.0d0 / mu + 0.5d0 ) )
              end if
          else 
              t = rho * y + 0.5d0 * ( log ( 0.5d0 * mu ) + 1.0d0 )
              if ( t .le. 0.5d0 ) then
                  dpdy = 1.0d0 / ( 1.0d0 - t )
              else
                  dpdy = mu * exp( 2.0d0 * rho * y )
              end if
          end if

      else if ( ptype .eq. 8 .and. ttype .eq. 7 ) then

          if ( mu .le. 4.0d0 ) then

              t = rho * y + 1.0d0 - 1.0d0 / sqrt ( mu )
              if ( t .le. 0.5d0 ) then
                  dpdy = 1.0d0 / ( 1.0d0 / sqrt( mu ) - rho * y ) ** 2
              else 
                  dpdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 )
              end if

         else 

              t = rho * y + 0.25d0 * log ( 0.25d0 * mu ) + 0.5d0
              if ( t .le. 0.5d0 ) then
                  dpdy = 1.0d0 / ( 1.0d0 - t ) ** 2
              else 
                  dpdy = exp ( 4.0d0 * rho * y ) * mu
              end if

          end if

      else if ( ptype .eq. 8 .and. ttype .eq. 9 ) then

          if ( mu .le. 0.5d0 ) then

              t = rho * y - 0.25d0 / mu

              if ( t .le. - 0.5d0 ) then
                  dpdy = - 0.25d0 / t
              else 
                  dpdy = t + 1.0d0
              end if

          else 

              tmp = rho * y + mu
              if ( tmp .le. 0.5d0 ) then
                  dpdy = - 0.125d0 / ( tmp - 1.0d0 )
              else 
                  dpdy = tmp
              end if
          end if

      else if ( ( ptype .eq. 8 .and. ttype .eq.  6 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq.  8 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 10 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 11 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 12 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 13 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 14 ) .or.
     +          ( ptype .eq. 8 .and. ttype .eq. 15 ) ) then

          call ytild(mu,zeta,ytildv,inform)
          call evaldtdy(rho * y + ytildv,dtdy,inform)
          dpdy = dtdy

      else if ( ( ptype .eq. 9 .and. ttype .eq. 16 ) .or.
     +          ( ptype .eq. 9 .and. ttype .eq. 17 ) .or.
     +          ( ptype .eq. 9 .and. ttype .eq. 18 ) ) then

          call ytild(mu,zeta,ytildv,inform)
          call evaldtdy(rho * y + ytildv,dtdy,inform)
          dpdy = zeta * dtdy

      else
          inform = -93
          return
      end if

c     This just true if p = NaN or |p| > 1.0d+20. In this case, 
c     project onto the box [-1.0d+20, 1.0d+20].
      if ( .not. dpdy .gt. -1.0d+20 .or. .not. dpdy .lt. 1.0d+20 ) 
     +then
          dpdy = max( -1.0d+20, min( dpdy, 1.0d+20 ) )
      end if

      end 

c ****************************************************
c ****************************************************
      subroutine theta(t,thetav,inform)

C     SCALAR ARGUMENTS
      integer inform
      double precision t,thetav

C     COMMON SCALARS
      integer ptype,ttype
      double precision tconst

C     LOCAL SCALARS
      double precision tpl,tmp

C     COMMON BLOCKS
      common /p/ ptype,ttype,tconst

      inform = 0

      if ( ttype .eq. 1 ) then
          thetav = 0.5d0 * t**2

      else if ( ttype .eq. 2 ) then
          thetav = 0.75d0 * abs(t) ** (4.0d0/3.0d0)

      else if ( ttype .eq. 3 ) then
          thetav = cosh(t) - 1.0d0

      else if ( ttype .eq. 4 ) then
          thetav = exp(t) - 1.0d0

      else if ( ttype .eq. 5 ) then
          if ( t .le. 0.5d0 ) then
              thetav = -log( 1.0d0 - t )
          else
              thetav = exp( 2.0d0 * t - 1.0d0 ) + log(2.0d0) - 1.0d0
          end if

      else if ( ttype .eq. 6 ) then
          if ( t .le. 0.5d0 ) then
              thetav = -log( 1.0d0 - t )
          else
              thetav = 2.0d0 * t ** 2 + log(2.0d0) - 0.5d0
          end if

      else if ( ttype .eq. 7 ) then
          if ( t .le. 0.5d0 ) then
              thetav = 1.0d0 / ( 1.0d0 - t ) - 1.0d0
          else
              thetav = exp ( 4.0d0 * t - 2.0d0 )
          end if

      else if ( ttype .eq. 8 ) then
          if ( t .le. 0.5d0 ) then
              thetav = 1.0d0 / ( 1.0d0 - t ) - 1.0d0
          else
              thetav = 8.0d0 * t ** 2 - 4.0d0 * t + 1.0d0
          end if

      else if ( ttype .eq. 9 ) then
          if ( t .le. -0.5d0 ) then
              thetav = -0.25d0 * log ( -2.0d0 * t ) - 0.375d0
          else
              thetav = 0.5d0 * t ** 2 + t
          end if

      else if ( ttype .eq. 10 ) then
          tpl = 1.0d0 + t
          if ( abs( tpl ) .le. 1.0d0 ) then
              tmp = tpl + sqrt( tpl ** 2 + 8.0d0 )
          else
              tmp = tpl + abs( tpl ) * sqrt( 1.0d0 + 8.0d0 / tpl ** 2 )
          end if
          thetav = 0.0625d0 * tmp ** 2 + log ( 0.25d0 * tmp ) - 1.0d0

      else if ( ttype .eq. 11 ) then
          if ( t .le. 0.5d0 ) then
              thetav = max(0.0d0, t + 0.5d0)**3 / 6.0d0 - 1.0d0 / 24.0d0
          else
              thetav = 0.5d0 * t ** 2
          end if

      else if ( ttype .eq. 12 ) then
          if ( t .le. 0.5d0 ) then
              thetav = exp (t)
          else
              thetav = exp(0.5d0) * (0.5d0*t**2 + 0.5d0*t + 0.625d0)
          end if

      else if ( ttype .eq. 13 ) then
          if ( t .le. -0.5d0 ) then
              thetav = - log ( -t ) - 1.0d0 
          else
              thetav = 2.0d0 * t**2 + 4.0d0 * t + 0.5d0 + log ( 2.0d0 )
          end if

      else if ( ttype .eq. 14 ) then
          if ( t .le. -0.5d0 ) then
              thetav = - 1.0d0 / t
          else
              thetav = 8.0d0 * t**2 + 12.0d0 * t + 6.0d0
          end if

      else if ( ttype .eq. 15 ) then
          if ( t .le. -1.0d0 ) then
              thetav = 4.0d0 / ( 1.0d0 - t ) - 2.0d0
          else if ( t .le. -0.25d0 ) then
              thetav = - log ( -t )
          else
              thetav = 8.0d0 * t**2 + 8.0d0 * t + 1.5d0 + log ( 4.0d0 )
          end if

      else if ( ttype .eq. 16 ) then
          if ( abs( t ) .le. 1.0d0 ) then
              tmp = sqrt( t ** 2 + 4.0d0 )
          else
              tmp = abs( t ) * sqrt( 1.0d0 + 4.0d0 / t ** 2 )
          end if
          thetav = 0.5d0 * ( t + tmp )

      else if ( ttype .eq. 17 ) then
          if ( t .le. 0.0d0 ) then
              thetav = log ( 1.0d0 + exp ( t ) )
          else
              thetav = t + log ( 1.0d0 + exp ( -t ) )
          end if

      else if ( ttype .eq. 18 ) then
          if ( t .le. 0.0d0 ) then
              thetav = 0.5d0 * exp ( t )
          else
              thetav = t + 0.5d0 * exp ( -t )
          end if
      end if

      thetav = thetav * tconst

      end


c ****************************************************
c ****************************************************
      subroutine evaldtdy(t,dtdy,inform)

C     SCALAR ARGUMENTS
      integer inform
      double precision t,dtdy

C     COMMON SCALARS
      integer ptype,ttype
      double precision tconst

C     LOCAL SCALARS
      double precision tpl,tmp,dtmpdc,aux

C     COMMON BLOCKS
      common /p/ ptype,ttype,tconst

      inform = 0

      if ( ttype .eq. 1 ) then
          dtdy = t

      else if ( ttype .eq. 2 ) then
          dtdy = sign(1.0d0,t) * abs(t) ** (1.0d0/3.0d0)

      else if ( ttype .eq. 3 ) then
          dtdy = sinh ( t )

      else if ( ttype .eq. 4 ) then
          dtdy = exp ( t )

      else if ( ttype .eq. 5 ) then
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t )
          else
              dtdy = 2.0d0 * exp( 2.0d0 * t - 1.0d0 )
          end if

      else if ( ttype .eq. 6 ) then
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t )
          else
              dtdy = 4.0d0 * t
           end if

      else if ( ttype .eq. 7 ) then
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2
          else
              dtdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 )
          end if

      else if ( ttype .eq. 8 ) then
          if ( t .le. 0.5d0 ) then
              dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2
          else
              dtdy = 16.0d0 * t - 4.0d0
          end if

      else if ( ttype .eq. 9 ) then
          if ( t .le. -0.5d0 ) then
              dtdy = - 0.25d0 / t
          else
              dtdy = 1.0d0 + t
          end if

      else if ( ttype .eq. 10 ) then
          tpl = 1.0d0 + t
          aux = sqrt( 1.0d0 + 8.0d0 / tpl ** 2 )
          tmp = tpl + abs( tpl ) * aux
          dtmpdc = 1.0d0 + sign(1.0d0,tpl) / aux
          dtdy = 0.125d0 * tmp * dtmpdc + dtmpdc / tmp

      else if ( ttype .eq. 11 ) then
          if ( t .le. 0.5d0 ) then
              dtdy = 0.5d0 * max(0.0d0, t + 0.5d0) ** 2 
          else
              dtdy = t
          end if

      else if ( ttype .eq. 12 ) then
          if ( t .le. 0.5d0 ) then
              dtdy = exp ( t )
          else
              dtdy = exp( 0.5d0 ) * ( t + 0.5d0 )
          end if

      else if ( ttype .eq. 13 ) then
          if ( t .le. -0.5d0 ) then
              dtdy = - 1.0d0 / t
          else
              dtdy = 4.0d0 * t + 4.0d0
          end if

      else if ( ttype .eq. 14 ) then
          if ( t .le. -0.5d0 ) then
              dtdy = 1.0d0 / t ** 2
          else
              dtdy = 16.0d0 * t + 12.0d0
          end if

      else if ( ttype .eq. 15 ) then
          if ( t .le. -1.0d0 ) then
              dtdy = 4.0d0 / ( 1.0d0 - t ) ** 2
          else if ( t .le. -0.25d0 ) then
              dtdy = - 1.0d0 / t
          else
              dtdy = 16.0d0 * t + 8.0d0
          end if

      else if ( ttype .eq. 16 ) then
          if ( abs( t ) .le. 1.0d0 ) then
              tmp = sqrt( t ** 2 + 4.0d0 )
          else
              tmp = abs( t ) * sqrt( 1.0d0 + 4.0d0 / t ** 2 )
          end if
          dtdy = 0.5d0 * ( 1.0d0 + t / tmp )

      else if ( ttype .eq. 17 ) then
          if ( t .le. 0.0d0 ) then
              dtdy = exp ( t ) / ( 1.0d0 + exp ( t ) )
          else
              dtdy = 1.0d0 / ( 1.0d0 + exp ( -t ) )
          end if

      else if ( ttype .eq. 18 ) then
          if ( t .le. 0.0d0 ) then
              dtdy = 0.5d0 * exp ( t )
          else
              dtdy = 1.0d0 - 0.5d0 * exp ( -t )
          end if
      end if

      dtdy = dtdy * tconst

      end


c ****************************************************
c ****************************************************
      subroutine ytild(mu,zeta,ytildv,inform)

C     SCALAR ARGUMENTS
      integer inform
      double precision mu,zeta,ytildv

C     COMMON SCALARS
      integer ptype,ttype
      double precision tconst

C     COMMON BLOCKS
      common /p/ ptype,ttype,tconst

      inform = 0

      if ( ttype .eq. 4 ) then
          ytildv = log ( mu )

      else if ( ttype .eq. 5 ) then
          if ( mu .le. 2.0d0 ) then
              ytildv = 1.0d0 - 1.0d0 / mu
          else 
              ytildv = 0.5d0 * ( log ( 0.5d0 * mu ) + 1.0d0 )
          end if

      else if ( ttype .eq. 6 ) then
          if ( mu .le. 2.0d0 ) then
              ytildv = 1.0d0 - 1.0d0 / mu
          else 
              ytildv = 0.25d0 * mu
          end if

      else if ( ttype .eq. 7 ) then
          if ( mu .le. 4.0d0 ) then
              ytildv = 1.0d0 - 1.0d0 / sqrt ( mu )
          else 
              ytildv = 0.25d0 * log ( 0.25d0 * mu ) + 0.5d0
          end if

      else if ( ttype .eq. 8 ) then
          if ( mu .le. 4.0d0 ) then
              ytildv = 1.0d0 - 1.0d0 / sqrt ( mu )
          else 
              ytildv = 0.0625d0 * mu + 0.25d0
          end if

      else if ( ttype .eq. 9 ) then
          if ( mu .le. 0.5d0 ) then
              ytildv = - 0.25d0 / mu
          else 
              ytildv = mu - 1.0d0
          end if

      else if ( ttype .eq. 10 ) then
          ytildv = 2.0d0 * mu - 1.0d0 / mu - 1.0d0

      else if ( ttype .eq. 11 ) then
          if ( mu .le. 0.5d0 ) then
              ytildv = sqrt ( 2.0d0 * mu ) - 0.5d0
          else 
              ytildv = mu
          end if

      else if ( ttype .eq. 12 ) then
          if ( mu .le. exp (0.5d0) ) then
              ytildv = log ( mu )
          else 
              ytildv = mu * exp( -0.5d0 ) - 0.5d0
          end if

      else if ( ttype .eq. 13 ) then
          if ( mu .le. 2.0d0 ) then
              ytildv = - 1.0d0 / mu
          else 
              ytildv = 0.25d0 * mu - 1.0d0
          end if

      else if ( ttype .eq. 14 ) then
          if ( mu .le. 4.0d0 ) then
              ytildv = - 1.0d0 / sqrt ( mu )
          else 
              ytildv = 0.0625d0 * mu - 0.75d0
          end if

      else if ( ttype .eq. 15 ) then
          if ( mu .le. 1.0d0 ) then
              ytildv = 1.0d0 - 2.0d0 / sqrt ( mu )
          else if ( mu .le. 4.0d0 ) then
              ytildv = - 1.0d0 / mu
          else
              ytildv = 0.0625d0 * mu - 0.5d0
          end if

      else if ( ttype .eq. 16 ) then
          ytildv = ( 2.0d0 * mu - zeta ) / sqrt( mu * zeta - mu ** 2  )

      else if ( ttype .eq. 17 ) then
          ytildv = log ( mu / ( zeta - mu )  )

      else if ( ttype .eq. 18 ) then
          if ( mu .le. 0.5d0 * zeta ) then
              ytildv = log ( 2.0d0 * mu / zeta )
          else
              ytildv = - log ( ( 2.0d0 * ( zeta - mu ) ) / zeta ) 
          end if
      end if

      end

c ****************************************************
c ****************************************************
      subroutine gotc(n,x,m,c,got)

C     SCALAR ARGUMENTS
      logical got
      integer n,m

C     ARRAY ARGUMENTS
      double precision c(n),x(n)

C     PARAMETERS
      integer nmax,mmax
      parameter (nmax=100000,mmax=100000)

C     LOCAL SCALARS
      integer k,i
      logical samep

C     COMMON SCALAR
      integer nc

C     COMMON ARRAYS
      double precision xs(nmax,2),cs(mmax,2)

C     COMMON BLOCKS
      common /constr/xs,cs,nc

      got= .false.

      k= 1

 100  if (k.le.min0(2,nc)) then

          i= 1
          samep= .true.
 110      if (i.le.n.and.samep) then
              if (x(i).ne.xs(i,k)) then
                  samep= .false.
              end if
              i= i + 1
              go to 110
          end if

          if (samep) then
              got= .true.
              do i= 1, m
                  c(i)= cs(i,k)
              end do
              return
          end if

          k= k + 1
          go to 100

      end if

      return

      end


c ****************************************************
c ****************************************************
      subroutine savec(n,x,m,c)

C     SCALAR ARGUMENTS
      integer n,m

C     ARRAY ARGUMENTS
      double precision c(n),x(n)

C     PARAMETERS
      integer nmax,mmax
      parameter (nmax=100000,mmax=100000)

C     LOCAL SCALARS
      integer i

C     COMMON SCALAR
      integer nc

C     COMMON ARRAYS
      double precision xs(nmax,2),cs(mmax,2)

C     COMMON BLOCKS
      common /constr/xs,cs,nc

      if (nc.lt.2) then
          nc= nc + 1
      end if

      do i= 1, n
          xs(i,2)= xs(i,1)
      end do
      do i= 1, m
          cs(i,2)= cs(i,1)
      end do

      do i= 1, n
          xs(i,1)= x(i)
      end do
      do i= 1, m
          cs(i,1)= c(i)
      end do

      return

      end


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

      double precision function norm2s(n,x)

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      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     May 14th, 2003.
c
c     1) The way amax2 and amax2n were being computing may casued a 
c     segmentation fault. Its initialization was changed from infty and
c     -infty to 1.0d+99 and -1.0d+99, respectively. Using infty, when
c     combined with a big trust region radius, the final value of amax2
c     or amax2n may cause the impression that a bound is being attained, 
c     when it is not. "Redundant" ifs inside the amax2 and anax2n calculation
c     were deleted. It should considered the possibility of using two 
c     conatants, namely, bignum = 1.0d+20 and infty = 1.0d+99, instead 
c     of just infty. 

c     Modifications introduced near April 1st 2003 in ocassion of the 
c     PHR and inequality-constraints Augemented Lagrangean methods 
c     development:
c
c     1) The use of iprint was redefined and iprint2 was deleted.
c
c     2) The way to detect no progress in the log of the projected 
c     gradient norm was changed. As it was, "no progess means no
c     reduction in the projected gradient norm over M iterations".
c     But this criterion implicitely assumed that the projected
c     gradient norm msut decrease monotonicaly. Is it is clearly not
c     true, the criterion was changed by a non-monotone decrease
c     criterion. Now, progress means that the projected gradient
c     norm is, at each iteration, smaller than the maximum over the
c     last M iterations. And "no porgess" means the it does not occurs
c     during  not smaller than the 
c
c     3 ) The computation of qamaxn inside cg subroutine was in the wrong
c     place (it was being used before computed) and it was maybe the 
c     reason for which the option nerlyq = .true. never worked properly. 
c     With this correction this option should be tested again.
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,x,d,g,n,xc,gtype,htvtype,hd,xtemp,
     *sterel,steabs,goth,inform)

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

C     ARRAY ARGUMENTS
      integer ind(nind)
      double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n)

C     PARAMETERS
      integer mmax,jcnzemax
      parameter (mmax=100000,jcnzemax=100000)

C     COMMON SCALARS
      integer m
 
C     COMMON ARRAYS
      double precision mubar(mmax),zeta(mmax),rho(mmax)

C     LOCAL SCALARS
      logical got
      integer i,jcnze,flag
      double precision xinfn,dinfn,step,f

C     LOCAL ARRAYS
      integer jcvar(jcnzemax),jcfun(jcnzemax)
      double precision c(mmax),ctemp(mmax),dpdc(mmax),dpdctemp(mmax),
     +                 jcval(jcnzemax)

C     COMMON BLOCKS
      common /auglag/mubar,zeta,rho,m

      inform = 0

c     Compute incremental quotients step

      xinfn = 0.0d0
      dinfn = 0.0d0
      do i = 1,n
          xinfn = max( xinfn, abs( x(i) ) )
          dinfn = max( dinfn, abs( d(i) ) )
      end do

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

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

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

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

C     COMPUTE CONSTRAINTS AT x AND xtemp

      call gotc(n,x,m,c,got)

      if ( .not. got ) then
          call evalfc(n,x,f,m,c,flag)
          if ( flag .ne. 0 ) then
              inform = -90
              return
          end if
      end if

      call evalfc(n,xtemp,f,m,ctemp,flag)
      if ( flag .ne. 0 ) then
          inform = -90
          return
      end if

C     COMPUTE dP/dc

      do i = 1,m
          dpdc(i)     = mubar(i) + rho(i) * c(i)
          dpdctemp(i) = mubar(i) + rho(i) * ctemp(i)
      end do

C     COMPUTE \nabla L(x) = \nabla f(x) + \sum_{i=1}^m dPdc * dcdx

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

      call evalgjc(n,xtemp,m,jcvar,jcfun,jcval,jcnze,flag)
      if ( flag .ne. 0 ) then
          inform = -91
          return
      end if

      do i = 1,jcnze
          if ( jcfun(i) .eq. 0 ) then
              hd(jcvar(i)) = hd(jcvar(i)) + jcval(i)
          else
              if ( dpdc(jcfun(i)) .gt. 0.0d0 ) then
                  hd(jcvar(i)) = 
     +            hd(jcvar(i)) + dpdctemp(jcfun(i)) * jcval(i)
              end if
          end if
      end do

c     Compute incremental quotients

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

      return

      end


