C     =================================================================
C     File: albetra.for
C     =================================================================

C     Last update of any of the component of this file: 
C
C     May 2nd, 2005.

C     =================================================================
C     File: albetrama.f
C     =================================================================

C     Users are encouraged to download periodically updated versions of 
C     this code at the TANGO home page:
C
C     www.ime.usp.br/~egbirgin/tango/

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

C     ALBETRA solves problems of the form
C     -----------------------------------
C
C     min f(x)
C
C     subject to
C
C             c_j(x)  = 0, j in E,
C             c_j(x) <= 0, j in I,
C             l <= x <= u,
C
C     where E is the set of indices of the equality constraints, I is
C     the set of indices of the inequality constraints, and there are
C     n variables and m constraints.
C
C     ALBETRA is an Augmented Lagrangian method that uses BETRA to
C     solve the bound-constrained problems.
C
C     ALBETRA is part of the TANGO Project.
C
C     Visit the TANGO home page in the web:
C
C     www.ime.usp.br/~egbirgin/tango/

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

C     TANGO LICENSE:
C     --------------
C
C     The codes of TANGO are free for non-commercial use. They may be 
C     distributed within a "Research Group" but not distributed to 
C     third parties without the authorization of E. G. Birgin or J. M. 
C     Martínez. The code must not be included in any commercial package 
C     without previous agreement. 
C
C     Every published work that uses ALBETRA must cite:
C
C     R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, 
C     "On Augmented Lagrangian methods with general lower-level 
C     constraints", Technical Report MCDO-050304, Department of Applied 
C     Mathematics, UNICAMP, Brazil, 2005
C
C     and/or
C
C     R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, 
C     "Augmented Lagrangian methods under the Constant Positive Linear 
C     Dependence constraint qualification", Technical Report 
C     MCDO-040806, Department of Applied Mathematics, UNICAMP, Brazil, 
C     2004.
C
C     Every published work that uses BETRA must cite:
C
C     M. Andretta, E. G. Birgin and J. M. Martínez, "Practical 
C     active-set Euclidian trust-region method with spectral projected 
C     gradients for bound-constrained minimization", to appear in 
C     Optimization, 2005.
C
C     (See other related works at the TANGO home page.)

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

C     HOW TO USE ALBETRA TO SOLVE YOUR OWN PROBLEM
C     --------------------------------------------
C
C     You will find below 8 subroutines and that YOU MUST MODIFY to 
C     solve your own problem. The modifications that you must do are:
C
C     1) In the SUBROUTINE INIP you must set the number of variables
C     of your problem (n), the initial point (x), the lower and upper
C     bounds on the variables (l and u), the number of constraints (m)
C     (set m equal to zero if your problems has no constraints), a
C     vector (equatn) that indicates, for each constraint, whether it 
C     is and equality constraint or an inequality (least-or-equal) 
C     constraint, the initial estimation of the Lagrange multipliers 
C     (you can set it equal to zero if you do not have a better 
C     estimative) and the initial penalty parameters (rho).
C
C     2) SUBROUTINES EVALF, EVALG, EVALH, EVALC, EVALJAC AND EVALHC, 
C     MUST also be MODIFIED to evaluate the objective function, its 
C     gradient, its Hessian, the constraints, their gradients and
C     their Hessians, respectively. You will find the way in which 
C     these subroutines must be modified below.
C
C     3) OPTIONALLY, in the SUBROUTINE PARAM you may modify tolerances 
C     and maximum of iterations related to the stopping criteria. There 
C     are also two arguments related to the output in the screen and
C     a parameter related to the derivatives calculation. If you 
C     prefer, leave them with the suggested values.
C
C     Besides the modifications that you must do to solve your own 
C     problem using ALBETRA, ALBETRA has many other arguments to
C     control a large variety of issues of the algorithmic behaviour. 
C     To see how to control these other arguments, see their detailed 
C     description at the begining of subroutine albetra (for arguments 
C     related to Augmented Lagrangians) and at the begining of 
C     subroutine betra (for arguments related to the bound-constraints 
C     internal solver).

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

      subroutine inip(nmax,mmax,n,x,l,u,m,lambda,rho,equatn)

C     This subroutine must set some problem data. For achieving this 
C     objective YOU MUST MODIFY it according to your problem. See below 
C     where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     nmax     integer,
C              maximum allowed number of variables,
C
C     mmax     integer,
C              maximum allowed number of constraints.
C
C     On Return
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              initial point,
C
C     l        double precision l(n),
C              lower bounds on x,
C
C     u        double precision u(n),
C              upper bounds on x,
C
C     m        integer,
C              number of constraints (excluding the bounds),
C
C     lambda   double precision lambda(m),
C              initial estimation of the Lagrange multipliers,
C
C     rho      double precision rho(m),
C              initial penalty parameters,
C
C     equatn   logical equatn(m),
C              for each constraint i, set equatn(i) = .true. if it is an 
C              equality constraint of the form c_i(x) = 0, and set 
C              equatn(i) = .false. if it is an inequality constraint of 
C              the form c_i(x) <= 0.


C     SCALAR ARGUMENTS
      integer m,mmax,n,nmax

C     ARRAY ARGUMENTS
      logical equatn(mmax)
      double precision l(nmax),lambda(mmax),rho(mmax),u(nmax),x(nmax)

C     LOCAL SCALARS
      integer i

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR PROBLEM
C     DATA:
C     ******************************************************************

C     Number of variables

      n = 2

C     Initial point

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

C     Lower and upper bounds

      l(1) = - 1.0d+99
      l(2) =     0.0d0

      u(1) =   1.0d+99
      u(2) =    10.0d0

C     Number of constraints (equalities plus inequalities)

      m = 1

C     Lagrange multipliers approximation. Most users prefer to use the 
C     null initial Lagrange multipliers estimates. However, if the 
C     problem that you are solving is "slightly different" from a 
C     previously solved problem of which you know the correct Lagrange 
C     multipliers, we encourage you to set these multipliers as initial 
C     estimates. Of course, in this case you are also encouraged to use 
C     the solution of the previous problem as initial estimate of the 
C     solution. Similarly, most users prefer to use rho = 10 as initial 
C     penalty parameters. But in the case mentioned above (good 
C     estimates of solution and Lagrange multipliers) larger values of 
C     the penalty parameters (say, rho = 1000) may be more useful. More 
C     warm-start procedures are being elaborated.

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

C     Initial penalty parameters

      do i = 1,m
          rho(i) = 10.0d0
      end do

C     For each constraint i, set equatn(i) = .true. if it is an equality
C     constraint of the form c_i(x) = 0, and set equatn(i) = .false. if 
C     it is an inequality constraint of the form c_i(x) <= 0.

      equatn(1) = .true.

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE INIP.
C     ******************************************************************

      end

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

      subroutine evalf(n,x,f,flag)

C     This subroutine must compute the objective function. For achieving 
C     this objective YOU MUST MODIFY it according to your problem. See 
C     below where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     On Return
C
C     f        double precision,
C              objective function value at x,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the objective 
C              function. (For example, trying to compute the square root 
C              of a negative number, dividing by zero or a very small 
C              number, etc.) If everything was o.k. you must set it 
C              equal to zero.

C     SCALAR ARGUMENTS
      integer flag,n
      double precision f

C     ARRAY ARGUMENTS
      double precision x(n)

C     Objective function

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR OBJECTIVE
C     FUNCTION:
C     ******************************************************************

      flag = 0

      f = - x(1) ** 2 - x(2) ** 2

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALF.
C     ******************************************************************

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalg(n,x,g,flag)

C     This subroutine must compute the gradient vector of the objective 
C     function. For achieving these objective YOU MUST MODIFY it in the 
C     way specified below. However, if you decide to use numerical 
C     derivatives (we dont encourage this option at all!) you dont need
C     to modify evalg.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     On Return
C
C     g        double precision g(n),
C              gradient vector of the objective function evaluated at x,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of any component 
C              of the gradient vector. (For example, trying to compute 
C              the square root of a negative number, dividing by zero or 
C              a very small number, etc.) If everything was o.k. you 
C              must set it equal to zero.

C     SCALAR ARGUMENTS
      integer flag,n

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

C     Gradient vector

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE GRADIENT
C     VECTOR OF YOUR OBJECTIVE FUNCTION: 
C     ******************************************************************

      flag = 0

      g(1) = - 2.0d0 * x(1)
      g(2) = - 2.0d0 * x(2)

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALG. 
C     ******************************************************************
 
      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalh(n,x,ldh,H,flag)

C     This subroutine must compute the Hessian matrx of the objective 
C     function. For achieving these objective YOU MUST MODIFY it in the 
C     way specified below. However, if you decide to use numerical 
C     derivatives (we dont encourage this option at all!) you dont need
C     to modify evalh.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     ldh      integer
C              leading dimension of H,
C
C     On Return
C
C     H        double precision H(ldh,n),
C              Hessian matrix of the objective function evaluated at x,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of any component 
C              of the Hessian matrix. (For example, trying to compute 
C              the square root of a negative number, dividing by zero or 
C              a very small number, etc.) If everything was o.k. you 
C              must set it equal to zero.

C     SCALAR ARGUMENTS
      integer flag,ldh,n

C     ARRAY ARGUMENTS
      double precision H(ldh,n),x(n)

C     LOCAL SCALARS
      integer i,j

C     Hessian matrix

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE HESSIAN
C     MATRIX OF YOUR OBJECTIVE FUNCTION: 
C     ******************************************************************

      flag = 0

      
      do j = 1,n
          do i = 1,n
              H(i,j) = 0.0d0
          end do
      end do

      H(1,1) = - 2.0d0
      H(2,2) = - 2.0d0

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALH. 
C     ******************************************************************
 
      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalc(n,x,i,c,flag)

C     This subroutine must compute the i-th constraint of your problem. 
C     For achieving this objective YOU MUST MOFIFY it according to your 
C     problem. See below the places where your modifications must be 
C     inserted.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     i        integer,
C              index of the constraint to be computed,
C
C     On Return
C
C     c        double precision,
C              i-th constraint evaluated at x,
C
C     flag     integer
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the 
C              constraint. (For example, trying to compute the square 
C              root of a negative number, dividing by zero or a very 
C              small number, etc.) If everything was o.k. you must set 
C              it equal to zero.
 
C     SCALAR ARGUMENTS
      integer i,flag,n
      double precision c

C     ARRAY ARGUMENTS
      double precision x(n)

C     Constraints

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR
C     CONSTRAINTS: 
C     ******************************************************************

      flag = 0

      if ( i .eq. 1 ) then

          c = x(1) - x(2) ** 2
          return

      end if

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALC. 
C     ******************************************************************
 
      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evaljac(n,x,i,indjac,valjac,nnzjac,flag)

C     This subroutine must compute the gradient of the constraint i. For 
C     achieving these objective YOU MUST MODIFY it in the way specified 
C     below.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     i        integer,
C              index of the constraint whose gradient will be computed,
C
C     On Return
C
C     nnzjac   integer,
C              number of perhaps-non-null elements of the computed 
C              gradient,
C
C     indjac   integer indjac(nnzjac),
C              see below,
C
C     valjac   double precision valjac(nnzjac),
C              the non-null value of the partial derivative of the i-th
C              constraint with respect to the indjac(k)-th variable must 
C              be saved at valjac(k).
C
C     flag     integer
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the 
C              constraint. (For example, trying to compute the square 
C              root of a negative number, dividing by zero or a very 
C              small number, etc.) If everything was o.k. you must set 
C              it equal to zero.

C     SCALAR ARGUMENTS
      integer flag,i,m,n,nnzjac

C     ARRAY ARGUMENTS
      integer indjac(n)
      double precision x(n),valjac(n)

C     Sparse gradient vector of the i-th constraint

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE GRADIENTS  
C     OF YOUR CONSTRAINTS: 
C     ******************************************************************

      flag = 0

      if ( i .eq. 1 ) then

          nnzjac = 2

          indjac(1) = 1
          valjac(1) = 1.0d0

          indjac(2) = 2
          valjac(2) = - 2.0d0 * x(2)

          return

      end if

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALJAC. 
C     ******************************************************************
 
      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalhc(n,x,i,ldhc,Hc,flag)

C     This subroutine must compute the Hessian matrix of the i-th 
C     constraint of your problem. For achieving this objective YOU MUST 
C     MOFIFY it according to your problem. See below the places where 
C     your modifications must be inserted.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     i        integer,
C              index of the constraint to be computed,
C
C     ldhc     integer,
C              leading dimension of matrix Hc,
C
C     On Return
C
C     Hc       double precision Hc(ldhc,n),
C              Hessian matrix of the i-th constraint evaluated at x,
C
C     flag     integer
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the Hessian
C              matrix. (For example, trying to compute the square 
C              root of a negative number, dividing by zero or a very 
C              small number, etc.) If everything was o.k. you must set 
C              it equal to zero.
 
C     SCALAR ARGUMENTS
      integer ldhc,i,flag,n

C     ARRAY ARGUMENTS
      double precision Hc(ldhc,n),x(n)

C     LOCAL SCALARS
      integer k,l

C     Hessian matrix of the Constraints

C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE HESSIAN 
C     MATRICES OF YOUR CONSTRAINTS: 
C     ******************************************************************

      flag = 0

      do l = 1,n
          do k = 1,n
              Hc(k,l) = 0.0d0
          end do
      end do

      if ( i .eq. 1 ) then

          Hc(2,2) = - 2.0d0

          return

      end if

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALC. 
C     ******************************************************************
 
      end

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

      subroutine param(epsfeas,epsopt,maxoutit,maxtotit,maxtotfc,iprint,
     +ncomp)

C     This subroutine must set some albetra arguments related to 
C     stopping criteria and output. OPTIONALLY, you may modify the
C     suggested default values.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     This subroutine has no input parameters.
C
C     On Return
C
C     epsfeas  double precision,
C              feasibility tolerance (for the sup-norm of the 
C              constraints),
C
C     epsopt   double precision,
C              optimality tolerance (for the sup-norm of the projected
C              gradient of the Augmented Lagrangian),
C
C     maxoutit integer,
C              maximum number of outer iterations,
C
C     maxtotit integer,
C              maximum total number of inner iterations (total means
C              summing up the inner iterations of each outer iteration),
C
C     maxtotfc integer,
C              maximum total number of Augmented Lagrangian functional
C              evaluations,
C
C     iprint   integer,                
C              controls the ammount of information of the output
C              according to the following convention:
C
C              = 0, no output at all,
C              = 1, information at each outer iteration but without any 
C                   information of how the subproblems are being solved,
C              = 2, the same as 1 plus information of each inner 
C                   iteration,
C              = 3, the same as 2 plus information of the line searches
C                   and the calculation of truncated Newton direction
C                   (using CG) of each inner iteration,
C
C              (In all cases, an output file named solution.txt with
C              the final point, Lagrange mutipliers and penalty 
C              parameters will be generated. Moreover, the same output
C              of the screen will be saved in a file named 
C              albetra.out.)
C
C     ncomp    integer,         
C              every time a vector is printed, just its first ncomp
C              component will be printed.

C     SCALAR ARGUMENTS
      integer iprint,maxoutit,maxtotfc,maxtotit,ncomp
      double precision epsfeas,epsopt
 
C     ******************************************************************
C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET SOME ALBETRA
C     ARGUMENTS RELATED TO DERIVATIVES, STOPPING CRITERIA AND OUTPUT:
C     ******************************************************************

      epsfeas = 1.0d-04
      epsopt  = 1.0d-04

      maxoutit = 50
      maxtotit = 1000000
      maxtotfc = 5 * maxtotit

      iprint = 3
      ncomp  = 5

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE PARAM.
C     ******************************************************************

      end

C     ******************************************************************
C     EXCEPT IF YOU ARE AN EXPERIENT USER, YOUR MODIFICATIONS FINISH 
C     HERE !!!!
C     ******************************************************************

C     ******************************************************************
C     ALL THE INFORMATION THAT ALBETRA PRINTS IN THE SCREEN WILL ALSO
C     BE SAVED IN TEXT FILE CALLED albetra.out. THE SOLUTION, AS WELL
C     AS THE FINAL LAGRANGE MULTIPLIERS APPROXIMATION AND THE FINAL
C     PENALTY PARAMETERS WILL BE SAVED IN A TEXT FILE CALLED 
C     solution.txt.
C     ******************************************************************

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

      program albetrama

C     ******************************************************************
C     NOW, WE SET THE MAXIMUM NUMBER OF VARIABLES (NMAX) AND THE MAXIMUM
C     NUMBER OF CONSTRAINTS (MMAX) THAT THE PROBLEM IS ALLOWED TO HAVE.
C     IF YOUR PROBLEM HAS MORE THAN THE DECLARED VALUES FOR NMAX AND
C     MMAX THEN YOU MUST MODIFY THESE VALUES IN SUCH A WAY THAT THE
C     NUMBER OF VARIABLES (N) OF YOUR PROBLEM BE SMALLER THAN OR EQUAL 
C     TO NMAX AND THE NUMBER OF CONSTRAINTS (M) OF YOUR PROBLEM BE 
C     SMALLER THAN OR EQUAL TO MMAX. 
C
C     ATTENTION: PARAMETERS NMAX AND MMAX ARE ALSO DECLARED IN TWO MORE
C     PLACES, SUBROUTINES EVALNAL AND EVALHD. IF YOU CHANGE THEIR VALUES
C     HERE THEN YOU MUST ALSO DO THE SAME MODIFICATION IN THOSE OTHER
C     PLACES TOO. THE VALUES OF NMAX AND MMAX MUST BE THE SAME IN ALL 
C     PLACES!
C     ******************************************************************

C     PARAMETERS
      integer mmax,nmax
      parameter ( mmax      =  1000 )
      parameter ( nmax      =  1000 )
 
C     LOCAL SCALARS
      integer i,inform,m,n,maxoutit,maxtotit,maxtotfc,outiter,totchcnt,
     +        totfcnt,totgcnt,tothcnt,totiter,iprint,ncomp
      double precision f,nalpsupn,epsfeas,epsopt,snorm
      real time,dtime

C     LOCAL ARRAYS
      integer wi1(nmax),wi2(nmax),wi3(nmax)
      logical equatn(mmax)
      double precision H(nmax,nmax),l(nmax),lambda(mmax),rho(mmax),
     +        u(nmax),wd1(mmax),wd2(mmax),wd3(nmax),wd4(mmax),wd5(mmax),
     +        wd6(nmax),wd7(nmax),wd8(nmax),wd9(nmax),wd10(nmax),
     +        wd11(nmax),wd12(nmax),wd13(nmax),x(nmax)
      real dum(2)

C     EXTERNAL SUBROUTINES
      external dtime,albetra

C     SET UP PROBLEM DATA

      call inip(nmax,mmax,n,x,l,u,m,lambda,rho,equatn)

C     SET SOME ALBETRA ARGUMENTS

      call param(epsfeas,epsopt,maxoutit,maxtotit,maxtotfc,iprint,ncomp)

C     CALL ALBETRA

      open(10, file='albetra.out')

      time = dtime(dum)

      call easyalbetra(n,x,l,u,m,lambda,rho,equatn,epsfeas,epsopt,
     +maxoutit,maxtotfc,maxtotit,iprint,ncomp,f,snorm,nalpsupn,nmax,H,
     +outiter,totiter,totfcnt,totgcnt,tothcnt,totchcnt,inform,wi1,wi2,
     +wi3,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8,wd9,wd10,wd11,wd12,wd13)

      time = dtime(dum)
      time = dum(1)

      close(10)

C     SAVE THE SOLUTION

      open(10, file='solution.txt')

C     Solution point

      write(10,9000)
      do i = 1,n
          write(10,9010) i,x(i)
      end do

C     Lagrange multipliers and penalty parameters

      write(10,9020)
      do i = 1,m
          write(10,9030) i,lambda(i),rho(i)
      end do

      close(10)

C     WRITE STATISTICS

      open(10,file='albetrama-tabline.out')

      write(10,9040) time,inform,n,m,outiter,totiter,totfcnt,totgcnt,
     +               tothcnt,totchcnt,f,snorm,nalpsupn

      close(10)

      stop

C     NON-EXECUTABLE STATEMENTS

 9000 format(/,'FINAL POINT:',//,2X,'INDEX',16X,'X(INDEX)')
 9010 format(I7,1P,D24.16)
 9020 format(/,'FINAL ESTIMATION OF THE LAGRANGE MULTIPLIERS AND ',
     +       'PENALTY PARAMETERS: ',//,2X,'INDEX',11X,'LAMBDA(INDEX)',
     +       15X,'RHO(INDEX)')
 9030 format(I7,1P,D24.16,1X,1P,D24.16)
 9040 format(F8.2,1X,I1,1X,I5,1X,I5,1X,I3,1X,I7,1X,I7,1X,I7,1X,I7,1X,
     +       I7,1X,1P,D24.16,1X,1P,D7.1,1X,1P,D7.1)

      end

C     ******************************************************************
C     STOP HERE ALL YOUR MODIFICATIONS !!!!!!
C     ******************************************************************

C     =================================================================
C     File: albetra.f
C     =================================================================

C     Last update of EASYALBETRA: March 22th, 2005.

      subroutine easyalbetra(n,x,l,u,m,lambda,rho,equatn,epsfc,epsop,
     +maxoutit,maxtotfc,maxtotit,iprint,ncomp,f,snorm,nalpsupn,ldh,H,
     +outiter,totiter,totfcnt,totgcnt,tothcnt,totchcnt,inform,wi1,wi2,
     +wi3,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8,wd9,wd10,wd11,wd12,wd13)

C     SCALAR ARGUMENTS
      integer inform,iprint,ldh,m,maxoutit,maxtotfc,maxtotit,n,
     +        ncomp,outiter,totchcnt,totfcnt,totgcnt,tothcnt,totiter
      double precision epsfc,epsop,f,nalpsupn,snorm

C     ARRAY ARGUMENTS
      integer wi1(n),wi2(n),wi3(n)
      logical equatn(m)
      double precision H(ldh,n),l(n),lambda(m),rho(m),u(n),wd1(m),
     +        wd2(m),wd3(n),wd4(m),wd5(m),wd6(n),wd7(n),wd8(n),wd9(n),
     +        wd10(n),wd11(n),wd12(n),wd13(n),x(n)

C     This subroutine aims to simplify the usage of albetra.

C     CONSTANTS FOR GENERAL USES

C     Steps: h = max( steabs, sterel * abs( x ) ) should be a number 
C     such that h is small ( relatively to x ) and x + h is different 
C     from x. So, h is something that can be used a a step for a finite 
C     differences approximation of a partial derivative relative to x.

C     Epsilons: something smaller than max( epsabs, epsrel * abs( x ) )
C     should be considered as ``zero'' when compared with x. It is used, 
C     for example, to detect that a step taken during a line search is 
C     too small.

C     Infinitys: infrel is a big number that may appear in the 
C     calculations. infabs is a number that should never be reached in 
C     the calculations and is used the represent ``infinite''. Detailed 
C     explanations of how are they used are rather cumbersome.

      double precision steabs,sterel,epsabs,epsrel,infabs,infrel
      parameter ( steabs    = 1.0d-10 )
      parameter ( sterel    = 1.0d-07 )
      parameter ( epsabs    = 1.0d-20 )
      parameter ( epsrel    = 1.0d-10 )
      parameter ( infabs    = 1.0d+99 )
      parameter ( infrel    = 1.0d+20 )

C     LOCAL SCALARS
      character * 3 pptype
      integer maxitncp
      double precision epsopfin,epsopini,epsopcon,rhomult,rhofrac,
     +        lammin,lammax

      lammin    = - 1.0d+20         
      lammax    =   1.0d+20         
   
      rhomult   =    10.0d0
      rhofrac   =     0.5d0
      pptype    =     'ONE'

      epsopini  =     epsop
      epsopfin  =     epsop
      epsopcon  =   1.0d-01

      maxitncp  =         9

C     Call albetra

      call albetra(n,x,l,u,m,lambda,rho,equatn,lammin,lammax,rhomult,
     +rhofrac,pptype,epsfc,epsopini,epsopfin,epsopcon,maxitncp,maxoutit,
     +maxtotfc,maxtotit,iprint,ncomp,steabs,sterel,epsabs,epsrel,
     +infabs,infrel,f,wd1,wd2,snorm,wd3,nalpsupn,ldh,H,outiter,totiter,
     +totfcnt,totgcnt,tothcnt,totchcnt,inform,wd4,wd5,wi1,wi2,wi3,wd6,
     +wd7,wd8,wd9,wd10,wd11,wd12,wd13)

      return

      end

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

C     Last update of ALBETRA or any of its dependencies: 
C
C     May 2nd, 2005.
C
C     See report of modifications at the end of this file.

      subroutine albetra(n,x,l,u,m,lambda,rho,equatn,lammin,lammax,
     +rhomult,rhofrac,pptype,epsfc,epsopini,epsopfin,epsopcon,maxitncp,
     +maxoutit,maxtotfc,maxtotit,iprint,ncomp,steabs,sterel,
     +epsabs,epsrel,infabs,infrel,f,c,sigma,snorm,nal,nalpsupn,ldh,H,
     +outiter,totiter,totfcnt,totgcnt,tothcnt,totchcnt,inform,lambar,
     +sigpre,wi1,wi2,wi3,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8)

C     SCALAR ARGUMENTS
      character * 3 pptype
      integer inform,iprint,ldh,m,maxitncp,maxoutit,maxtotfc,
     +        maxtotit,n,ncomp,outiter,totchcnt,totfcnt,totgcnt,tothcnt,
     +        totiter
      double precision epsabs,epsfc,epsopcon,epsopfin,epsopini,epsrel,f,
     +        infabs,infrel,lammax,lammin,nalpsupn,rhofrac,rhomult,
     +        snorm,steabs,sterel

C     ARRAY ARGUMENTS
      integer wi1(n),wi2(n),wi3(n)
      logical equatn(m)
      double precision c(m),H(ldh,n),l(n),lambar(m),lambda(m),nal(n),
     +        rho(m),sigma(m),sigpre(m),u(n),wd1(n),wd2(n),wd3(n),
     +        wd4(n),wd5(n),wd6(n),wd7(n),wd8(n),x(n)

C     Solves the general-constrained minimization problem
C
C                  Minimize f(x)
C
C                  subject to 
C
C                  c_i(x)  = 0, i \in E,
C
C                  c_i(x) <= 0, i \in I,
C
C                  l <= x <= u
C     
C     using the method of multipliers described in
C 
C     R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, 
C     "On Augmented Lagrangian methods with general lower-level 
C     constraints", Technical Report MCDO-050304, Department of Applied 
C     Mathematics, UNICAMP, Brazil, 2005.
C
C     Parameters of the subroutine:
C
C     On entry:
C
C     n        integer 
C              number of variables
C
C     m        integer
C              number of constraints (equations plus inequations
C
C     x        double precision x(n)
C              initial estimation of the solution
C
C     lambda   double precision lambda(m)
C              initial Lagrange multipliers estimation
C
C     rho      double precision rho(m)
C              initial penalty parameters
C
C     equatn   logical equatn(m)
C              equatn is logical array that, for each constraint, 
C              indicates whether the constraint is an equality constraint
C              (.true.) or an inequality constraint (.false.)
C
C     lammin   double precision
C              see below
C
C     lammax   double precision
C              the Lagrange multipliers approximations used for the 
C              definition of the subproblems must be inside the box
C              [lammin,lammax].
C
C     rhomult  double precision
C              when increased, the penalty parameters will be multiplied
C              by rhomult
C
C     rhofrac  double precision
C              the penalty parameters will be, basically, increased when
C              the infeasibility is larger than a fraction of the 
C              previous infeasibility. This fraction is rhofrac.
C
C     epsfc    double precision 
C              required precision for feasibility and complementarity
C
C     epsopini double precision 
C              required precision for the first subproblem optimality
C              condition (sup-norm)
C
C     epsopfin double precision 
C              required precision for optimality
C
C     epsopcon double precision 
C              the precision used for each subproblem is allways 
C              computed as epsopk = max( epsopfin, epsopcon * 
C              epsop_k-1 )
C
C     maxitncp integer 
C              maximum allowed number of outer iterations without 
C              progress in feasibility
C
C     maxoutit integer 
C              maximum allowed number of outer iterations
C
C     maxtotit integer 
C              maximum allowed number of inner iterations
C
C     maxtotfc integer 
C              maximum allowed number of objective function evaluations
C
C     lambar   double precision lambar(m)
C     sigpre   double precision sigpre(m)
C              working vectors
C
C     On return:
C
C     x        double precision x(n)
C              final solution estimation
C
C     lambda   double precision lambda(m)
C              final Lagrange multipliers estimation
C
C     rho      double precision rho(m)
C              final penalty parameters
C
C     f        double precision
C              objective functional value at the solution
C
C     c        double precision c(m)
C              constraints at the solution
C
C     s        double precision c(m)
C              constraints values at the solution
C
C     snorm    double precision 
C              constraints violation is defined as abs( c(i) for equality 
C              constraints and abs( max( c(i), - lambar(i) / rho(i) ) ) 
C              for inequalities, where lambar(i) is the projection of the 
C              Lagrange multiplier approximation lambda(i) onto the 
C              interval [lammin,lammax]. snorm represents the maximum of 
C              these violations.
C
C     nal      double precision nal(n)
C              gradient vector of the Augmented Lagrangian at the 
C              solution
C
C     nalpsupn double precision 
C              sup-norm of the continuous projected gradient of the 
C              Augmented Lagrangian at the solution
C
C     outiter  integer 
C              number of outer iterations
C
C     totiter  integer 
C              total number of inner iterations
C
C     totfcnt  integer 
C              total number of objective functional evaluations
C
C     totgcnt  integer 
C              total number of gradient evaluations
C
C     totcgcnt integer 
C              total number of Conjugate Gradients iterations
C
C     inform   integer
C
C              This output parameter tells what happened in this 
C              subroutine, according to the following conventions:
C 
C              0 = convergence with feasibility, optimality and 
C                  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 
C                  progress'', that means that the current point is 
C                  infeasible (the constraints violation is larger than 
C                  the tolerance epsfc) and the constraint violation has 
C                  not even simple decrease during maxitncp consecutive 
C                  iterations; in this case, the problem may be 
C                  infeasible;
C
C            -90 = error in evalfc subroutine;
C
C            -91 = error in evalgjc subroutine;
C 
C            -92 = error in evalhd subroutine.

C     LOCAL SCALARS
      integer chcnt,fcnt,gcnt,hcnt,i,icrit,iter,mprint,nprint,maxit,
     +        maxfc,rhoincr
      double precision al,snormprev,epsopk,nalpi

C     SET UP SOME CONSTANTS

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

C     INTIALIZATION

      icrit    = 0

      outiter  = 0
      totiter  = 0
      totfcnt  = 0
      totgcnt  = 0
      tothcnt  = 0
      totchcnt = 0

      if ( m .ne. 0 ) then
          epsopk = epsopini
      else
          epsopk = epsopfin
      end if

C     COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS

      call evalfc(n,x,f,m,c,inform)
      totfcnt = totfcnt + 1

      if ( inform .lt. 0 ) then

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

          go to 500
      end if

C     COMPUTE LAGRANGE MULTIPLIERS

      do i = 1,m
          lambar(i) = max( lammin, min( lammax, lambda(i) ) )
      end do

C     COMPUTE COMPLEMENTARITY AND FEASIBILITY VIOLATIONS
 
      snormprev = infabs

      snorm = 0.0d0
      do i = 1,m
          if ( equatn(i) ) then
              sigma(i) = c(i)
          else
              sigma(i) = max( c(i), - lambar(i) / rho(i) )
          end if

          snorm = max( snorm, abs( sigma(i) ) )
      end do

C     COMPUTE CONTINUOUS PROJECTED GRADIENT NORM

      call evalnal(n,x,m,lambar,rho,equatn,nal,inform)
      totgcnt = totgcnt + 1

      if ( inform .lt. 0 ) then

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

          go to 500
      end if

      nalpsupn = 0.0d0
      do i = 1,n
          nalpi    = max( l(i), min( x(i) - nal(i), u(i) ) ) - x(i)
          nalpsupn = max( nalpsupn, abs( nalpi ) )
      end do

C     PRINT INITIAL INFORMATION

      if  ( iprint .ge. 1 ) then
          write(*,  900) n,m
          write(*,  910) nprint,(l(i),i=1,nprint)
          write(*,  920) nprint,(u(i),i=1,nprint)

          write(10, 900) n,m
          write(10, 910) nprint,(l(i),i=1,nprint)
          write(10, 920) nprint,(u(i),i=1,nprint)
      end if

C     MAIN LOOP

 100  continue

C     PRINT ITERATION INFORMATION

      if ( iprint .ge. 1 ) then
          write(*,  970) outiter
          write(*,  980) nprint,(x(i),i=1,nprint)
          write(*,  990) mprint,(lambda(i),i=1,mprint)
          write(*, 1000) mprint,(rho(i),i=1,mprint)
          write(*, 1020) f,snorm,nalpsupn,totiter,totfcnt,totgcnt

          write(10, 970) outiter
          write(10, 980) nprint,(x(i),i=1,nprint)
          write(10, 990) mprint,(lambda(i),i=1,mprint)
          write(10,1000) mprint,(rho(i),i=1,mprint)
          write(10,1020) f,snorm,nalpsupn,totiter,totfcnt,totgcnt
      end if

C     SAVING INTERMEDIATE DATA FOR CRASH REPORT

      open(20, file='albetra-tabline.out')

      write(20,3000) n,m,outiter,totiter,totfcnt,totgcnt,tothcnt,
     +               totchcnt,f,snorm,nalpsupn

      close(20)

C     TEST STOPPING CRITERIA

      if ( snorm .le. epsfc .and. nalpsupn .le. epsopfin ) then
        ! THE POINT IS FEASIBLE AND OPTIMAL

          inform = 0

          if ( iprint .ge. 1 ) then
              write(*, 1060) inform,epsfc,epsfc,epsopfin
              write(10,1060) inform,epsfc,epsfc,epsopfin
          end if

          go to 500
      end if

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

          if ( iprint .ge. 1 ) then
              write(*, 1070) inform,maxoutit
              write(10,1070) 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(*, 1080) inform,maxtotit
              write(10,1080) 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(*, 1090) inform,maxtotfc
              write(10,1090) inform,maxtotfc
          end if

          go to 500
      end if

C     TEST IF THE INFEASIBILITY DECREASED SUFFICIENTLY

      if ( snorm .gt. epsfc .and. snorm .ge. snormprev ) then

          icrit = icrit + 1

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

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

              go to 500
          end if

      else
          icrit = 0
      end if

C     DO AN OUTER ITERATION

      outiter = outiter + 1      

C     OPTIMALITY REQUERIMENT FOR THE SUBPROBLEM

      if ( m .ne. 0 ) then
          if ( outiter .eq. 1 ) then
              epsopk = epsopini
          else
              epsopk = max( epsopfin, epsopk * epsopcon )
          end if
      else
          epsopk = epsopfin
      end if

      maxit =       5000
      maxfc = 10 * maxit

C     CALL THE INNER-SOLVER

      call easybetra(n,x,l,u,m,lambar,rho,equatn,epsopk,maxit,maxfc,
     +iprint,ncomp,sterel,steabs,epsrel,epsabs,infrel,infabs,al,nal,
     +nalpsupn,ldh,H,iter,fcnt,gcnt,hcnt,chcnt,inform,wi1,wi2,wi3,wd1,
     +wd2,wd3,wd4,wd5,wd6,wd7,wd8)

      totiter  = totiter  + iter
      totfcnt  = totfcnt  + fcnt
      totgcnt  = totgcnt  + gcnt
      tothcnt  = tothcnt  + hcnt
      totchcnt = totchcnt + chcnt

      if ( inform .lt. 0 ) then

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

          go to 500
      end if

C     COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS

      call evalfc(n,x,f,m,c,inform)
      totfcnt = totfcnt + 1

      if ( inform .lt. 0 ) then

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

          go to 500
      end if

C     UPDATE LAGRANGEAN MULTIPLIERS APPROXIMATION 

      do i = 1,m
          call evaldpdy(c(i),rho(i),lambar(i),equatn(i),lambda(i))
          lambar(i) = max( lammin, min( lammax, lambda(i) ) )
      end do

C     COMPUTE COMPLEMENTARITY AND FEASIBILITY VIOLATIONS
 
      snormprev = snorm

      snorm = 0.0d0
      do i = 1,m
          sigpre(i) = sigma(i)

          if ( equatn(i) ) then
              sigma(i) = c(i)
          else
              sigma(i) = max( c(i), - lambar(i) / rho(i) )
          end if

          snorm = max( snorm, abs( sigma(i) ) )
      end do

C     UPDATE PENALTY PARAMETERS

      if ( pptype .eq. 'ONE' ) then

          if ( snorm .gt. rhofrac * snormprev ) then

              do i = 1,m
                  rho(i) = rhomult * rho(i)
              end do

              if ( iprint .ge. 1 ) then
                  write(*, 1030) rhofrac * snormprev,rhomult
                  write(10,1030) rhofrac * snormprev,rhomult
              end if

          else

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

          end if

      else ! if ( pptype .eq. 'DIF' ) then

          rhoincr = 0
          do i = 1,m
              if ( abs(sigma(i)) .gt. rhofrac * abs(sigpre(i)) ) then
                  rho(i)  = rhomult * rho(i)
                  rhoincr = rhoincr + 1
              end if
          end do

          if ( iprint .ge. 1 ) then

              if ( rhoincr .ne. 0 ) then
                  write(*, 1040) rhoincr,m,rhomult
                  write(10,1040) rhoincr,m,rhomult
              else
                  write(*, 1050)
                  write(10,1050)
              end if

          end if

      end if

C     ITERATE

      go to 100

C     STOP

 500  continue

      return

C     NON-EXECUTABLE STATEMENTS 

  900 format( /,' Entry to ALBETRA.',/,
     +        /,' Number of variables: ',I7,
     +          ' Number of constraints: ',I7)
  910 format( /,' Lower bounds (first ',I7, ' components): ',
     +        /,6(1X,1PD11.4))
  920 format( /,' Upper bounds (first ',I7, ' components): ',
     +        /,6(1X,1PD11.4))
  970 format( /,' ALBETRA outer iteration: ',I7)
  980 format( /,' Current point (first ',I7, ' components): ',
     +        /,6(1X,1PD11.4))
  990 format( /,' Updated Lagrange multipliers (first ',I7,
     +          ' components): ',
     +        /,6(1X,1PD11.4))
 1000 format( /,' Updated penalty parameters (first ',I7,
     +          ' components): ',
     +        /,6(1X,1PD11.4))
 1020 format( /,' Objective function value',38X,' = ',1PD11.4,
     +        /,' Sup-norm of the constraints',35X,' = ',1PD11.4,
     +        /,' Sup-norm of the projected gradient of the',
     +          ' Lagrangian',10X,' = ',1PD11.4,/,
     +        /,' Up-to-now number of inner iterations',26X,' = ',4X,I7,
     +        /,' Up-to-now number of Augmented Lagrangian function',
     +          ' evaluations',1X,' = ',4X,I7,
     +        /,' Up-to-now number of Augmented Lagrangian gradient',
     +          ' evaluations',1X,' = ',4X,I7)
 1030 format( /,' The desired infeasibility improvement was not',
     +          ' achieved (',1PD11.4,').',
     +        /,' The penalty parameter will be increased multiplying',
     +          ' by rhofrac = ',1PD11.4,'.')
 1040 format( /,' The desired infeasibility improvement was not',
     +          ' achieved in ',I7,' over ',I7,
     +        /,' constraints.',
     +        /,' Penalty parameters will be increased multiplying',
     +          ' by rhofrac = ',1PD11.4,'.')
 1050 format( /,' Desired feasibility improvement was achieved,',
     +        /,' so, penalty parameter(s) will not be modified.')
 1060 format( /,' Flag of ALBETRA = ',I2,
     +          ' (Convergence with Sup-norm of the Lagrangian',
     +          ' projected',/,' gradient smaller than ',1PD11.4,',',
     +          ' constraints Sup-norm smaller than',
     +        /,' ',1PD11.4,', and largest complementarity violation',
     +          ' smaller than ',1PD11.4,')')
 1070 format( /,' Flag of ALBETRA = ',I2,
     +          ' (It was exceeded the maximum allowed number of outer',
     +        /,' iterations (maxoutit =',I7,').)')
 1080 format( /,' Flag of ALBETRA = ',I2,
     +          ' (It was exceeded the maximum allowed number of inner',
     +        /,' iterations (maxinnit =',I7,').)')
 1090 format( /,' Flag of ALBETRA = ',I2,
     +          ' (It was exceeded the maximum allowed number of',
     +        /,' functional evaluations (maxfc =',I7,').)')
 1100 format( /,' Flag of ALBETRA = ',I2,
     +          ' (Constraint violations have not decreased',
     +          ' substantially over',I2,' outer iterations. ',
     +        /,' Problem possibly infeasible.)')
 2000 format( /,' Flag of ALBETRA = ',I3,' Fatal Error',/,
     +        /,' The following codes means: ',/,
     +        /,' -90 : error in evalf   subroutine', 
     +        /,' -91 : error in evalc   subroutine', 
     +        /,' -92 : error in evalg   subroutine', 
     +        /,' -93 : error in evaljac subroutine',/) 
 3000 format(I5,1X,I5,1X,I3,1X,I7,1X,I7,1X,I7,1X,I7,1X,I7,1X,1P,D24.16,
     +       1X,1PD7.1,1X,1PD7.1)

      end

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

      subroutine evalfc(n,x,f,m,c,inform)

C     This subroutine computes the objective function and the 
C     constraints.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     m        integer,
C              number of constraints,
C
C     On Return
C
C     f        double precision,
C              objective function value at x,
C
C     c        double precision c(m),
C              constraints at x,
C
C     inform   integer
C              = 0 means that the evaluation was successfuly done,
C              < 0 means that some error occured during the evaluation.

C     SCALAR ARGUMENTS
      double precision f
      integer inform,m,n

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

C     LOCAL SCALARS
      integer flag,i

      inform = 0

C     COMPUTE OBJECTIVE FUNTION

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

C     COMPUTE CONSTRAINTS

      do i = 1,m

C         COMPUTE THE i-TH CONSTRAINT
          call evalc(n,x,i,c(i),flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

      end do

      end

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

      subroutine evalal(n,x,m,lambda,rho,equatn,al,inform)

C     This subroutine computes the Augmented Lagrangian function.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              point at which the Augmented Lagrangian function will be 
C              evaluated.
C
C     m        integer,
C              number of constraints (equalities plus inequalities),
C
C     lambda   double precision lambdae(m),
C              current estimation of the Lagrange multipliers,
C
C     rho      double precision rho(m)
C              penalty parameters,
C
C     equatn   logical equatn(m)
C              equatn is logical array that, for each constraint, 
C              indicates whether the constraint is an equality 
C              constraint (.true.) or an inequality constraint 
C              (.false.)
C
C     On Return
C
C     al       double precision al,
C              value of the Augmented Lagrangian function,
C
C     inform   integer
C              = 0 means that the evaluation was successfuly done,
C              < 0 means that some error occured during the evaluation.

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

C     ARRAY ARGUMENTS
      logical equatn(m)
      double precision lambda(m),rho(m),x(n)

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

      inform = 0

C     COMPUTE OBJECTIVE FUNTION

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

C     COMPUTES AL(x) = f(x) + \sum_j=1^m P(c(j),rho(j),lambda(j))

      al = f

      do i = 1,m

C         COMPUTE i-TH CONSTRAINT
          call evalc(n,x,i,c,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

C         ADD P(c(i),rho(i),lambda(i))
          call evalp(c,rho(i),lambda(i),equatn(i),p)
          al = al + p

      end do

      end

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

      subroutine evalp(y,rho,lambda,equatn,p)

C     SCALAR ARGUMENTS
      logical equatn
      double precision lambda,p,rho,y

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

      return

      end

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

      subroutine evalnal(n,x,m,lambda,rho,equatn,nal,inform)

C     This subroutine computes the gradient of the Augmented Lagrangian 
C     function.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              point at which the gradient will be evaluated.
C
C     m        integer,
C              number of constraints,
C
C     lambda   double precision lambda(m),
C              current estimation of the Lagrange multipliers,
C
C     rho      double precision rho(m)
C              penalty parameters,
C
C     equatn   logical equatn(m)
C              equatn is logical array that, for each constraint, 
C              indicates whether the constraint is an equality 
C              constraint (.true.) or an inequality constraint 
C              (.false.)
C
C     On Return
C
C     nal      double precision nal(n),
C              gradient of the Augmented Lagrangian function,
C
C     inform   integer
C              = 0 means that the evaluation was successfuly done,
C              < 0 means that some error occured during the evaluation.

C     ******************************************************************
C     HERE WE SET THE MAXIMUM NUMBER OF VARIABLES (NMAX) AND THE MAXIMUM
C     NUMBER OF CONSTRAINTS (MMAX) THAT THE PROBLEM IS ALLOWED TO HAVE.
C     IF YOUR PROBLEM HAS MORE THAN THE DECLARED VALUES FOR NMAX AND
C     MMAX THEN YOU MUST MODIFY THESE VALUES IN SUCH A WAY THAT THE
C     NUMBER OF VARIABLES (N) OF YOUR PROBLEM BE SMALLER THAN OR EQUAL 
C     TO NMAX AND THE NUMBER OF CONSTRAINTS (M) OF YOUR PROBLEM BE 
C     SMALLER THAN OR EQUAL TO MMAX. 
C
C     ATTENTION: PARAMETERS NMAX AND MMAX ARE ALSO DECLARED IN TWO MORE
C     PLACES, THE MAIN PROGRAM AND SUBROUTINE EVALHD. IF YOU CHANGE 
C     THEIR VALUES HERE THEN YOU MUST ALSO DO THE SAME MODIFICATION IN 
C     THOSE OTHER PLACES TOO. THE VALUES OF NMAX AND MMAX MUST BE THE 
C     SAME IN ALL PLACES!
C     ******************************************************************

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

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF PARAMETERS NMAX AND MMAX. 
C     ******************************************************************
 
C     SCALAR ARGUMENTS
      integer inform,m,n

C     ARRAY ARGUMENTS
      logical equatn(m)
      double precision lambda(m),nal(n),rho(m),x(n)

C     LOCAL SCALARS
      integer flag,i,j,jcnnz
      double precision c,dpdc

C     LOCAL ARRAYS
      integer jcvar(nmax)
      double precision jcval(nmax)
 
      inform = 0

C     COMPUTE THE GRADIENT OF THE OBJECTIVE FUNCTION

      call evalg(n,x,nal,flag)
      if ( flag .ne. 0 ) then
          inform = - 92
          return
      end if

C     COMPUTE \nabla L(x) = \nabla f(x) + \sum_{j=1}^m dPdcj * dcjdx

      do j = 1,m

C         COMPUTE THE j-TH CONSTRAINT
          call evalc(n,x,j,c,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

C         COMPUTE dP/dc
          call evaldpdy(c,rho(j),lambda(j),equatn(j),dpdc)

          if ( dpdc .ne. 0.0d0 ) then

C             COMPUTE THE GRADIENT OF THE j-TH CONSTRAINT
              call evaljac(n,x,j,jcvar,jcval,jcnnz,flag)
              if ( flag .ne. 0 ) then
                  inform = - 93
                  return
              end if

C             ADD dPdc * dcdx
              do i = 1,jcnnz
                  nal(jcvar(i)) = nal(jcvar(i)) + dpdc * jcval(i)
              end do

          end if

      end do

      end

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

      subroutine evaldpdy(y,rho,lambda,equatn,dpdy)

C     SCALAR ARGUMENTS
      logical equatn
      double precision y,rho,lambda,dpdy

      if ( equatn ) then
          dpdy = lambda + rho * y
      else
          dpdy = max( 0.0d0, lambda + rho * y )
      end if

      return

      end

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

      subroutine evalhal(n,x,m,lambda,rho,equatn,ldhal,Hal,diag,inform)

C     This subroutine computes the Hessian of the Augmented Lagrangian 
C     function.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              point at which the gradient will be evaluated.
C
C     m        integer,
C              number of constraints,
C
C     lambda   double precision lambda(m),
C              current estimation of the Lagrange multipliers,
C
C     rho      double precision rho(m)
C              penalty parameters,
C
C     equatn   logical equatn(m)
C              equatn is logical array that, for each constraint, 
C              indicates whether the constraint is an equality 
C              constraint (.true.) or an inequality constraint 
C              (.false.)
C
C     ldhal    integer,
C              leading dimension of matrix Hal,
C
C     On Return
C
C     Hal      double precision Hal(ldhal,n),
C              Hessian of the Augmented Lagrangian function,
C
C     diag     double precision diag(n),
C              diagonal of the Hessian matrix,
C
C     inform   integer
C              = 0 means that the evaluation was successfuly done,
C              < 0 means that some error occured during the evaluation.

C     ******************************************************************
C     HERE WE SET THE MAXIMUM NUMBER OF VARIABLES (NMAX) AND THE MAXIMUM
C     NUMBER OF CONSTRAINTS (MMAX) THAT THE PROBLEM IS ALLOWED TO HAVE.
C     IF YOUR PROBLEM HAS MORE THAN THE DECLARED VALUES FOR NMAX AND
C     MMAX THEN YOU MUST MODIFY THESE VALUES IN SUCH A WAY THAT THE
C     NUMBER OF VARIABLES (N) OF YOUR PROBLEM BE SMALLER THAN OR EQUAL 
C     TO NMAX AND THE NUMBER OF CONSTRAINTS (M) OF YOUR PROBLEM BE 
C     SMALLER THAN OR EQUAL TO MMAX. 
C
C     ATTENTION: PARAMETERS NMAX AND MMAX ARE ALSO DECLARED IN TWO MORE
C     PLACES, THE MAIN PROGRAM AND SUBROUTINE EVALHD. IF YOU CHANGE 
C     THEIR VALUES HERE THEN YOU MUST ALSO DO THE SAME MODIFICATION IN 
C     THOSE OTHER PLACES TOO. THE VALUES OF NMAX AND MMAX MUST BE THE 
C     SAME IN ALL PLACES!
C     ******************************************************************

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

C     ******************************************************************
C     STOP HERE YOUR MODIFICATIONS OF PARAMETERS NMAX AND MMAX. 
C     ******************************************************************
 
C     SCALAR ARGUMENTS
      integer inform,ldhal,m,n

C     ARRAY ARGUMENTS
      logical equatn(m)
      double precision diag(n),Hal(ldhal,n),lambda(m),rho(m),x(n)

C     LOCAL SCALARS
      integer flag,i,j,k,l,jcnnz
      double precision c,dpdc

C     LOCAL ARRAYS
      integer jcvar(nmax)
      double precision jcval(nmax),Hc(nmax,nmax)
 
      inform = 0

C     COMPUTE THE HESSIAN OF THE OBJECTIVE FUNCTION

      call evalh(n,x,ldhal,Hal,flag)
      if ( flag .ne. 0 ) then
          inform = - 94
          return
      end if

C     COMPUTE \nabla^2 L(x) = \nabla^2 f(x) + \sum_{j=1}^m dPdcj * dcjdx

      do j = 1,m

C         COMPUTE THE j-TH CONSTRAINT
          call evalc(n,x,j,c,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

C         COMPUTE dP/dc
          call evaldpdy(c,rho(j),lambda(j),equatn(j),dpdc)

          if ( equatn(j) .or. dpdc .gt. 0.0d0 ) then

              if ( dpdc .ne. 0.0d0 ) then

C                 COMPUTE THE HESSIAN OF THE j-TH CONSTRAINT
                  call evalhc(n,x,j,nmax,Hc,flag)
                  if ( flag .ne. 0 ) then
                      inform = - 95
                      return
                  end if

C                 ADD dPdc * \nabla^2 c_j(x)
                  do k = 1,n
                      do l = 1,n
                          Hal(l,k) = Hal(l,k) + dpdc * Hc(l,k)
                      end do
                  end do

              end if

C             COMPUTE THE GRADIENT OF THE j-TH CONSTRAINT
              call evaljac(n,x,j,jcvar,jcval,jcnnz,flag)
              if ( flag .ne. 0 ) then
                  inform = - 93
                  return
              end if

C             ADD rho_j * \nabla c_j(x)^T \nabla c_j(x)
              do k = 1,jcnnz
                  do l = 1,jcnnz
                      Hal(jcvar(l),jcvar(k)) = Hal(jcvar(l),jcvar(k)) + 
     +                rho(j) * jcval(l) * jcval(k)
                  end do
              end do

          end if

      end do

C     SET THE DIAGONAL VECTOR

      do i = 1,n
          diag(i) = Hal(i,i)
      end do

      end

C     =================================================================
C     File: betra.f
C     =================================================================

C     Last update of EASYBETRA: April 4th, 2005.

      subroutine easybetra(n,x,l,u,m,lambda,rho,equatn,epsgpsn,maxit,
     +maxfc,iprint,ncomp,sterel,steabs,epsrel,epsabs,infrel,infabs,f,g,
     +gpsupn,ldh,H,iter,fcnt,gcnt,hcnt,chcnt,inform,wi1,wi2,wi3,wd1,wd2,
     +wd3,wd4,wd5,wd6,wd7,wd8)

C     SCALAR ARGUMENTS
      integer chcnt,fcnt,gcnt,hcnt,m,maxfc,maxit,n,ncomp,inform,iprint,
     +        iter,ldh
      double precision epsabs,epsgpsn,epsrel,f,gpsupn,infabs,infrel,
     +        steabs,sterel

C     ARRAY ARGUMENTS
      integer wi1(n),wi2(n),wi3(n)
      logical equatn(m)
      double precision g(n),H(ldh,n),l(n),lambda(m),rho(m),u(n),wd1(n),
     +        wd2(n),wd3(n),wd4(n),wd5(n),wd6(n),wd7(n),wd8(n),x(n)

C     This subroutine aims to simplify the use of BETRA. For this 
C     purpose it gives values to most of the BETRA arguments and 
C     leaves to the user those arguments which he/she may would like to 
C     set by him/herself.

C     PARAMETERS

      logical extrp2,extrp4,extrp5
      parameter ( extrp2   = .false. )
      parameter ( extrp4   =  .true. )
      parameter ( extrp5   =  .true. )

      double precision eta,phieps,mseps,rhoms,sigma,delini,next,beta,
     +        gamma,delmin,alpha,lspgmi,lspgma,sigma1,sigma2

      parameter ( phieps   = 1.0d-05 )
      parameter ( mseps    = 1.0d-05 )

      parameter ( delini   = 1.0d+02 )
      parameter ( delmin   = 1.0d-04 )
      parameter ( sigma1   = 1.0d-01 )
      parameter ( sigma2   = 5.0d-01 )
      parameter ( lspgmi   = 1.0d-10 )
      parameter ( lspgma   = 1.0d+10 )
      parameter ( eta      = 1.0d-01 )
      parameter ( alpha    = 1.0d-01 )
      parameter ( beta     = 5.0d-01 )
      parameter ( gamma    = 1.0d-04 )
      parameter ( sigma    = 2.0d-01 )
      parameter ( next     =   4.0d0 )
      parameter ( rhoms    = 9.0d-01 )

C     LOCAL SCALARS
      integer nitspgg,nitspgi,nefspgg,nefspgi

C     CALL BETRA

      call betra(n,x,l,u,m,lambda,rho,equatn,epsgpsn,maxit,maxfc,iprint,
     +ncomp,sterel,steabs,epsrel,epsabs,infrel,infabs,f,g,gpsupn,ldh,H,
     +iter,fcnt,gcnt,hcnt,chcnt,nitspgg,nitspgi,nefspgg,nefspgi,inform,
     +wi1,wi2,wi3,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8,eta,gamma,lspgmi,
     +lspgma,sigma1,sigma2,next,rhoms,mseps,phieps,sigma,delini,beta,
     +alpha,delmin,extrp2,extrp4,extrp5)

      end

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

C     Last update of BETRA or any of its dependencies: 
C
C     April 4th, 2005.
C
C     See report of modifications at the end of this file.

c     ------------------------------------------------------------------

c     Box-Euclidian-trust-region algorithm (BETRA): algorithm to find a 
c     minimizer of a box-constrained function.

c     The problem to be solved by this algorithm is:

c                            min  f(x)
c                            s.t. l <= x <= u,

c     with l and u in R^n and f(x) twice-continuously differentiable.

c     This is a active set algorithm. The box of constraints is divided 
c     into open faces and, at each iteration k, the current point x^k 
c     belongs to a open face F.

c     The first step of iteration k is to decide if, in the current 
c     iteration, BETRA will continue to minimize function f at face F or 
c     if it will leave this face and go to a point x^{k+1} in another 
c     face F'.

c     This decision is taken based on the relation between the 
c     Euclidian-norm of function projected gradient at x^k (g_P(x^k)) 
c     and internal gradient at x^k (g_I(x^k)). 

c     When ||g_I(x^k)|| < eta * ||g_P(x^k)||, with eta in (0,1), the 
c     decision is to leave face F. The idea is that, when this relation 
c     is true, the gradient at x^k projected on the box of constraints 
c     (g_P(x^k)) says that there is a good chance of finding a smaller 
c     value of the function if the working face is changed. Is this 
c     case, a Spectral Projected Gradient algorithm iteration (SPG) is 
c     used to calculate x^{k+1} that belongs to a new face F'.

c     When ||g_I(x^k)|| >= eta * ||g_P(x^k)||, the decision is to 
c     continue working on the current face, because it seems that there 
c     is a point in the current face that has better (smaller) function 
c     value. In this case, an inner algorithm iteration (INNER) is 
c     performed to calculate x^{k+1} in face F. The inner algorithm uses 
c     SPG when x^k is close to the boundary and trust-region like 
c     algorithm to find x^{k+1} that gives sufficient function value 
c     decrease or hits the boundary with simple function value decrease 
c     (see comments of INNER).

c     The advantage of working on this faces is that, once the algorithm 
c     is in a face, only the free variables at the current face are used, 
c     reducing the size of the problem.

c     The main algorithm stops when the inner algorithm stops declaring 
c     that a first- or second-order stationary point is found or when 
c     the current point is a vertex with null projected gradient.

c     OBS: the user must provide subroutines to evaluate function, 
c     gradient and Hessian. The prototypes of this subroutines are as 
c     follows:

c     ------------------------------------------------------------------

c     subroutine evalF(n, x, fx, flag)
      
c     integer n, flag
c     double precision x(n), fx

c     Evaluation of objective function at x. Function value stored at fx.
c     If some problem occurred, flag should be set to a negative value.

c     end

c     ------------------------------------------------------------------

c     subroutine evalG(n, x, g, flag)

c     integer n, flag
c     double precision x(n), g(n)

c     Evalutation of function gradient at x. Gradient stored at g.
c     If some problem occurred, flag should be set to a negative value.

c     end

c     ------------------------------------------------------------------

c     subroutine evalH(n, x, NMAX, H, diag, flag)

c     integer n, NMAX, flag
c     double precision x(n), H(NMAX,n), diag(n)

c     Evaluation of function Hessian matrix. Hessian matrix should be 
c     stored at H, and its diagonal must be stored at diag.
c     If some problem occurred, flag should be set to a negative value.

c     end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     On entry:
c     n        - integer
c                dimension of full space;

c     x        - double precision (n)
c                initial point;

c     l        - double precision (n)
c                lower bounds on x;

c     u        - double precision (n)
c                upper bounds on x;

c     NMAX     - integer
c                leading dimension of matrix H;

c     eta      - double precision
c                used to decide if the current face should be left; 
c                <recommended: 1.0d-1>

c     epsg     - double precision
c                allowed error for projected gradient norm;
c                <recommended: 1.0d-5>

c     rhoms      - double precision
c                distance from the boundary (to prevent that the error 
c                in the lenght of step calculated by More-Sorensen 
c                algorithm make the solution candidate lay outside the 
c                box of constraints) (INNER);
c                <recommended: 9.0d-1>

c     mseps    - double precision
c                tolerance for More-Sorensen algorithm (MEQB);
c                <recommended: 1.0d-5>

c     phieps   - double precision
c                allowed error in comparing the quadratic model value to 
c                zero (INNER); 
c                <recommended: 1.0d-5>

c     sigma    - double precision
c                allowed error to trust-region subproblem solution 
c                (MEQB);
c                <recommended: 2.0d-1>

c     deltain  - double precision
c                used to calculate the initial trust-region radius 
c                (INNER);
c                <recommended: 1.0d2>

c     ene      - double precision
c                size of extrapolation (INNER);
c                <recommended: 4.0d0>

c     beta     - double precision
c                used to decide whether or not extrapolation must be 
c                performed (INNER);
c                <recommended: 5.0d-1>

c     gamma    - double precision
c                sufficient decrease occurs when ared >= gamma*pred 
c                (INNER);
c                <recommended: 1.0d-1>

c     dmin     - double precision
c                minimum value for trust-region radius (delta) in the 
c                beggining of each iteration (INNER);
c                <recommended: 1.0d-4>

c     alpha    - double precision
c                Armijo condition constant (SPG);
c                <recommended: 1.0d-4>

c     sigmamin - double precision
c                lower bound on spectral step (SPG);
c                <recommended: 1.0d-10>

c     sigmamax - double precision
c                upper bound on spectral step (SPG);
c                <recommended: 1.0d10>

c     mu1      - double precision
c                safeguard for line search (lower) (SPG);
c                <recommended: 1.0d-1>

c     mu2      - double precision
c                safeguard for line search (upper) (SPG);
c                <recommended: 5.0d-1>

c     epsrel   - double precision
c                used to calculate the first spectral step (SPG);
c                <recommended: 1.0d-7>

c     epsabs   - double precision
c                used to calculate the first spectral step (SPG) and in 
c                the extrapolation procedure;
c                <recommended: 1.0d-10>

c     sterel   - double precision
c                not used;

c     steabs   - double precision
c                not used;

c     infrel   - double precision
c                not used;

c     infabs   - double precision
c                not used;

c     extrp2   - logical
c                true iff extrapolation is wanted after step 2 
c                (candidate solution calculated by an inner SPG 
c                iteration) (INNER);
c                <recommended: .false.>

c     extrp4   - logical
c                true iff extrapolation is wanted after step 4
c                (candidate solution is on the boundary);
c                iteration) (INNER);
c                <recommended: .true.>

c     extrp5   - logical
c                true iff extrapolation is wanted after step 5
c                (candidate solution provides sufficient decrease of 
c                objective function and is interior) (INNER);
c                <recommended: .true.>

c     ncomp    - integer
c                number of components to be printed;

c     maxit    - integer
c                maximum allowed number of iterations;

c     maxfc    - integer
c                maximum allowed number of function evaluations;

c     iprint   - integer
c                data printing options:
c                -1, nothing is printed;
c                 0, just the data from final iteration is printed;
c                >0, data obtained at each iprint iterations is printed;
                           
c     On return:
c     x        - double precision (n)
c                solution given by the algorithm;

c     g        - double precision (n)
c                gradient at x;

c     nigp     - double precision
c                infinite-norm of projected gradient at x;

c     fx       - double precision
c                objective function value at x;

c     nit      - integer
c                total number of iterations;

c     nitSPGi  - integer
c                number of inner SPG iterations;

c     nitSPGg  - integer
c                number of global SPG iterations;

c     nef      - integer
c                total number of function evaluations;

c     nefSPGi  - integer
c                number of function evaluations performed by inner SPG;

c     nefSPGi  - integer
c                number of function evaluations performed by global SPG;

c     nch      - integer
c                number of Cholesky decompositions;

c     status   - integer
c                0, if maximum allowed number of iterations is achieved;
c                1, if x is close to the boundary and is a first-order 
c                   stationary point (INNER);
c                2, if x is a second-order stationary point (INNER);
c                3, if x is a vertex with projected gradient null;
c                4, if decides to leave the current face with projected 
c                   gradient null;
c                5, if trust-region radius is too small (INNER);
c                6, if maximum allowed number of function evaluations is 
c                   achieved;
c              -90, if an error occurred during function evaluation;
c              -92, if an error occurred during gradient evaluation;
c              -94, if an error occurred during hessian evaluation;

c     Auxiliaries:
c     xant     - double precision (n)
c                point calculated in previous iteration;

c     free     - integer (nr)
c                indices of free variables;

c     gant     - double precision (n)
c                gradient at xant, used to calculate the spectral step 
c                sigma when SPG is called. Also used as auxiliary 
c                vector;

c     H       - double precision (NMAX,n)
c               Hessian matrix at certain point (stored at upper 
c               triangular part);

c     diag    - double precision (n)
c               Hessian matrix diagonal;

c     d       - double precision (n)
c               descent direction;

c     y       - double precision (n)
c               solution candidate;

c     tmp1    - double precision (n)
c               auxiliary vector;

c     tmp2    - double precision (n)
c               auxiliary vector;

c     tmp3    - double precision (n)
c               auxiliary vector;

c     rbdind  - integer (nr)
c               indices of x that hit the boundary;

c     rbdtype - integer (nr)
c               for each index i of x that hit the boundary, rbdtype(i) 
c               is set to
c               1, if component hit lower bound;
c               2, if component hit upper bound.
          
      subroutine BETRA(n, x, l, u, m, lambda, rho, equatn, epsg, maxit, 
     *maxfc, iprint, ncomp, sterel, steabs, epsrel, epsabs, infrel, 
     *infabs, fx, g, nigp, NMAX, H, nit, nef, neg, neh, nch, nitSPGg, 
     *nitSPGi, nefSPGg, nefSPGi, status, free, rbdind, rbdtype, xant, 
     *gant, diag, d, y, tmp1, tmp2, tmp3, eta, alpha, sigmamin, 
     *sigmamax, mu1, mu2, ene, rhoms, mseps, phieps, sigma, deltain, 
     *beta, gamma, dmin, extrp2, extrp4, extrp5)
                      
      integer m, n, maxit, maxfc, iprint, ncomp, NMAX, nit, nef, neg, 
     *neh, nch, nitSPGg, nitSPGi, nefSPGg, nefSPGi, status, free(n), 
     *rbdind(n), rbdtype(n)
      double precision x(n), l(n), u(n), lambda(m), rho(m), epsg, 
     *sterel, steabs, epsrel, epsabs, infrel, infabs, fx, g(n), nigp, 
     *H(NMAX,n), xant(n), gant(n), diag(n), d(n), y(n), tmp1(n), 
     *tmp2(n), tmp3(n), eta, alpha, sigmamin, sigmamax, mu1, mu2, ene,
     *rhoms, mseps, phieps, sigma, deltain, beta, gamma, dmin
      logical extrp2, extrp4, extrp5
      logical equatn(m)

c     Local variables:
c     nr      - integer
c               dimension of the reduced space, that is, the number of 
c               free variables at each iteration;

c     neg     - integer
c               number of gradient evaluations;

c     neh     - integer
c               number of Hessian evaluations;

c     nextrb  - integer
c               number of well-succeeded extrapolations (INNER);

c     nextrm  - integer
c               number of not well-succeeded extrapolations (INNER);

c     i       - integer
c               counter;

c     nprint  - integer
c               number of components to actually printed;

c     lamms  - double precision
c               value that More-Sorensen algorithm calculates to find 
c               the trust-region subproblem solution (MEQB);

c     delta   - double precision
c               trust-region radius (INNER);

c     nigi    - double precision
c               infinite-norm of internal gradient (at x);

c     nx      - double precision
c               Euclidian-norm of x;

c     fy      - double precision
c               function evaluated at y;

c     t       - double precision
c               scalar used to calculate y = x + t*d (SPG);

c     ittype  - character *4
c               type of iteration used to calculate y:
c                 TR, for trust-region;
c               iSPG, for inner SPG;
c               gSPG, for global SPG;
c               frst, if y is the initial point;

c     calculG - logical
c               true iff gradient at x must be calculated at the 
c               begining of next iteration.

      integer nr, nextrb(3), nextrm(3), i, nprint
      double precision lamms, delta, nigi, nx, fy, t
      character*4 ittype
      logical calculG

c     Initialize values.
      
      lamms    = 0.0d0
      status    = -1
      nit       = 0
      nitSPGi   = 0
      nitSPGg   = 0
      nef       = 0
      nefSPGi   = 0
      nefSPGg   = 0
      neg       = 0
      neh       = 0
      nch       = 0
      nextrb(1) = 0
      nextrb(2) = 0
      nextrb(3) = 0
      nextrm(1) = 0
      nextrm(2) = 0
      nextrm(3) = 0
      nprint    = min(ncomp, n)
      ittype    = 'frst'

      calculG = .true.

c     Initial point must lay inside the box.

      call project(n, x, l, u)

c     Calculate x Euclidian-norm (used to calculate the initial 
c     trust-region radius).

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

c     Calculate the initial trust-region radius. 

      delta = dmax1(dmin, deltain * dmax1(1.0d0, nx))
      
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     Evaluate function at x (entry parameter for either SPG or inner 
c     algorithm).

      call evalal(n,x,m,lambda,rho,equatn,fx,status)
      if (status .lt. 0) then
         if ( iprint .ge. 2 ) then
            write(*,994) status
            write(10,994) status
         end if
         return
      end if

      nef = nef + 1

c     Verify in which face the point x lays, that is, verifies which 
c     variables are free.

 1    nr = 0

      do i = 1, n
         if ((x(i) .gt. l(i)) .and. (x(i) .lt. u(i))) then
            nr = nr + 1
            free(nr) = i
         end if
      end do
      
c     Evaluate the function gradient at x.

      if (calculG) then
         call evalnal(n,x,m,lambda,rho,equatn,g,status)
         if (status .lt. 0) then
            if ( iprint .ge. 2 ) then
               write(*,994) status
               write(10,994) status
            end if
            return
         end if
         neg = neg + 1
      end if
      calculG = .true.

c     Calculate the projected gradient at x.

      do i = 1, n
         tmp1(i) = x(i) - g(i)
      end do

      call project(n, tmp1, l, u) 

      do i = 1, n
         tmp1(i) = tmp1(i) - x(i)
      end do

c     Calculate the projected and internal gradient infinite-norm.

      nigp = 0.0d0
      nigi = 0.0d0
      do i = 1, n
         nigp = dmax1(nigp, dabs(tmp1(i)))
         if ((x(i) .gt. l(i)) .and. (x(i) .lt. u(i))) then
            nigi = dmax1(nigi, dabs(tmp1(i)))
         end if
      end do

c     Print problem informations (current iteration).

      if ( iprint .ge. 2 ) then
         write(*, 982) nit,ittype
         write(*, 983) nprint,(x(i),   i=1,nprint)
         write(*, 984) nprint,(g(i),   i=1,nprint)
         write(*, 985) nprint,(tmp1(i),i=1,nprint)
         write(*, 986) min0(nprint,nr),nr,(free(i),i=1,min(nprint,nr))
         write(*, 995) fx,nigp,nigi,nr,n,nitSPGi,nitSPGg,
     +                 nit - nitSPGi - nitSPGg,nef,neg,neh,nch

         write(10,982) nit,ittype
         write(10,983) nprint,(x(i),   i=1,nprint)
         write(10,984) nprint,(g(i),   i=1,nprint)
         write(10,985) nprint,(tmp1(i),i=1,nprint)
         write(10,986) min0(nprint,nr),nr,(free(i),i=1,min(nprint,nr))
         write(10,995) fx,nigp,nigi,nr,n,nitSPGi,nitSPGg,
     +                 nit - nitSPGi - nitSPGg,nef,neg,neh,nch
      end if

c     If the current point is a vertex whose projected gradient norm is 
c     null, the algorithm stops with this vertex as a solution.

      if ((nr .eq. 0) .and. (nigp .le. epsg)) then
         status = 3

         if ( iprint .ge. 2 ) then
            write(*, 990) status,epsg
            write(10,990) status,epsg
         end if

         go to 999
      end if

c     Main test. If the condition is satisfied, the point x is reduced 
c     to the current face, as well as its gradient and bounds. Then, an 
c     iteration of the inner algorithm is performed and everything is 
c     restored to the original space (dimension n). If the condition is 
c     not satisfied, a iteration of SPG is performed to decide whether 
c     or not it should continue working on the current face.

      if (nigi .ge. eta*nigp) then

         call shrink(n, nr, free, x)
         call shrink(n, nr, free, xant)
         call shrink(n, nr, free, g)
         call shrink(n, nr, free, gant)
         call shrink(n, nr, free, l)
         call shrink(n, nr, free, u)

         call INNER(n, nr, x, m, lambda, rho, equatn,
     *   xant, free, l, u, g, gant, NMAX, H, diag,
     *   d, y, tmp1, tmp2, tmp3, fx, fy, rbdind, rbdtype, delta, lamms,
     *   rhoms, mseps, phieps, sigma, ene, beta, gamma, dmin, epsg, 
     *   alpha, sigmamin, sigmamax, mu1, mu2, epsrel, epsabs, 
     *   nit, nitSPGi, nef, nefSPGi, neg, neh, nch, nextrb, nextrm, 
     *   ittype, status, calculG, extrp2, extrp4, extrp5, maxfc)

         if (status .lt. 0) then
            if ( iprint .ge. 2 ) then
               write(*,994) status
               write(10,994) status
            end if
            return
         end if

c        OBS: because the function has already been evaluated at y, its 
c             upper portion already contains the fixed elements of x. So, 
c             y can be expanded.

         do i = nr+1, n
            y(i) = x(i)
         end do

         call expand(n, nr, free, x)   
         call expand(n, nr, free, y)
         call expand(n, nr, free, g)
         call expand(n, nr, free, l)
         call expand(n, nr, free, u)

c        If the inner algorithm stops declaring "first-order stationary 
c        point close to the boundary" (status = 1) or "second-order 
c        stationary point" (status = 2) then the main algorithm also 
c        stops (when the projected gradient is null).

         if ((status .eq. 1) .and. (nigp .le. epsg)) then

            if ( iprint .ge. 2 ) then
               write(*, 988) status,epsg,2*dmin
               write(10,988) status,epsg,2*dmin
            end if

            go to 999
         end if

         if ((status .eq. 2) .and. (nigp .le. epsg)) then

            if ( iprint .ge. 2 ) then
               write(*, 989) status 
               write(10,989) status
            end if

            go to 999
         end if

c        If the inner algorithm stops because the trust-region radius 
c        became too small without obtaining sufficient decrease 
c        (status = 5) then the main algorithm also stops.

         if (status .eq. 5) then
            if ( iprint .ge. 2 ) then
               write(*, 992) status 
               write(10,992) status
            end if
            go to 999
         end if

c        If the inner algorithm stops because maximum number of function 
c        evaluations is achieved, the main algorithm also stops.

         if (status .eq. 6) then
            if ( iprint .ge. 2 ) then
               write(*,993) status
               write(10,993) status
            end if
            go to 999
         end if

c        There are cases when the inner algorithm stops declaring 
c        "first-order stationary point close to the boundary" 
c        (status = 1) or "second-order stationary point" (status = 2), 
c        but the projected gradient is not null. This happens when eta 
c        is too small. In such cases, the main algorithm continues to 
c        work in the current face (using the inner algorithm) when it 
c        would be better to leave this face (by using global SPG).
c        Because of this, when the inner algorithm stops with 
c        status = 1 or status = 2 but the projected gradient is not
c        null, an iteration of the global SPG is performed.

         if ((status .eq. 1) .or. (status .eq. 2)) then
            status = -1

            call expand(n, nr, free, xant)
            call expand(n, nr, free, gant)

            call SPG(n, n, x, m, lambda, rho, equatn,
     *      xant, free, l, u, g, gant, d, y, t, fx,
     *      fy, alpha, sigmamin, sigmamax, mu1, mu2, epsrel, epsabs,  
     *      nit, nef, nefSPGg, neg, ittype, .false., maxfc, status)

            if (status .lt. 0) then
               if ( iprint .ge. 2 ) then
                  write(*,994) status
                  write(10,994) status
               end if
               return
            end if

c           If the SPG iteration stops because maximum number of function 
c           evaluations is achieved, the main algorithm also stops.

            if (status .eq. 6) then
               if ( iprint .ge. 2 ) then
                  write(*,993) status
                  write(10,993) status
               end if
               go to 999
            end if
            
            nitSPGg = nitSPGg + 1
         end if
         
      else
         if (nigp .le. epsg) then
            status = 4

            if ( iprint .ge. 2 ) then
               write(*, 991) status,epsg 
               write(10,991) status,epsg
            end if

            go to 999
         end if

         call SPG(n, n, x, m, lambda, rho, equatn, 
     *   xant, free, l, u, g, gant, d, y, t, fx, fy,
     *   alpha, sigmamin, sigmamax, mu1, mu2, epsrel, epsabs, nit, nef, 
     *   nefSPGg, neg, ittype, .false., maxfc, status)

         if (status .lt. 0) then
            if ( iprint .ge. 2 ) then
               write(*,994) status
               write(10,994) status
            end if
            return
         end if

c        If the SPG iteration stops because maximum number of function 
c        evaluations is achieved, the main algorithm also stops.

         if (status .eq. 6) then
            if ( iprint .ge. 2 ) then
               write(*,993) status
               write(10,993) status
            end if
            go to 999
         end if

         nitSPGg = nitSPGg + 1
      end if

c     Updates xant, gant and fx, just in case the next iteration is a 
c     global or inner SPG iteration.
      
      do i = 1, n
         xant(i) = x(i)
         x(i) = y(i)
      end do

c     When the gradient at y has already been calculated by the inner 
c     algorithm (in this case, g already contains the gradient at the 
c     point to be set as the current one), the value of gant has also 
c     been updated  (although the latter must be rearranged).

      if (calculG) then
         do i = 1, n
            gant(i) = g(i)
         end do
      else
         call expand(n, nr, free, gant)         
      end if

      fx = fy
      
c     Start new iteration.

      if (nit .ge. maxit) then
         status = 0

         if ( iprint .ge. 2 ) then
            write(*, 987) status
            write(10,987) status
         end if

         go to 999
      end if

      nit = nit + 1
      go to 1

c     Print final informations.

 999  if ( iprint .ge.  2 ) then
         write(*, 981) nit,ittype         
         write(*, 983) nprint,(x(i),   i=1,nprint)
         write(*, 984) nprint,(g(i),   i=1,nprint)
         write(*, 985) nprint,(tmp1(i),i=1,nprint)
         write(*, 986) min(nprint,nr),nr,(free(i),i=1,min(nprint,nr))
         write(*, 995) fx,nigp,nigi,nr,n,nitSPGi,nitSPGg, 
     *                 nit - nitSPGi - nitSPGg,nef,neg,neh,nch

         write(10,981) nit,ittype         
         write(10,983) nprint,(x(i),   i=1,nprint)
         write(10,984) nprint,(g(i),   i=1,nprint)
         write(10,985) nprint,(tmp1(i),i=1,nprint)
         write(10,986) min(nprint,nr),nr,(free(i),i=1,min(nprint,nr))
         write(10,995) fx,nigp,nigi,nr,n,nitSPGi,nitSPGg, 
     *                 nit - nitSPGi - nitSPGg,nef,neg,neh,nch
      end if

c     Printing formats.

 977  format(/,1X,'Entry to BETRA. 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,'BETRA iteration: ', I6, 
     +            ' (Final point, obtained by a ',A4,' iteration.)')
 982  format(/,1X,'BETRA iteration: ', I6,
     +            ' (This point was obtained by a ',A4,' iteration.)')
 983  format(  1X,'Current point (first ', I6, ' components): ',
     +       /,6(1X,1PD11.4))
 984  format(  1X,'Current gradient (first ', I6, ' components): ',
     +       /,6(1X,1PD11.4))
 985  format(  1X,'Current continuous projected gradient (first ',I6, 
     +            ' components): ',
     +       /,6(1X,1PD11.4))
 986  format(  1X,'Current free variables (first ',I6, ', total number',
     +             I6,'): ',/,10(1X,I6))
 987  format(/,1X,'Status of BETRA = ', I2,
     +       /,1X,'(Number of iterations exceeded the maximum allowed)')
 988  format(/,1X,'Status of BETRA = ',I2,
     +       /,1X,'(Convergence with internal gradient Euclidian-norm ',
     +       /,1X,'less than or equal to ',1PD11.4,' and ', 
     +       /,1X,'distance to boundary less than or equal to ',1PD11.4,
     +            ')')
 989  format(/,1X,'Status of BETRA = ', I2,
     +       /,1X,'(Convergence to a second-order stationary point)')
 990  format(/,1X,'Status of BETRA = ',I2,
     +       /,1X,'(Convergence with projected gradient sup-norm less',
     +       /,1X,'than or equal to ', 1PD11.4,'at a vertex)')
 991  format(/,1X,'Status of BETRA = ',I2,
     +       /,1X,'(Convergence with projected gradient sup-norm less '
     +       /,1X,'than or equal to ', 1PD11.4,')')
 992  format(/,1X,'Status of BETRA = ',I2,
     +       /,1X,'(Doubtful convergence, trust-region radius too',
     +         1X,'small)')
 993  format(/,1X,'Status of BETRA = ', I2,
     +       /,1X,'(Number of function evaluations exceeded the',
     +         1X,'maximum allowed)')
 994  format(/,1X,'Status of BETRA = ', I2,
     +       /,1X,'(Error in function evaluation, if status = -90;',
     +       /,1X,'(Error in gradient evaluation, if status = -92;',
     +       /,1X,'(Error in hessian evaluation,  if status = -94)')
 995  format(  1X,'Functional value: ', 1PD11.4,
     +       /,1X,'Sup-norm of the projected gradient: ', 1PD11.4,
     +       /,1X,'Sup-norm of the internal projected gradient: ', 
     +             1PD11.4,
     +       /,1X,'Free variables at the current point: ', I7, 
     +            ' (over a total of ',I7,')',
     +       /,1X,'Inner SPG iterations: ', I7,
     +       /,1X,'Global SPG iterations: ', I7,
     +       /,1X,'Trus-region iterations: ', I7,
     +       /,1X,'Function evaluations: ', I7,
     +       /,1X,'Gradient evaluations: ', I7,
     +       /,1X,'Hessian evaluations: ', I7,
     +       /,1X,'Cholesky decompositions: ', I7)

      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Calculate the projection of x onto the box bounded by l and u.

c     On entry:
c     n - integer
c         dimension of x;

c     x - double precision (n)
c         vector in R^n to be projected;

c     l - double precision (n)
c         lower bounds;

c     u - double precision (n)
c         upper bounds;

c     On return:
c     x - double precision (n)
c         projected vector.

      subroutine project(n, x, l, u)

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

      integer i

      do i = 1, n
         if (x(i) .lt. l(i)) then
            x(i) = l(i)
         else 
            if (x(i) .gt. u(i)) then
               x(i) = u(i)
            end if
         end if
      end do

      end

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

c     Last update of INNER (or any of its dependencies): 
c     April 4th, 2005.
c     See report of modifications at the end of this file.

c     ------------------------------------------------------------------

c     Algorithm that finds a unconstrained minimizer of objective 
c     function inside the box of constraints, hits the boundary 
c     (obtaining function decrease), or finds an interior point where 
c     the objective function has sufficient decrease (compared to its 
c     value at x). Extrapolation may be done.

c     When the current point x is "close to" the boundary, a Spectral 
c     Projected Gradient (SPG) iteration is used to calculate the new 
c     point. If this new point is at the boundary, the algorithm stops. 
c     Otherwise, a new iteration begins.

c     When x is "far from" the boundary, trust-region radius is 
c     determined and d is calculated using More-Sorensen algorithm to 
c     solve the trust-region subproblem (which is to find a minimizer a 
c     to a function quadratic model provided that the minimizer's 
c     Euclidian-norm is smaller than a given delta). The new point is
c     y = x + d. 

c     If y lays outside the box of constraints, it is truncated on the 
c     boundary. This new y on the boundary will be candidate to be a 
c     solution. If function value at new y is smaller than function 
c     value at x, inner algorithm stops with y. Otherwise, the 
c     trust-region radius is decreased so that the new solution d' to 
c     the trust-region subproblem makes x + d' be interior to the box.
c     More-Sorensen algorithm is used to calculate d' too.

c     If y lays inside the box, sufficient decrease of objective 
c     function is tested. If it is true, y is accepted as a solution
c     candidate. If y in on the boundary, inner algorithm stops and if 
c     it is interior, a new iteration begins. If sufficient decrease is 
c     not obtained, trust-region radius is decreased and a new quadratic 
c     model minimizer is calculated (as in a classical trust-region 
c     algorithm for unconstrained minimization).

c     If the user wants, after calculating the candidate solution y, 
c     extrapolation may be performed. For this, set extrpi to true, 
c     where i is the step of the inner algorithm that can call 
c     extrapolation procedure. If extrp2 is true, extrapolation will be 
c     tried after y is calculated by an inner SPG iteration. If extrp4 
c     is true, extrapolation will be tried after y hit the boundary. And 
c     if extrp5 is true, extrapolation will be tried after y is 
c     calculated by trust-region algorithm, when y is interior and 
c     provides sufficient decrease of objective function.

c     If gradient at current point is null, inner algorithm stops 
c     declaring "first-order stationary point". If quadratic model 
c     minimum is 0, inner algorithm stops declaring "second-order 
c     stationary point".

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     On entry:
c     n        - integer
c                dimension of full space;

c     nr       - integer
c                dimension of reduced space;

c     x        - double precision (n)
c                initial point, interior to the current face;

c     xant     - double precision (n)
c                candidate solution calculated by previous iteration. It 
c                is used to calculate the spectral step sigma when 
c                subroutine SPG is called. This vector is altered in 
c                this subroutine;

c     free     - integer (nr)
c                indices of free variables;

c     l        - double precision (n)
c                lower bounds on x;

c     u        - double precision (n)
c                upper bounds on x;

c     g        - double precision (n)
c                gradient at x;

c     gant     - double precision (n)
c                gradient at xant, used to calculate the spectral step 
c                sigma when SPG is called. Also used as auxiliary 
c                vector;

c     NMAX     - integer
c                leading dimension of matrix H;

c     fx       - double precision
c                objective function value at x;

c     delta    - double precision 
c                trust-region radius;

c     lamms   - double precision
c                value that More-Sorensen algorithm calculates to find 
c                the trust-region subproblem solution (MEQB);

c     rhoms      - double precision
c                distance from the boundary (to prevent that the error 
c                in the lenght of step calculated by More-Sorensen 
c                algorithm make the solution candidate lay outside the 
c                box of constraints);

c     mseps    - double precision
c                tolerance for More-Sorensen algorithm (MEQB);

c     phieps   - double precision
c                allowed error in comparing the quadratic model value to 
c                zero; 

c     sigma    - double precision
c                allowed error to trust-region subproblem solution 
c                (MEQB);

c     ene      - double precision
c                size of extrapolation;

c     beta     - double precision
c                used to decide whether or not extrapolation must be 
c                performed;

c     gamma    - double precision
c                sufficient decrease occurs when ared >= gamma*pred;

c     dmin     - double precision
c                minimum value for trust-region radius (delta) in the 
c                beggining of each iteration;

c     epsg     - double precision
c                allowed error for projected gradient norm;

c     alpha    - double precision
c                Armijo condition constant (SPG);

c     sigmamin - double precision
c                lower bound on spectral step (SPG);

c     sigmamax - double precision
c                upper bound on spectral step (SPG);

c     mu1      - double precision
c                safeguard for line search (lower) (SPG);

c     mu2      - double precision
c                safeguard for line search (upper) (SPG);

c     epsrel   - double precision
c                used to calculate the first spectral step (SPG);

c     epsabs   - double precision
c                used to calculate the first spectral step (SPG) and in 
c                the extrapolation procedure;

c     nit      - integer
c                number of current interation;

c     extrp2   - logical
c                true iff extrapolation is wanted after step 2 
c                (candidate solution calculated by an inner SPG 
c                iteration);

c     extrp4   - logical
c                true iff extrapolation is wanted after step 4
c                (candidate solution is on the boundary);
c                iteration);

c     extrp5   - logical
c                true iff extrapolation is wanted after step 5
c                (candidate solution provides sufficient decrease of 
c                objective function and is interior);

c     maxfc    - integer
c                maximum allowed number of function evaluations;

c     On return:
c     x        - double precision (n)
c                solution candidate, with status as described bellow;

c     fy       - double precision
c                function value at y (solution candidate);

c     nitSPGi  - integer
c                number of inner SPG iterations;

c     nef      - integer
c                total number of function evaluations;

c     nefSPGi  - integer
c                number of function evaluations performed by inner SPG;

c     neg      - integer
c                total number of gradient evaluations;

c     neh      - integer
c                total number of Hessian evaluations;

c     nch      - integer
c                number of Cholesky decompositions;

c     nextrb   - integer
c                number of well-succeeded extrapolations;

c     nextrm   - integer
c                number of not well-succeeded extrapolations;

c     delta    - double precision
c                updated trut-region radius;

c     lamms   - double precision
c                updated value for next iteration (see entry paramater);

c     ittype   - character *4
c                iteration type: 
c                  TR, for trust-region;
c                iSPG, for inner SPG;

c     status   - integer
c                1, if x is close to the boundary and is a first-order 
c                   stationary point;
c                2, if x is a second-order stationary point;
c                3, if x hit the boundary;
c                5, if trust-region radius is too small;
c                6, if maximum allowed number of function evaluations is 
c                   achieved;
c              -90, if an error occurred during function evaluation;
c              -92, if an error occurred during gradient evaluation;
c              -94, if an error occurred during hessian evaluation;

c     calculG  - logical
c                false iff the point calculated already has its gradient 
c                calculated;

c     Auxiliaries:
c     xant    - double precision (n)
c               point calculated in previous iteration;

c     gant    - double precision (n)
c               gradient at xant;

c     H       - double precision (NMAX,n)
c               Hessian matrix at certain point (stored at upper 
c               triangular part);

c     diag    - double precision (n)
c               Hessian matrix diagonal;

c     d       - double precision (n)
c               descent direction;

c     y       - double precision (n)
c               solution candidate;

c     tmp1    - double precision (n)
c               auxiliary vector;

c     tmp2    - double precision (n)
c               auxiliary vector;

c     tmp3    - double precision (n)
c               auxiliary vector;

c     rbdind  - integer (nr)
c               indices of x that hit the boundary;

c     rbdtype - integer (nr)
c               for each index i of x that hit the boundary, rbdtype(i) 
c               is set to
c               1, if component hit lower bound;
c               2, if component hit upper bound;

      subroutine INNER(n, nr, x, m, lambda, rho, equatn,
     *xant, free, l, u, g, gant, NMAX, H,
     *diag, d, y, tmp1, tmp2, tmp3, fx, fy, rbdind, rbdtype, delta, 
     *lamms, rhoms, mseps, phieps, sigma, ene, beta, gamma, dmin, epsg, 
     *alpha, sigmamin, sigmamax, mu1, mu2, epsrel, epsabs,  nit, 
     *nitSPGi, nef, nefSPGi, neg, neh, nch, nextrb, nextrm, ittype, 
     *status, calculG, extrp2, extrp4, extrp5, maxfc)
                                                
      integer m, n, nr, free(nr), NMAX, rbdind(nr), rbdtype(nr), nit,
     *nitSPGi, nef, nefSPGi, neg, neh, nch, nextrb(3), nextrm(3), 
     *status, maxfc
      double precision lambda(m), rho(m),
     *x(n), xant(n), l(nr), u(nr), g(n), gant(n), H(NMAX,n), diag(n), 
     *d(nr), y(n), tmp1(n), tmp2(n), tmp3(n), fx, fy, 
     *delta, lamms, rhoms, mseps, phieps, sigma, ene, beta, gamma, dmin, 
     *epsg, alpha, sigmamin, sigmamax, mu1, mu2, epsrel, epsabs 
      character*4 ittype
      logical calculG, extrp2, extrp4, extrp5
      logical equatn(m)

c     Local variables:
c     i        - integer
c                counter;

c     nrbd     - integer
c                number of x components that hit the boundary;

c     type     - integer
c                1, if a x component hit the lower bound and
c                2, if a x component hit the upper bound;

c     extr     - integer
c                1, if extrapolation is requested at step 2;
c                2, if extrapolation is requested at step 4;
c                3, if extrapolation is requested at step 5;

c     convCrit - integer
c                store which convergence criterion was satisfied by 
c                MEQB. If both gradient and Hessian are null, cinvCrit 
c                is set to 0;

c     dbound   - double precision
c                biggest radius for a ball centered in x to stay 
c                interior to the box os constraints;

c     tmax     - double precision
c                biggest scalar such that x + tmax*d in interior to the 
c                box of constraints;

c     pred     - double precision
c                expected reduction of objective function (difference 
c                between value of quadratic model at 0 and y);

c     ared     - double precision
c                actual objective function reduction (fx - fy); 

c     nd       - double precision
c                Euclidian-norm of step d calculated by More-Sorensen 
c                algorithm;

c     nig      - double precision
c                gradient infinite-norm;

c     temp     - double precision
c                temporary variable;

c     phi      - double precision
c                quadratic model evaluated at a certain point;

c     dTg      - double precision
c                vectors d and g dot product;

c     dTgt     - double precision
c                dot product of vector d and gradient at x + g;

c     niy      - double precision
c                infinite-norm of y (= x + t*d);

c     nidy     - double precision
c                infinite-norm of y - (x + tt*d);

c     fyt      - double precision
c                function evaluated at x + tt*d;

c     a        - double precision
c                scalar that multiplies d when y is calculated;

c     aa       - double precision
c                scalar to multiply vector d;

c     amax     - double precision
c                biggest positive scalar such that x + amax*d is still 
c                interior to the box of constraints;

c     SPGiter  - logical
c                true iff a SPG iteration is performed.

      integer i, itemp, nrbd, type, extr, convCrit
      double precision dbound, tmax, pred, ared, nd, nig, temp, phi, 
     *dTg, dTgt, niy, nidy, fyt, a, aa, amax
      logical SPGiter

c     step 1: calculate the distance between x and the boundary.
c             dbound is the largest positive delta such that the ball 
c             centered at x with radius delta is still inside the box 
c             (set of constraints).

      dbound = x(1) - l(1)

      do i = 2, nr
         dbound = dmin1(dbound, x(i) - l(i))
      end do

      do i = 1, nr
         dbound = dmin1(dbound, u(i) - x(i))
      end do
      
c     Calculate infinite-norm of gradient.

      nig = 0.0d0
      do i = 1, nr
         nig = dmax1(nig, dabs(g(i)))
      end do
            
      SPGiter = .false.

c      step 2: close to the boundary: if the gradient is null then 
c              algorithm stops. Else, performe an inner SPG iteration.

      if (dbound .lt. 2*dmin) then
         if (nig .le. epsg) then
            status = 1
            return
         else
            call SPG(n, nr, x, m, lambda, rho, equatn,
     *      xant, free, l, u, g, gant, d, y, tmax,
     *      fx, fy, alpha, sigmamin, sigmamax, mu1, mu2, epsrel, epsabs,  
     *      nit, nef, nefSPGi, neg, ittype, .true., maxfc, status)

            if ((status .lt. 0) .or. (status .eq. 6)) then
               return
            end if

            nitSPGi = nitSPGi + 1
            SPGiter = .true.

            if (.not. extrp2) then
               go to 8
            end if

            extr = 1
            go to 6
         end if
      end if

      ittype = 'TR'

c     step 3: far from the boundary, solve trust-region subproblem. 
c             Evaluate function Hessian at x.

      call calcH(nr,free,x,n,x,m,lambda,rho,equatn,NMAX,H,diag,status)
      if (status .lt. 0) then
         return
      end if

      neh = neh + 1

c     step 4: solve the trust-region subproblem using More-Sorensen's 
c             algorithm to minimize "exactly" quadratics subjected to 
c             balls. 

c     If trust-region radius is too small, the inner algorithm stops.

 4    if (delta .lt. epsabs) then
         status = 5
         return
      end if

      call MEQB(nr, g, NMAX, H, diag, d, tmp1, tmp2, tmp3, lamms,
     *delta, sigma, 0.0d0, mseps, -1, nch, convCrit)
                        
c     If both internal gradient and Hessian matrix are null, subroutine 
c     MEQB stops with convCrit = 0 and then the inner algorithm stops 
c     declaring "second-order stationary point".

      if (convCrit .eq. 0) then
         status = 2
         return
      end if

c     Evaluate the quadratic model of the objective function at d.
      
      call quad(nr, d, g, NMAX, H, diag, tmp1, phi)
      
c     If the value of the quadratic model at d is 0 it means that x is a 
c     second-order stationary point. In this case, inner algorithm stops
c     declaring this.
      
      if ((dabs(phi) .le. phieps) .and. (nig .le. epsg)) then
         status = 2
         return
      end if

c     Calculate point y = x + d.

      do i = 1, nr
         y(i) = x(i) + d(i)
      end do

c     Verify if y is inside de box. If not, interior is set to false and 
c     tmax is set to the biggest positive scalar such that x + tmax*d is 
c     inside the box.
 
      tmax = 1.0d0
      nrbd = 0

      do i = 1, nr
         if (dabs(d(i)) .ge. dmax1(epsrel * dabs(x(i)), epsabs)) then
            if (y(i) .le. l(i)) then
               temp = (l(i) - x(i)) / d(i)
               if (temp .lt. tmax) then
                  tmax = temp
                  nrbd = 1
                  rbdind(nrbd) = i
                  rbdtype(nrbd) = 1
               else 
                  if (temp .le. tmax + dmax1(epsrel*tmax, epsabs)) then
                     nrbd = nrbd + 1
                     rbdind(nrbd) = i
                     rbdtype(nrbd) = 1
                  end if
               end if
            else 
               if (y(i) .ge. u(i)) then
                  temp = (u(i) - x(i)) / d(i)
                  if (temp .lt. tmax) then
                     tmax = temp
                     nrbd = 1
                     rbdind(nrbd) = i
                     rbdtype(nrbd) = 2
                  else
                     if (temp .le. tmax+dmax1(epsrel*tmax, epsabs)) then
                        nrbd = nrbd + 1
                        rbdind(nrbd) = i
                        rbdtype(nrbd) = 2
                     end if
                  end if
               end if
            end if
         end if
      end do
      
c     If y is not interior to the box, y = x + d is replaced by 
c     y = x + tmax*d. Now y is definitly interior. Actually, it is in 
c     the boundary. If the objective function decreases in y, the inner 
c     algorithm stops with y as a solution. Otherwise, a new 
c     trust-region radius delta is chosen (smaller than dbound) and a 
c     new quadratic model minimizer is calculated (which is necessarilly 
c     interior because of the choise of delta).

      if (nrbd .gt. 0) then
         do i = 1, nr
            y(i) = x(i) + tmax * d(i)
         end do
         
c     Set x(i) to l(i) or u(i) for the indices i that got to the 
c     boundary (to prevent errors).

         do i = 1, nrbd
            itemp = rbdind(i)
            if (rbdtype(i) .eq. 1) then
               y(itemp) = l(itemp)
            else
               y(itemp) = u(itemp)
            end if
         end do

         call calcF(nr,free,y,n,x,m,lambda,rho,equatn,fy,status)
         if (status .lt. 0) then
            return
         end if

         nef = nef + 1
         if (nef .ge. maxfc) then
            status = 6
            return
         end if
         
         if(fy .lt. fx) then
            status = 3

            if (.not. extrp4) then
               go to 8
            end if

            extr = 2
            go to 6
         else
            temp = dmin + rhoms * ((dbound / (1 + sigma)) - dmin)
            delta = dmax1(dmin, temp)
            go to 4
         end if
      end if

c     step 5: in this case y is inside the box. Acceptance or rejection 
c             of the trust-region subproblem solution.

      pred = -phi

      call calcF(nr,free,y,n,x,m,lambda,rho,equatn,fy,status)
      if (status .lt. 0) then
         return
      end if

      nef = nef + 1
      if (nef .ge. maxfc) then
         status = 6
         return
      end if

      ared = fx - fy
      
c     If there is not sufficient decrease of the function, the 
c     trust-region radius is decreased and the new quadratic model 
c     minimizer will be calculated. 

      nd = 0.0d0
      do i = 1, nr
         nd = nd + d(i)**2
      end do
      nd = dsqrt(nd)
      
      if (((pred .gt. phieps) .and. (ared .lt. gamma*pred)) .or. 
     *   ((pred .le. phieps) .and. (fx .lt. fy))) then
         delta = nd/4
         go to 4
      end if

c     If extrapolation at step 5 is not to be performed, point y is 
c     accepted.

      if (.not. extrp5) then
         go to 8
      end if

      extr = 3

c     step 6: extrapolation. 
      
 6    call calcG(nr,free,y,n,x,m,lambda,rho,equatn,tmp2,status)
      if (status .lt. 0) then
         return
      end if

      neg = neg + 1
      
      calculG = .false.

      dTg = 0.0d0
      do i = 1, nr
         dTg = dTg + d(i) * g(i)
      end do

      dTgt = 0.0d0
      do i = 1, nr
         dTgt = dTgt + d(i) * tmp2(i)
      end do

c     Extrapolation should not be performed. Point y calculated by SPG 
c     or by the trust-region iteration will be used.

      if (dTgt .ge. beta*dTg) then
         nextrm(extr) = nextrm(extr) + 1
         go to 7
      end if
            
c     Point y was calculated as x + tmax*d. Even when the iteration 
c     performed to calculate y was a inner SPG one (when it is not 
c     verified which is the scalar tmax that makes x + tmax*d hit the 
c     boundary), tmax was set to the scalar used to calculate 
c     y = x + t*d.

      a = tmax
      
c     When y was not calculated by inner SPG and a = tmax is smaller 
c     than 1, the old y was outside the box and tmax made the new y hit 
c     the boundary. In this case, by the definition of amax shown below, 
c     amax is set to tmax.

      if ((a .lt. 1.0d0) .and. (.not. SPGiter)) then
         amax = tmax
         go to 61
      end if

c     Calculate amax, which is the biggest positive scalar such that 
c     x + amax*d is still in the box delimited by l and u.
      
      amax = 1.0d30
      nrbd = 0
            
      do i = 1, nr
         if (d(i) .lt. 0.0d0) then
            temp = (l(i) - x(i)) / d(i)
            type = 1
         else  
            if (d(i) .gt. 0.0d0) then
               temp = (u(i) - x(i)) / d(i)
               type = 2
            else
               temp = 1.0d30
               type = 3
            end if
         end if
         
         if (temp .lt. amax) then
            amax = temp
            nrbd = 1
            rbdind(nrbd) = i
            rbdtype(nrbd) = type
         else 
            if (temp .le. amax + dmax1(epsrel * amax, epsabs)) then
               nrbd = nrbd + 1
               rbdind(nrbd) = i
               rbdtype(nrbd) = type
            end if
         end if
      end do
      
c     Choose the scalar that will multiply d to get the new point.
      
 61   if ((a .lt. amax) .and. (ene*a .gt. amax)) then
         aa = amax
      else 
         aa = ene * a
      end if
      
c     There already is a point y that satisfies the conditions for the 
c     iteration termination. Now a new candidate point (tmp1) is calculated.

      do i = 1, nr
         tmp1(i) = x(i) + aa * d(i)
      end do
      
c     If aa is equal to amax, set x(i) to l(i) or u(i) for the indices i 
c     that hit the boundary (to avoid errors in future calculations).

      if (aa .ge. amax - dmax1(epsrel * amax, epsabs)) then
         do i = 1, nrbd
            itemp = rbdind(i)
            if (rbdtype(i) .eq. 1) then
               tmp1(itemp) = l(itemp)
            else
               tmp1(itemp) = u(itemp)
            end if
         end do
      end if

      call project(nr, tmp1, l, u)
      
c     Calculate infinite-norm of the point tmp1 - y.

      nidy = 0.0d0
      do i = 1, nr
         temp = dabs(tmp1(i) - y(i))
         nidy = dmax1(nidy, temp)
      end do
      
c     Calculate infinite-norm of y.

      niy = 0.0d0
      do i = 1, nr
         niy = dmax1(niy, dabs(y(i)))
      end do
      
      temp = dmax1(epsabs, epsrel * niy)
     
      if ((a .lt. amax) .or. (nidy .ge. temp)) then

c        Evaluate function at tmp1.
 
         call calcF(nr,free,tmp1,n,x,m,lambda,rho,equatn,fyt,status)
         if (status .lt. 0) then
            return
         end if

         nef = nef + 1
         if (nef .ge. maxfc) then
            status = 6
            return
         end if

c        If the objective function decreases from y to tmp1, tmp1 is 
c        accepted (and is called y, since it also satisfies a stopping 
c        criterion) and extrapolation is tried again. 
c        Note that, since tmp1 is accepted, the new gradient has not yet 
c        been calculated.

         if (fyt .lt. fy) then
            calculG = .true.
            a = aa
           
            do i = 1, nr
               y(i) = tmp1(i)
            end do
            
            fy = fyt
            nextrb(extr) = nextrb(extr) + 1

            go to 61
         else
            nextrm(extr) = nextrm(extr) + 1
         end if

      else
         nextrm(extr) = nextrm(extr) + 1
      end if
            
c     If extrapolation was not performed, y gradient has already been 
c     calculated (when the extrapolation was tried). In this case, the 
c     new gradient is in tmp2. Here, the valeus of g and gant are 
c     updated, so that g contains y gradient (new point) and gant 
c     constains x gradient (old point).

 7    if (.not. calculG) then
         do i = 1, n
            gant(i) = g(i) 
            g(i) = tmp2(i)
         end do
      end if

c     If new point y is in the boundary, the inner algorithm stops and 
c     returns y as solution.

 8    nrbd = 0

      do i = 1, nr
         if (y(i) .le. l(i) + dmax1(epsrel * dabs(l(i)), epsabs)) then
            y(i) = l(i)
            nrbd = nrbd + 1
            rbdind(nrbd) = i
            rbdtype(nrbd) = 1
         else if (y(i) .ge. u(i)-dmax1(epsrel*dabs(u(i)), epsabs)) then
            y(i) = u(i)
            nrbd = nrbd + 1
            rbdind(nrbd) = i
            rbdtype(nrbd) = 2
         end if
      end do
      
      if (nrbd .gt. 0) then
         status = 3
         return
      end if
      
c     Update the trust-region radius (which may or may not be used). 
c     This update can only be done when the current iteration was a 
c     trust-region iteration (and not an inner SPG one).

      if (.not. SPGiter) then
         if (ared .lt. pred/4) then
            delta = dmax1(nd/4, dmin)
         else
            if ((ared .ge. pred/2) .and. 
     *         (nd .ge. delta - dmax1(epsrel*delta, epsabs))) then
               delta = dmax1(2 * delta, dmin)
            else
               delta = dmax1(delta, dmin)
            end if
         end if
      end if

      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Reduct point x from space of dimension n to a space of dimension 
c     nr (space of free variables). To do so, in the first nr components 
c     of x are stored the values corresponding to the free variables 
c     and in the last n - nr components are stored the values of the 
c     fixed variables.

c     On entry:
c     n    - integer
c            dimension of full space;

c     nr   - integer
c            dimension of reduced space;

c     free - integer (nr)
c            free variables indices on original space;

c     x    - double precision (n)
c            point to be reduced to a space of smaller dimension 
c            (reduced space);

c     On return:
c     x    - double precision (n)
c            point to be reduced to a space of smaller dimension 
c            (reduced space). The first nr positions store the free 
c            variable values and the last n - nr store the fixed 
c            variable values.


      subroutine shrink(n, nr, free, x)

      integer n, nr, free(nr)
      double precision x(n)

      integer i, aux
      double precision temp

      do i = 1, nr
         aux = free(i)
         if (i .ne. aux) then
            temp = x(i)
            x(i) = x(aux)
            x(aux) = temp
         end if
      end do

      end 

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Expand point x from space of dimension nr to space of dimension n.
c     This is done by rearranging vector x so it is exactly the way it
c     was before using the subroutine shrink.

c     On entry:
c     n    - integer
c            dimension of full space;

c     nr   - integer
c            dimension of reduced space;

c     free - integer (nr)
c            free variables indices on original space;

c     x    - double precision (n)
c            point to be expanded to a space of higher dimension (full 
c            space);

c     On return:
c     x    - double precision (n)
c            point expanded to a space of higher dimension (full space).

      subroutine expand(n, nr, free, x)

      integer n, nr, free(nr)
      double precision x(n)

      integer i, aux
      double precision temp

      do i = nr, 1, -1
         aux = free(i)
         if (i .ne. aux) then
            temp = x(i)
            x(i) = x(aux)
            x(aux) = temp
         end if
      end do
      
      end 

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Reduct n x n symmetric matrix H to reduced space of dimension nr.
c     This is done by rearranging rows and columns of H in a way that 
c     all the indices corresponding to the free variables are stored in 
c     the submatrix H(1:nr;1:nr).

c     On entry:
c     n    - integer
c            dimension of full space;

c     nr   - integer
c            dimension of reduced space;

c     free - integer (nr)
c            free variables indices on original space;

c     NMAX - integer
c            leading dimension of matrix H;

c     H    - double precision (NMAX,n)
c            matrix to be reduced to the space of smaller dimension 
c            (reduced space);

c     On return:
c     H    - double precision (NMAX,n)
c            matrix reduced to space of dimension nr (reduced space).

      subroutine shrinkM(n, nr, free, NMAX, H)

      integer n, nr, free(nr), NMAX
      double precision H(NMAX,n)

      integer i, j, aux
      double precision temp

      do j = 1, n
         do i = 1, nr
            aux = free(i)
            if (i .ne. aux) then
               temp = H(i,j)
               H(i,j) = H(aux,j)
               H(aux,j) = temp
            end if
         end do
      end do

      do j = 1, nr
         aux = free(j)
         if (j .ne. aux) then
            do i = 1, nr
               temp = H(i,j)
               H(i,j) = H(i,aux)
               H(i,aux) = temp
            end do
         end if
      end do
      
      end 

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Evaluate objective function at x. Because x is in the reduced 
c     space, it is not possible to simply call subroutine evalF. It is 
c     necessary to copy to the end of x the values of the n - nr fixed 
c     variables (which are in the last n - nr positions of xc), 
c     rearrange vector x so it is in the full space, use subroutine 
c     evalF and then shrink x back.

c     On entry:
c     n    - integer
c            dimension of full space;

c     nr   - integer
c            dimension of reduced space;

c     free - integer (nr)
c            free variables indices on original space;

c     x    - double precision (n)
c            point in the reduced space where the Hessian evaluation is 
c            required;

c     xc   - double precision (n)
c            point that constains x fixes variables in its n - nr last 
c            positions;

c     On return:
c     fx   - double precision
c            function value evaluated at x.

      subroutine calcF(nr,free,x,n,xc,m,lambda,rho,equatn,fx,status)

      integer m, n, nr, free(nr), status
      double precision x(n), xc(n), lambda(m), rho(m), fx
      logical equatn(m)

      integer i

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

      call expand(n, nr, free, x)
      call evalal(n,x,m,lambda,rho,equatn,fx,status)
      call shrink(n, nr, free, x)

      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Evaluate gradient at x. Because x is in the reduced space, it is 
c     not possible to simply call subroutine evalG. It is necessary to 
c     copy to the end of x the values of the n - nr fixed variables 
c     (which are in the last n - nr positions of xc), rearrange vector x 
c     so it is in the full space, use subroutine evalG and then shrink 
c     both x and g (so both are in the reduced space).

c     On entry:
c     n    - integer
c            dimension of full space;

c     nr   - integer
c            dimension of reduced space;

c     free - integer (nr)
c            free variables indices on original space;

c     x    - double precision (n)
c            point in the reduced space where the Hessian evaluation is 
c            required;

c     xc   - double precision (n)
c            point that constains x fixes variables in its n - nr last 
c            positions;

c     On return:
c     g    - double precision (n)
c            function gradient evaluated at x in the reduced space.

      subroutine calcG(nr,free,x,n,xc,m,lambda,rho,equatn,g,status)

      integer m, n, nr, free(nr), status
      double precision x(n), xc(n), g(n), lambda(m), rho(m)
      logical equatn(m)

      integer i

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

      call expand(n, nr, free, x)
      call evalnal(n,x,m,lambda,rho,equatn,g,status)
      call shrink(n, nr, free, x)
      call shrink(n, nr, free, g)

      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Evaluate Hessian at x. Because x is in the reduced space, it is
c     not possible to simply call subroutine evalH. It is necessary to 
c     copy to the end of x the values of the n - nr fixed variables 
c     (which are in the last n - nr positions of xc), rearrange vector x 
c     so it is in the full space, use subroutine evalH and then shrink
c     x, H and diag (so they are all in the reduced space).

c     On entry:
c     n    - integer
c            dimension of full space;

c     nr   - integer
c            dimension of reduced space;

c     free - integer (nr)
c            free variables indices on original space;

c     x    - double precision (n)
c            point in the reduced space where the Hessian evaluation is 
c            required;

c     xc   - double precision (n)
c            point that constains x fixes variables in its n - nr last 
c            positions;

c     NMAX - integer
c            leading dimension of matrix H;

c     On return:
c     H    - double precision (NMAX,n)
c            function Hessian evaluated at x in the reduced space (H is 
c            stored in the lower triangular part);

c     diag - double precision (n)
c            Hessian matrix diagonal in the reduced space.

      subroutine calcH(nr,free,x,n,xc,m,lambda,rho,equatn,NMAX,H,diag,
     +status)

      integer m, n, nr, free(nr), NMAX, status
      double precision x(n), xc(n), H(NMAX,n), diag(n), lambda(m),rho(m)
      logical equatn(m)

      integer i

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

      call expand(n, nr, free, x)
      call evalhal(n,x,m,lambda,rho,equatn,nmax,H,diag,status)
      call shrink(n, nr, free, x)
      call shrink(n, nr, free, diag)
      call shrinkM(n, nr, free, NMAX, H)

      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Evaluate the quadratic model phi(x) = 1/2 x^T H x + g^T x.

c     On entry:
c     nr   - integer
c            dimension of the reduced space;

c     x    - double precision (nr)
c            point where the quadratic model is to be calculated;

c     g    - double precision (nr)
c            vector used to define the quadratic model (see above);

c     NMAX - integer
c            leading dimension of matrix H;

c     H    - double precision (NMAX,nr)
c            symmetric lower triangular matrix used to define the 
c            quadratic model (see above);

c     diag - double precision (nr)
c            diagonal of matrix H;

c     On return:
c     phi  - double precision
c            quadratic model value at x;

c     Auxiliary:
c     tmp  - double precision (nr)
c            vector in R^nr.

      subroutine quad(nr, x, g, NMAX, H, diag, tmp, phi)
      
      integer nr, NMAX
      double precision x(nr), g(nr), H(NMAX,nr), diag(nr), tmp(nr), phi

      integer i, j
      
      do i = 1, nr
         tmp(i) = 0.0d0
      end do
      
      do j = 1, nr
         do i = 1, j-1
            tmp(i) = tmp(i) + H(j,i) * x(j)
         end do
         
         tmp(i) = tmp(i) + diag(i) * x(j)
         
         do i = j+1, nr
            tmp(i) = tmp(i) + H(i,j) * x(j)
         end do
      end do
      
      phi = 0.0d0
      do i = 1, nr
         phi = phi + tmp(i) * x(i)
      end do
      
      phi = phi / 2
      
      do i = 1, nr
         phi = phi + g(i) * x(i)
      end do
      
      end
      
c     ******************************************************************
c     ******************************************************************

c     Last update of SPG: March 22nd, 2005.
c     See report of modifications at the end of this file.

c     ------------------------------------------------------------------

c     An iteration of the Spectral Projected Gradient method. Calculate 
c     a descent direction d and a point y = x + t*d such that the value 
c     of the objective function at y is sufficiently smaller than its 
c     value at initial point x. If inner = true, the working space is 
c     the full one. If inner = false, the working space is the reduced 
c     space of the free variables.

c     On entry:
c     n        - integer
c                dimension of full space;

c     nr       - integer
c                dimension of reduced space;

c     x        - double precision (n)
c                initial point;

c     xant     - double precision (n)
c                candidate solution calculated by previous iteration. It 
c                is used to calculate the spectral step sigma. This 
c                vector is altered in this subroutine;

c     free     - integer (nr)
c                indices of free variables;

c     l        - double precision (n)
c                lower bounds on x;

c     u        - double precision (n)
c                upper bounds on x;

c     g        - double precision (n)
c                gradient at x;

c     gant     - double precision (n)
c                gradient at xant, used to calculate the spectral step 
c                sigma;

c     fx       - double precision
c                function value at x;

c     alpha    - double precision
c                Armijo condition constant;

c     sigmamin - double precision
c                lower bound on spectral step;

c     sigmamax - double precision
c                upper bound on spectral step;

c     mu1      - double precision
c                safeguard for line search (lower);

c     mu2      - double precision
c                safeguard for line search (upper);

c     epsrel   - double precision
c                used to calculate the first spectral step;

c     epsabs   - double precision
c                used to calculate the first spectral step;

c     nit      - integer
c                number of current interation;

c     inner    - logical
c                true iff x is in the reduced space;

c     maxfc    - integer
c                maximum allowed number of function evaluations;

c     On return:
c     y        - double precision (n)
c                objetive function at y is sufficiently smaller than at 
c                initial point x; 

c     t        - double precision
c                scalar that multiplies d so that y = x + t*d;

c     fy       - double precision
c                function value at y; 

c     nef      - integer
c                total number of function evaluations;

c     nefSPG   - integer
c                number of function evaluations performed by SPG;

c     neg      - integer
c                total number of gradient evaluations;

c     ittype   - character *4
c                iteration type:
c                iSPG, for inner SPG;
c                gSPG, for global SPG;

c     status   - integer
c                6, if maximum allowed number of function evaluations is 
c                   achieved;
c              -90, if an error occurred during function evaluation;
c              -92, if an error occurred during gradient evaluation;
c              -94, if an error occurred during hessian evaluation;
 
c     Auxiliaries:
c     xant     - double precision (n)
c                used to calculate difference between vectors;

c     d        - double precision (n)
c                descent direction.

      subroutine SPG(n, nr, x, m, lambda, rho, equatn,
     *xant, free, l, u, g, gant, d, y, t, fx,
     *fy, alpha, sigmamin, sigmamax, mu1, mu2, epsrel, epsabs, nit, nef,
     *nefSPG, neg, ittype, inner, maxfc, status)
      
      integer m, n, nr, free(nr), nit, nef, nefSPG, neg, maxfc, status
      double precision x(n), xant(n), l(n), u(n), g(n), gant(n), d(n),
     *y(n), t, fx, fy, alpha, sigmamin, sigmamax, mu1, mu2, epsrel, 
     *epsabs, lambda(m), rho(m)
      character*4 ittype
      logical inner, equatn(m)

      integer i
      double precision sts, sty, sigma, tt, gtd, nix, tsmall, temp
      
      if (inner) then
         ittype = 'iSPG'
      else
         ittype = 'gSPG'
      end if

c     To calculate the spectral step sigma, SPG needs the products s^Ts 
c     and s^Ty, where s_k = x_k - x_{k-1} and y_k = g(x_k) - g(x_{k-1}) 
c     (g is the gradient). s_k will be placed in xant and y_k in y.
 
c     Calculate s and y for the first iteration. As in the first 
c     iteration x_{k-1} is not defined, we calculate x' as x - tsmall*g 
c     and, then, set s_0 = x' - x_0 and y_0 = g(x') - g(x_0).
      
      if (nit .eq. 0) then
         nix = 0.0d0
         do i = 1, n
            nix = dmax1(nix, dabs(x(i)))
         end do
         
         tsmall = dmax1(epsrel*nix, epsabs)

         do i = 1, n
            xant(i) = x(i) - tsmall * g(i)
         end do
         
         if (inner) then
            call calcG(nr,free,xant,n,xant,m,lambda,rho,equatn,gant,
     +      status)
         else 
            call evalnal(n,xant,m,lambda,rho,equatn,gant,status)
         end if

         if (status .lt. 0) then
            return
         end if

         neg = neg + 1
         
         do i = 1, n
            xant(i) = - tsmall * g(i)
         end do

         do i = 1, n
            y(i) = gant(i) - g(i)
         end do

      else
         do i = 1, n
            xant(i) = x(i) - xant(i)
         end do
         
         do i = 1, n
            y(i) = g(i) - gant(i)
         end do
      end if

      sts = 0.0d0
      sty = 0.0d0

      do i = 1, n
         sts = sts + xant(i)**2
         sty = sty + xant(i)*y(i)
      end do

c     Calculate the spectral step sigma.

      if (sty .le. 0.0d0) then
         sigma = sigmamax
      else
         temp = sts / sty
         sigma = dmin1(sigmamax, dmax1(sigmamin, temp))
      end if

c     Initial value of t, used in the line search.
      
      t = 1.0d0
      
c     Calculate the first solution candidate. Note that, as y = x + t*d, 
c     t = 1 and d = P(x - sigma*g) - x, we have 
c     y = x + P(x - sigma*g) - x = P(x - sigma*g).

      do i = 1, nr
         y(i) = x(i) - sigma * g(i)
      end do

      call project(nr, y, l, u)

c     Calculate the descent direction d = P(x - sigma*g) - x. As the 
c     first solution candidate was set to y = P(x - sigma*g), we have 
c     d = y - x.

      do i = 1, nr
         d(i) = y(i) - x(i)
      end do

c     Calculate the product g^Td.

      gtd = 0.0d0
      do i = 1, nr
         gtd = gtd + g(i) * d(i)
      end do
      
 2    if (inner) then
         call calcF(nr,free,y,n,x,m,lambda,rho,equatn,fy,status)
      else
         call evalal(n,y,m,lambda,rho,equatn,fy,status)
      end if

      if (status .lt. 0) then
         return
      end if
      
      nef = nef + 1
      nefSPG = nefSPG + 1
      if (nef .ge. maxfc) then
         status = 6
         return
      end if
      
c     If Armijo is satisfied, the point is accepted. Otherwise, a line 
c     search using d as the search direction is made until Armijo is 
c     satisfied.
 
      if (fy .gt. fx + alpha * t * gtd) then
         if (t .lt. mu1) then
            t = t / 2.0d0
         else
            tt = - ((t**2) * gtd) / (2.0d0 * (fy - fx - t * gtd))
            if ((tt .lt. mu1) .or. (tt .gt. mu2*t)) then
               tt = t / 2.0d0
            end if
            t = tt
         end if
         
         do i = 1, nr
            y(i) = x(i) + t * d(i)
         end do

         go to 2
      end if

      end

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

c     Last update of MEQB (or any of its dependencies): 
c     March 22nd, 2005.
c     See report of modifications at the end of this file.

c     ------------------------------------------------------------------

c     Method to minimize quadratic functions subjected to ||w|| <= delta

c     minimize    psi(w) = 1/2 w^TBw + g^Tw 
c     subject to  ||w|| <= delta

c     Method described in "Computing a trust region step", by More and 
c     Sorensen.

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     The main ideia of this method is to find a positive scalar \lamms 
c     that is a zero of the function
 
c     phi(\lamms) = 1/||p|| - 1/delta,

c     where p is the solution of the linear system

c     (B + \lamms I)p = -g.

c     Note that symmetric matrix B, vector g and positive real number 
c     delta are presented in the minimization problem above. I is the 
c     identity matrix.

c     The method used to find the zero of that function is basically the 
c     Newton method to find roots.

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     On entry:
c     n        - integer
c                dimension;

c     g        - double precision (n)
c                vector used to define the quadratic function;

c     NMAX     - integer
c                leading dimension of matrix B;

c     B        - double precision (NMAX,n)
c                symmetric n x n matrix used to define the quadratic 
c                function. Just the lower triangular portion must be 
c                defined. Both the upper triangular portion and the 
c                diagonal will be modified to store the Cholesky 
c                decomposition of (B + lI) = R^TR;

c     diag     - double precision (n)
c                diagonal of matrix B;

c     l        - double precision
c                initial value for \lamms;

c     delta    - double precision
c                radius;

c     sigma1   - double precision
c                allowed error for convergence criteria 1 and 2;

c     sigma2   - double precision
c                allowed error for convergence criteria 3 (hard case);

c     eps      - double precision
c                allowed error;

c     iprint   - integer
c                data printing options:
c                -1, nothing is printed;
c                 0, just the data from final iteration is printed;
c                >0, data obtained at each iprint iterations is printed;

c     On return:
c     p        - double precision (n)
c                solution to problem
c                minimize     psi(w)
c                subjected to ||w|| <= delta;

c     l        - double precision
c                \lamms value that gives p as a solution to the 
c                minimization problem, because p is also solution to 
c                (B + \lamms I)p = -g;

c     nch      - integer
c                number of Cholesky decompositions;

c     convCrit - integer
c                stores which of the 3 convergence criteria was 
c                satisfied; 0 if both g and B are null;

c     Auxiliaries:
c     t        - double precision (n)
c                used to calculate new value of l;

c     z        - double precision (n)
c                \tau * z is the approximation of the eigenvector
c                associated with the smallest eigenvalue of B (called 
c                \lamms _1).
c                z Euclidian-norm is 1.

c     tmp      - double precision (n)
c                used to store matrix-vector products.

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------
            
      subroutine MEQB(n, g, NMAX, B, diag, p, t, z, tmp, l, delta, 
     *sigma1, sigma2, eps, iprint, nch, convCrit)
                                 
      integer n, NMAX, iprint, nch, convCrit
      double precision g(n), B(NMAX,NMAX), diag(n), p(n), t(n), z(n), 
     *tmp(n), l, delta, sigma1, sigma2, eps

c     Local variables:
c     i, j     - integer
c                counters;

c     sgn      - integer
c                store a number sign;

c     idx      - integer
c                index of B where Cholesky decomposition of (B + lI) 
c                fails;

c     nit      - integer
c                number of iterations;

c     ll       - double precision
c                lower bound on l in each iteration;

c     lu       - double precision
c                upper bound on l in each iteration;

c     ls       - double precision
c                lower bound on -\lamms _1 (minus the smallest 
c                eigenvalue of B);

c     nt2      - double precision
c                squared Euclidian-norm of vector t;

c     np       - double precision
c                Euclidian-norm of vector p;

c     np2      - double precision
c                squared Euclidian-norm of vector p;

c     ng       - double precision
c                Euclidian-norm of vector g;

c     n1B      - double precision
c                1-norm of matrix B;

c     temp     - double precision
c                temporary variable used in some calculations;

c     d        - double precision
c                when B + lI is not positive definite, the submatrix S 
c                of (B + lI) with indices from 1 to idx and with d added 
c                to the element in (idx,idx), is a singular matrix.
c                This value is used to update ls;

c     nu2      - double precision
c                squared Euclidian-norm of u. Vector u is such that 
c                submatrix S (see variable d) multiplied by u is equal 
c                to the null vector;

c     nRp2     - double precision
c                squared Euclidian-norm of the product of matrix R and 
c                vector p;

c     nRz2     - double precision
c                squared Euclidian-norm of the product of matrix R and 
c                vector z;

c     tau      - double precision
c                scalar such that tau * z is the approximation of the 
c                eigenvector associated with the smallest eigenvalue of 
c                matrix B.

c     pTz      - double precision
c                result of the dot product of vectors p and z;

c     delta2   - double precision
c                squared value of delta;

c     dp       - logical
c                true iff matrix (B + lI) is positive definite;

      integer i, j, sgn, idx, nit
      double precision ll, lu, ls, nt2, np, np2, ng, n1B, temp, d, nu2, 
     *nRp2, nRz2, tau, pTz, delta2
      logical dp

      delta2 = delta**2
      convCrit = -1

c     step 1: initialize ls (lower bound on l) with max{-bii}, where bii 
c             are the elements of the diagonal of B.

      ls = -diag(1)
      do i = 2, n
         if (-diag(i) .gt. ls) then
            ls = -diag(i)
         end if
      end do

c     step 2: initialize ll (lower bound on l) with 
c             max{0, ls, ||g||/delta - ||B||1}, where ||B||1 is the 
c             1-norm of the matrix B.

      ng = 0.0d0
      do i = 1, n
         ng = ng + g(i)**2
      end do
      ng = dsqrt(ng)
      
      n1B = 0.0d0
      do j = 1, n
         temp = 0.0d0
         do i = 1, j-1
            temp = temp + dabs(B(j,i))
         end do

         temp = temp + dabs(diag(i))

         do i = j+1, n
            temp = temp + dabs(B(i,j))
         end do

         if (temp .gt. n1B) then
            n1B = temp
         end if
      end do

      ll = (ng / delta) - n1B 
      ll = dmax1(0.0d0, ll)
      ll = dmax1(ls, ll)

c     step 3: initialize lu (upper bound on l) with ||g||/delta + ||B||1.

      lu = (ng / delta) + n1B

c     If the matrix is null, there is nothing to be done.

      if ((ll .eq. 0.0d0) .and. (lu .eq. 0.0d0) .and. 
     *    (ls .eq. 0.0d0)) then
         convCrit = 0
         go to 21
      end if

c     step 4: initialize iteration counter.

      nit = 0

c     step 5: safeguard of l (ensures that l is bigger than ll).

 5    l = dmax1(l, ll)

c     step 6: safeguard of l (ensures that l is smaller than lu).

      l = dmin1(l, lu)

c     step 7: safeguard of l.

      if (l .le. ls) then
         l = dmax1(1.0d-3 * lu, dsqrt(ll*lu))
      end if

c     step 8: try to use the Cholesky decomposition: (B +lI) = R^TR.
c             If the decomposition is successfull, R is stored in the 
c             upper triangular portion of B (including the diagonal) and 
c             dp is set to true.
c             If the decomposition fails, d and idx are set as explained 
c             before, dp is set to false, and the Euclidian-norm of u is 
c             calculated (see explanation of variable nu2).

      call Cholesky(n, NMAX, B, l, diag, d, idx, dp)
      nch = nch + 1

c     In this case (B + lI) is not positive definite, and d and idx are 
c     calculated. Because p cannot be calculated (it is not possible to 
c     solve the system using Cholesky factorization), the values of l, 
c     ll and ls are updated, the iteration counter is increased and a
c     new iteration is started.

      if (.not. dp) then

c        Print information (current iteration).

         if ( iprint .ge. 3 ) then
            write(*, 989) 
            write(*, 980) nit
            write(*, 981) ls,ll,lu,l
            write(*, 988)
 
            write(10,989) 
            write(10,980) nit
            write(10,981) ls,ll,lu,l
            write(10,988) 
         end if
         
         call calcU(n, idx, NMAX, B, p, nu2)
         ll = dmax1(l, ll)
         ls = dmax1(l + (d / nu2), ls)
         ll = dmax1(ll, ls)
         l = ls
         nit = nit + 1
         go to 5
      end if

c     step 9: solve R^TRp = -g for p and calculate the squared 
c             Euclidian-norm of p. 

      call solveSys(n, NMAX, B, g, p, nRp2)

      np2 = 0.0d0
      do i = 1, n
         np2 = np2 + p(i)**2
      end do
      np = dsqrt(np2)

      
c     step 10: calculate z and tau, where tau * z is the approximation 
c              of the eigenvector associated with the smallest 
c              eigenvalue of B.

      if (np .lt. delta) then

c        Calculate z.

         call calcZ(n, NMAX, B, z)

c        Calculate z Euclidian-norm.

         temp = 0.0d0
         do i = 1, n
            temp = temp + z(i)**2
         end do
         temp = dsqrt(temp)

c        Divide z by its norm.

         do i = 1, n
            z(i) = z(i) / temp
         end do

c        Calculate the product of matrix R by the vector z.

         do i = 1, n
            tmp(i) = 0.0d0
         end do
         
         do j = 1, n
            do i = 1, j
               tmp(i) = tmp(i) + B(i,j) * z(j)
            end do
         end do
         
c        Calculate the squared Euclidian-norm of the product Rz.

         nRz2 = 0.0d0
         do i = 1, n
            nRz2 = nRz2 + tmp(i)**2
         end do

c        Calculate tau.

         pTz = 0.0d0
         do i = 1, n
            pTz = pTz + p(i) * z(i)
         end do

         if (pTz .lt. 0.0d0) then
            sgn = -1
         else
            sgn = 1
         end if
         
         temp = delta2 - np2

         tau = (pTz**2) + temp
         tau = temp / (pTz + sgn * dsqrt(tau))         
 
      end if

c     Print informations (current iteration).

      if ( iprint .ge. 3 ) then
         write(*, 989) 
         write(*, 980) nit
         write(*, 981) ls,ll,lu,l
         write(*, 982) np

         write(10,989) 
         write(10,980) nit
         write(10,981) ls,ll,lu,l
         write(10,982) np

         if ( np .lt. delta ) then
            write(*, 983) sqrt( np2 + 2 * tau * pTz + tau ** 2 )
            write(10,983) sqrt( np2 + 2 * tau * pTz + tau ** 2 )
         else
            write(*, 983) np
            write(10,983) np
         end if

         write(*, 985) delta
         write(10,985) delta
      end if

c     steps 11 and 12: update ll, lu and ls

      if (np .lt. delta) then
         lu = dmin1(l, lu)
         ls = dmax1(l - nRz2, ls)
      else
         ll = dmax1(l, ll)
      end if

c     step 13: update ls when B + lI is not positive definite.
c              This was done right after the Cholesky decomposition 
c              failure.

c     step 14: update ll.

      ll = dmax1(ll, ls)

c     step 15: convergence test.

      if ((dabs(l) .le. eps) .and. (np .le. delta)) then
         convCrit = 1
         go to 21
      end if

c     step 16: second convergence test.

      if (dabs(delta - np) .le. sigma1*delta) then
         convCrit = 2
      end if

c     step 17: convergence test for the hard case.

      temp = nRp2 + l * delta2
      temp = dmax1(sigma2, temp)
      temp = temp * sigma1 * (2 - sigma1)
      
      if ((np .lt. delta) .and. ((tau*tau*nRz2 .le. temp) 
     *.or. ((lu - ll) .le. eps))) then     
         convCrit = convCrit + 3
         go to 20
      end if

      if (convCrit .eq. 2) then
         go to 21
      end if

c     step 21: Calculate l to be used in the next iteration.

      if ((dabs(ng) .gt. eps) .or. 
     *   ((l .eq. ll) .and. (ls .le. ll))) then 

c        Solve R^Tt = p for t

         do i = 1, n
            t(i) = p(i)
         end do

         do j = 1, n           
            do i = 1, j-1
               t(j) = t(j) - B(i,j) * t(i)
            end do
            t(j) = t(j) / B(j,j)
         end do

c        Calculate t squared Euclidian-norm.
         
         nt2 = 0.0d0
         do i = 1, n
            nt2 = nt2 + t(i)**2
         end do

c        Update l using Newton's method update.
        
         l = l + (np2 / nt2) * ((np - delta) / delta)
      else
         l = ls
      end if

c     step 22: update iteration counter.

      nit = nit + 1

c     step 23: start a new iteration.

      go to 5

c     steps 18, 19 and 20:

c     The solution is given by p in 3 cases:
c     - if the first convergence criterion is satisfied;
c     - if only the second convergence criterion is satisfied;
c     - if both the second and the third convergence criteria are 
c       satisfied, but the squared Euclidian-norm of R(tau*z) is 
c       strictly bigger than l*(delta2 - np2)

c     The solution is given by p + tau*z when:
c     - just the third convergence criterion is satisfied;
c     - both the second and the third convergence criteria are 
c       satisfied, but the squared Euclidian-norm of R(tau*z) is smaller 
c       or equal to l*(delta2 - np2)

 20   temp = (tau * tau * nRz2) - l * (delta2 - np2)

      if ((convCrit .eq. 3) .or. 
     *((convCrit .eq. 5) .and. (temp .le. 0.0d0))) then
         
         np2 = 0.0d0
         do i = 1, n
            p(i) = p(i) + tau * z(i)
            np2 = np2 + p(i)**2
         end do
         np = dsqrt(np2)

         convCrit = 3
      else
         convCrit = 2
      end if

c     Print informations.

 21   if ( iprint .ge. 3 ) then
         write(*, 990) 
         write(*, 986) nit
         write(*, 981) ls,ll,lu,l
         write(*, 987) convCrit
         write(*, 984) np
         write(*, 985) delta
         write(*, 990) 

         write(10,990) 
         write(10,986) nit
         write(10,981) ls,ll,lu,l
         write(10,987) convCrit
         write(10,984) np
         write(10,985) delta
         write(10,990) 
      end if

c     Printing formats.

 980  format(/,1X,'Iteration: ',I7)
 981  format(/,1X,'ls = ',1PD11.4,
     +       /,1X,'ll = ',1PD11.4,
     +       /,1X,'lu = ',1PD11.4,
     +       /,1X,'l  = ',1PD11.4)
 982  format(/,1X,'Euclidian-norm of p: ',1PD11.4)
 983  format(/,1X,'Euclidian-norm of p + z: ',1PD11.4)
 984  format(/,1X,'Euclidian-norm of s (solution): ',1PD11.4)
 985  format(/,1X,'delta: ',1PD11.4)
 986  format(/,1X,'Number of iterations: ',I7)
 987  format(/,1X,'Convergence criterion: ',I7)
 988  format(/,1X,'B + lI is not positive definite!')
 989  format(/,1X,'*****')
 990  format(/,1X,'-------------------------------------------')
      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Cholesky decomposition.
c     If B + lI is a symmetric positive definite n x n matrix, B will 
c     store R, where R is such that B + lI = R^TR and dp is set to true.
c     Because R is upper triangular, just this portion of B will be 
c     changed. (remember the original diagonal of B is stored at diag).
c     If B + lI is not positive definite, dp is set to false and both d 
c     and idx can be set (look the definition of the latter two below).
c   
c     On entry:
c     n    - integer
c            dimension;

c     NMAX - integer
c            leading dimension of matrix B;

c     B    - double precision (NMAX,n)
c            symmetric n x n matrix (just the lower triangular portion 
c            must be defined);

c     l    - double precision
c            used to calculate B + lI;

c     diag - double precision (n)
c            diagonal of matrix B;

c     On return:
c     B    - double precision (NMAX,n)
c            Cholesky decomposition of B + lI, when it is possible to 
c            calculate, stored in the upper triangular part of B 
c            (including the diagonal);

c     idx  - integer
c            when B + lI is not positive definite, idx is the index 
c            where the decomposition fails;

c     d    - double precision
c            when B + lI is not positive definite, the submatrix S of 
c            (B + lI) with indices from 1 to idx and with d added to the 
c            element in (idx,idx), is a singular matrix;

c     dp   - true iff B + lI is positive definite.

      subroutine Cholesky(n, NMAX, B, l, diag, d, idx, dp)
                         
      integer n, NMAX, idx
      double precision B(NMAX,NMAX), l, diag(n), d
      logical dp

      integer i, j, k
      double precision sum

      do j = 1, n
         do i = j, n
            if (i .eq. j) then
               sum = diag(j) + l
            else
               sum = B(i,j)
            end if

            do k = j - 1, 1, -1
               sum = sum - B(k,j) * B(k,i)
            end do

            if (i .eq. j) then
               if (sum .le. 0.0d0) then
                  d = - sum
                  idx = j
                  dp = .false.
                  return
               else 
                  B(j,j) = dsqrt(sum)
               end if
            else 
               B(j,i) = sum / B(j,j)
            end if
         end do
      end do
      dp = .true.
      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Solve system Ax = -b, where A is the Cholesky decomposition of a 
c     matrix B + lI. Just the upper triangular portion of A will be used
c     (as well as its diagonal).

c     On entry:
c     n    - integer
c            dimension;

c     NMAX - integer
c            leading dimension of matrix A;

c     A    - double precision (NMAX,n)
c            matrix n x n, contaning R in its upper triangular portion;

c     b    - double precision (n)
c            vector;

c     On return:
c     x   - double precision (n)
c           system solution;

c     ny2 - double precision
c           y squared Euclidian-norm, where y = Rx.

      subroutine solveSys(n, NMAX, A, b, x, ny2)
                            
      integer n, NMAX
      double precision A(NMAX,NMAX), b(n), x(n), ny2

      integer i, j

      do i = 1, n
         x(i) = -b(i)
      end do
      
      do j = 1, n           
         do i = 1, j-1
            x(j) = x(j) - A(i,j) * x(i)
         end do
         x(j) = x(j) / A(j,j)
      end do
      
      ny2 = 0.0d0
      do i = 1, n
         ny2 = ny2 + x(i)**2
      end do

      do j = n, 1, -1
         x(j) = x(j) / A(j,j)
         do i = 1, j-1
            x(i) = x(i) - A(i,j) * x(j)
         end do
      end do

      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------

c     Solve a linear system of the form (A + ekek^Td)u = 0, where A is, 
c     in fact, the Cholesky decomposition of B + lI until entry k-1 in
c     its upper portion (row number >= column number). 
c     u is a vector of n positions with u(k) = 1 and u(j) = 0, for j > k.

c     On entry:
c     n    - integer
c            dimension of A;

c     k    - integer
c            dimension of e;

c     NMAX - integer
c            leading dimension if matrix A;

c     A    - double precision (NMAX,n)
c            symmetric matrix n x n;

c     On return:
c     u    - double precision (n)
c            system solution;

c     nu2  - double precision
c            u squared Euclidian-norm.

      subroutine calcU(n, k, NMAX, A, u, nu2)
                           
      integer n, k, NMAX
      double precision A(NMAX,n), u(n), nu2

      integer i, j

      do i = 1, k-1
         u(i) = -A(i,k)
      end do
      
      do j = 1, k-1           
         do i = 1, j-1
            u(j) = u(j) - A(i,j) * u(i)
         end do
         u(j) = u(j) / A(j,j)
      end do
      
      do j = k-1, 1, -1
         u(j) = u(j) / A(j,j)
         do i = 1, j-1
            u(i) = u(i) - A(i,j) * u(j)
         end do
      end do
 
      u(k) = 1

      nu2 = 0.0d0
      do i = 1, k
         nu2 = nu2 + u(i)**2
      end do

      end

c     ------------------------------------------------------------------
c     ------------------------------------------------------------------
      
c     Implementation of the tecnique presented by Cline, Moler, Stewart 
c     and Wilkinson to estimate the condition number of a matrix.
c     This tecnique is used by the More-Sorensen method to calculate an
c     approximation to the eigenvector associated to the smallest 
c     eigenvalue of a matrix B (\lamms_1). 
c     In this tecnique, when \lamms approaches -\lamms_1, \|Rz\| 
c     approaches 0. This garantees that z is an approximation to the 
c     wanted eigenvector.
c     Basically, it solves R^Ty = e, choosing e(k) as 1 or -1 (whatever 
c     gives maximum local growth of y). Then, it solves Rz = y.

c     On entry:
c     n    - integer
c            dimension of R;
c     NMAX - integer
c            leading dimension of R;
c     R    - double precision (NMAX,n)
c            upper triangular matrix (Cholesky decomposition of B + lI);

c     On return:
c     z    - double precision (n)
c            approximation of the eigenvector of B associated to \lamms_1.

      subroutine calcZ(n, NMAX, R, z)
      
      integer NMAX, n
      double precision R(NMAX,NMAX), z(n)
      
      double precision w, wk, wkm, ek
      double precision s, sm
      integer i, k
      
      ek = 1.0d0
      
      do i = 1, n
         z(i) = 0.0d0
      end do
      
      do k = 1, n
         if (z(k) .ne. 0.0d0) then
            ek = dsign(ek,-z(k))
         end if
         
         if (dabs(ek-z(k)) .gt. dabs(R(k,k))) then
            s = dabs(R(k,k)) / dabs(ek-z(k))            
            do i = 1, n
               z(i) = s * z(i)
            end do
            ek = s * ek
         end if
         
         wk = ek - z(k)
         wkm = -ek - z(k)
         s = dabs(wk)
         sm = dabs(wkm)
         
         if (R(k,k) .eq. 0.0d0) then 
            wk = 1
            wkm = 1
         else
            wk = wk / R(k,k)
            wkm = wkm / R(k,k)
         end if
         
         if (k .eq. n) then
            go to 10
         end if
         
         do i = k+1, n
            sm = sm + dabs(z(i) + wkm * R(k,i))
            z(i) = z(i) + wk * R(k,i)
            s = s + dabs(z(i))
         end do
         
         if (s .lt. sm) then
            w = wkm - wk
            wk = wkm
            
            do i = k+1, n
               z(i) = z(i) + w * R(k,i)
            end do
         end if
         
 10      z(k) = wk
      end do
      
c     Divide z by its 1-norm to avoid overflow.
      
      s = 0.0d0
      do i = 1, n
         s = s + dabs(z(i))
      end do
      
      do i = 1, n
         z(i) = z(i) / s
      end do
      
c     Solve Rz = y.
      
      do k = n, 1, -1
         z(k) = z(k) / R(k,k)                  
         do i = 1, k-1
            z(i) = z(i) - R(i,k) * z(k)
         end do
      end do
      
      end
      
c     ******************************************************************
c     ******************************************************************

c     Modifications report:

c     March 22nd, 2005.
c     1. Translation of all comments, subroutine names and printed 
c        messages to English;
c     2. Modification of MEQB printing format.

c     April 4th, 2005.
c     1. Change of initial value of iteration counter in BETRA 
c        (from 1 to 0);
c     2. Verification of error after function, gradient and hessian 
c        evaluation added;
c     3. New stopping criterion added: maximum number of function 
c        evaluations;
c     4. Parameters eps and epsdelta were dropped;
c     5. Relative comparison was introduced to decide if a point is 
c        interior to the box of constraints. Parameter epsbd was dropped;
c     6. New parameters were added: mseps and phieps. The first one is a 
c        tolerance for More-Sorensen algorithm (MEQB). The second is 
c        used to decided if the quadratic model at a certain point is 
c        too small.
