c*********************************************************************
c* SUBROUTINE BOX9903 *
c*********************************************************************
c
c
c Developed by the Optimization Group at the State University of
c Campinas
c
c Date: March 1999 (99/03)
C
C This subroutine solves the problem
C
C Minimize f(x)
c s. t. h(x) = 0
c A x = b (1)
C l <= x <= u
c
c where
c
c f: \R^n ---> \R
c
c h: \R^n ----> \R^{nonlc}
c
c A \in \R^{m x n}
c
c and f and h are differentiable functions. That is to say, the
c smooth nonlinear programming problem with m linear equality
c constraints, nonlc nonlinear equality constraints,
c and bounds on the variables.
c
c This code was written with the aim of dealing with large-scale
c problems, for this reason it does not impose storing or factorization
c of matrices of any kind.
c The outer algorithm for solving (1) is the Augmented Lagrangian
c method. Namely, at each outer-iteration the following
c box-constrained problem is solved
c
c
c Minimize L(x)
c (2)
c subject to l <= x <= u
c
c
c where L(x) =
c
c f(x) + mu^T h(x) +
c
c + (1/2) Sum_{i=1}^{nonlc} ra(i) [h(x)]_i^2 +
c
c + lambda^T (A x - b) + (ro/2) ||A x - b||_2^2
c
c L(x) is called ``the Augmented Lagrangian'' associated to (1).
c
c
c
c Therefore, mu is the vector of Lagrange multipliers associated
c to nonlinear equality constraints, lambda is the vector of Lagrange
c multipliers associated to linear constraints, ra is a vector of
c penalty parameters associated to nonlinear equality constraints,
c and ro is a single penalty parameter associated to linear constraints.
c
c A description of the method used in box9903 can be found in
c
c N. Krejic, J. M. Martinez, M. P. Mello and E. A. Pilotta
c [1999]: "Validation of an Augmented Lagrangian algorithm with a
c Gauss-Newton Hessian approximation using a set of hard-spheres
c problems", to appear in
c Computational Optimization and Applications.
c
c
c
c For solving (2) at each outer iteration it is used the subroutine
c BOX, enclosed within this piece of software. BOX solves (2) using
c the method described in:
c
C A. Friedlander, J. M. Martinez and S. A. Santos [1994]: "A new
C trust-region strategy for box-constrained minimization",
C Applied Mathematics and Optimization 30 (1994) 235-266.
c
c
c At each outer iteration BOX will be called trying to find
c a minimizer of (2) up to the precision eps (or epsd).
c
c
c The main work performed by box9903 consists in calling BOX
c and modifying the penalty parameters ro (associated to the
c linear constraints) and ra (associated to the nonlinear
c equality constraints). Moreover, the Lagrange multipliers mu
c (associated to the nonlinear equality constraints), and lambda
c (associated to the linear constraints) are also updated in box9903
c at each outer iteration.
c
c*******************************************************************
c* *
c* *
c* PARAMETERS OF BOX9903 *
c* *
c* *
c*******************************************************************
c
c n: number of variables.
c
c m: number of linear constraints.
c
c nonlc: number of nonlinear equality constraints.
c
c roinic: (input) if linear constraints are present (m > 0) you must
c set roinic equal to the initial value that you wish for
c the penalty parameter associated to the linear constraints.
c
c rainic: (input) is a vector of nonlc positions. If nonlc = 0
c you must dimension rainic with only one position. If nonlinear
c equality constraints are present (nonlc > 0) you
c must set rainic equal to the initial values that you wish for
c the penalty parameters associated to nonlinear constraints.
c
c ra: (auxiliar) is a vector of nonlc positions. If nonlc = 0
c you must dimension ra with only one position. If nonlinear
c equality constraints are present (nonlc > 0) ra is used by
c this subroutine to store the penalty parameters associated to
c each nonlinear equality constraint.
c
c facro: (input) if linear constraints are present (m > 0) you must
c set facro >= 1 equal to the factor by which you want to multiply
c penalty parameter ro after each outer iteration.
c
c facra: (input) if nonlinear equality constraints are present
c (nonlc > 0) you must set facra >= 1 equal to the factor by which
c you want to multiply the penalty parameters ra(i) after
c each outer iteration.
c
c indivi: (input) if you set indivi = 1, penalty parameters
c associated to each nonlinear equality constraint
c will be increased separately.
c If indivi = 0, penalty parameters associated to
c nonlinear equality constraints will be increased
c together, so that the process will work as if it existed only
c one penalty parameter associated to nonlinear equality constraints.
c
c epslin: (input) if linear constraints are present you must
c set epslin equal to the level of feasibility that you want associated
c to these constraints. In other words, a point x will be considered
c ``feasible'' with respect to linear constraints when
c
c ||A x - b||_\infty <= epslin
c
c epsnon: (input) if nonlinear equality constraints
c are present you must set epsnon equal to the level of feasibility
c that you want associated to these constraints.
c In other words, a point x will be considered
c ``feasible'' with respect to nonlinear equality constraints when
c
c || h(x) ||_\infty <= epsnon
c
c ratfea: (input) If the feasibility at the new outer iteration
c is better than ratfea*(the feasibility at the previous iteration)
c the penalty parameter will not be increased. ratfea must be greater
c than or equal to zero.
c
c epsinic: (input) This is the convergence criterion that will be
c used by Box at the first outer iteration.
c
c epsdiv: (input) If epsinic > epsfin, the convergence criterion
c used by Box at successive outer iterations is the previous one
c over epsdiv. For example, if Box used eps = 10^{-2} at the outer
c iteration kon, it will use eps/epsdiv at the following one.
c However if the new eps is less than epsfin, it is replaced by
c epsfin. epsdiv must be greater than 1.
c
c epsfin: (input) This is the convergence criterion that will
c be used by Box at the final outer iteration.
c If you set epsinic = epsfin
c the same convergence criterion will be used at all the iterations.
c
c
c maxout: (input) maximum number of outer iterations (calls to BOX).
c
c iprint: (input) parameter by means of which you can say what
c you want to be printed automatically by box9903, both
c in the screen and in a file (unity 10). If iprint < 0,
c nothing will be printed. If iprint = 0, only information
c at the initial point and at the final point will be
c printed. In all the other cases, information will also
c be printed every iprint outer iterations.
c
c konout: (output) effective number of outer iterations performed.
c konint: (output) total number of inner iterations (iterations of BOX)
c nefint: (output) total number of evaluations of f(x)
c itqint: (output) total number of quacan iterations
c mvpint: (output) total number of matrix-vector products performed
c in quacan
c
c hresan: (input) auxiliar double precision vector of at least
c nonlc positions (1 position if nonlc = 0)
c
c ierla: output parameter that indicates what happened in the
c execution of box9903. If ierla = 0, a feasible point, according
c to the prescribed tolerances epslin and epsnon, has been found.
c If there are no constraints in the problem ierla will be necessarily
c equal to 0. The final diagnostic concerning the execution of
c box9903 must be complemented looking at the BOX output parameter
c istop. In fact, if istop = 0 and ierla = 0 you can guarantee
c that a point that satisfies first-order optimality
c conditions of (1) has been found. Other situations are doubtful.
c If ierla = 1, the allowed number of outer iterations (calls to
c BOX) has been exhausted without achieving feasibility.
c
c ilag: (input) can take the values 0 and 1. If ilag = 0 the estimates
c of the Lagrange multipliers at each iteration are zero. In
c other words, in this case one is using a pure penalty method.
c If ilag = 1 first order Lagrange multipliers are estimated
c at each outer iteration.
c
c iguess: (input) can take the values 0 and 1. If iguess = 1
c this means that you have a procedure to guess a good approxima-
c tion to the solution of the quadratic problem that must be
c solved at each iteration of BOX. In this case, you must write
c the way to calculate your guess in the subroutine guess. More
c about this parameter and subroutine guess in the comments of
c the subroutine BOX.
c
c x: double precision vector with n positions. On input, initial
c approximation to the solution. On output, best solution obtained.
c
c l, u : on input, double precision vectors that define the
c bounds of Problem (1).
c
c hres: double precision auxiliar vector with at least nonlc positions.
c On output it contains the values of h(x) mentioned in (1).
c
c res: double precision auxiliar vector with m positions.
c On exit, when m is not 0, residual vector A x - b.
c
c grad: double precision auxiliar vector with n positions.
c On exit, g is the gradient of L (the augmented Lagrangian)
c at the final point.
c
c fx: on output, functional value f(x).
c
C
C iprbox : control for printing in subroutine BOX. If iprbox< 0 nothing
c is printed in BOX. If iprbox = 0, only a
C final summary is printed. If iprbox > 0, information every iprbox
C iterations is printed.
C
c
c
C mu: auxiliar double precision input vector with nonlc components,
c (1 component if nonlc = 0)
c that contributes to define the objective function. It stores the
c Lagrange multipliers associated to nonlinear equality constraints.
c
c lambda: auxiliar double precision input vector with m components
c (1 component if m = 0)
c that contributes to define the Augmented Lagrangian. It stores the
c Lagrange multipliers associated to linear equality constraints.
c
c
c---------------------------------------------------------------------------
c
c SUBROUTINES phi AND MAPRO, THAT DEFINE Phi(x),
c
c
c ITS GRADIENT, AND ITS HESSIAN APPROXIMATION
c
c
C---------------------------------------------------------------------------
c
C The user must provide a subroutine named phi for computing
C the function value
c
c Phi(x) =
c
c
c = f(x) + mu^T h(x) +
c
c + (1/2) Sum_{i=1}^{nonlc} ra(i) h_i(x)^2 +
c
c the gradient of Phi(x) and the information relative to the
c Hessian approximation of Phi(x).
c
C The calling sequence of phi is
C
C call phi (n, nonlc, x, mu, ra, flagra, fx, grad, hres, hess, inhess, modo )
C
C where n is the number of variables, nonlc is the number of
c components of h,
c x is the point where the function is calculated within phi,
c mu and ra are the other input parameters necessary
c for computing Phi(x).
c As output, phi must produce flagra (Phi(x)), fx (f(x)),
C grad, the gradient of the function Phi computed at x, hres
c (the nonlinear vector of residuals h(x)).
c Finally, hess and inhess inhess are arrays that can contain
c the information relative to the Hessian approximation.
C You can choose the name for the subroutine phi. For this
C reason, phi is a parameter of box9903. In the calling
C program, the name actually used must be declared EXTERNAL.
C If modo = 1, phi must compute ONLY the function
c values flagra and fx and the
c vector of nonlinear residuals hres,
c which, of course, are by-products of the computation of the
c function value flagra.
C If modo = 2, phi must compute ONLY the gradient at
c the point x and the information relative to the Hessian.
c Now, you can take
c advantage of the common calculations to function and gradient-
c Hessian learning that a call to phi with modo=2 is always
c preceeded, in BOX, by a call to phi with modo=1 at the same
c point. So, you can prepare information inside phi in such a way
c that calculations for modo=2 are already done, and you do not
c need to begin from scratch. In particular, for computing the
c gradient, you will need the nonlinear residuals hres.
c So, you can use it freely.
C You can choose the data structure for the Hessian
c information provided that it fits in the arrays inhess and
c hess. This information must be compatible with the coding of
c user supplied subroutine mapro, commented below.
c In previous versions of BOX it was assumed that the communication
c of Hessian information between phi and mapro was coded through
c COMMON statements. We think that passage through parameters is
c less prone to errors, so we decided to introduce the new
c parameters hess and inhess in 1995. However, if you
c find easier to make the communication by COMMON, you can do it
c disregarding the existence of hess and inhess. However, do not
c forget to dimension hess and inhess.
c
c Helping your memory:
c
c The gradient of the objective function Phi(x) considered here is
c
c \nabla f(x) + Sum_{i=1}^{nonlc} [mu(i) + ra(i) h_i(x)] \nabla h_i(x)
c
c
c The Hessian is
c
c \nabla^2 f(x) + Sum_{i=1}^{nonlc} ra(i) \nabla h_i(x) \nabla h_i(x)^T +
c
c + Sum_{i=1}^{nonlc} [mu(i) + ra(i) h_i(x)] \nabla^2 h_i(x)
C
c
C You must also provide the subroutine mapro. As in the case
C of phi , you can choose the name of this subroutine, which, for
C this reason, must be declared EXTERNAL in the calling program.
C mapro is called only from quacan. As a result, it is also
C a parameter of quacan, and it is declared EXTERNAL in BOX
C subroutine. The calling sequence of mapro is
C
C call mapro (n, nind, ind, u, v, hess, inhess )
C
C where n (the number of variables) and u (a double precision vector
C of n positions) are the main inputs, and v, a double precision
c vector of n positions is the output. v must be the product H * u,
c where H is the current Hessian approximation of the function Phi(x).
c That is, H is the last computed Hessian approximation of Phi(x)
c within the subroutine phi.
c Therefore, mapro must be coded in such a way that the structure
c given for this matrix within phi is compatible with formulae for
c the product.
c Moreover, if nind < n (both are input parameters), the integer vector
c ind, with nind positions, contains the indices where the input
c u is different from zero. So, you can take advantage of this
c information to write the matrix-vector product, but if you do
c not know how to do it, simply ignore nind and ind and write the
c matrix vector product as if u were a dense vector. The algorithm
c will be the same, but taking advantage of nind and ind makes it
c faster.
c Many times, you will find the task of coding the information
c relative to the Hessian very cumbersome. You can use a
c ``directional finite difference'' version of mapro instead of
c products with true Hessians. The idea is that, at the current
c point x, the product H u is the limit of
c
c [Gradient(x + t u) - Gradient(x) ] / t
c
c Using this directional derivative idea, you can code mapro
c passing, within hess, the current point x and the current
c gradient g to mapro. Then, you use, say,
c
c t = max ( 1.d-20, 1.d-8 ||x||_\infty ) / || d ||_\infty
c
c provided that d \neq 0, of course.
c
c So, in this case you evaluate the Gradient
c at the auxiliary point x + t u and
c finally, you use the quotient above to approximate H u. There
c are some risks using this device, but you can try.
c
C--------------------------------------------------------------------------
c PARAMETERS hess AND inhess
C--------------------------------------------------------------------------
C
C From the point of view of box9903, hess and inhess are two auxiliar
c arrays. hess is double precision and inhess is integer. They can
c be used in the subroutines phi and mapro to transmit Hessian information
c between them (from phi to mapro) as explained above. The
c number of positions reserved to hess and inhess must be sufficient
c for the information that you want to transmit.
c
c----------------------------------------------------------------------------
c
c SUBROUTINE GUESS
c
c----------------------------------------------------------------------------
C
C At each iteration of BOX, bound-constrained quadratic problems of
c the form
c
c Minimize grad(x^k) d + (1/2) d^T B_k d
c subject to lint <= d <= uint
c
c are solved by the subroutine Quacan. Here, B_k is an approximation
c of the Hessian of the objective function used in BOX (the whole
c augmented Lagrangian). Sometimes, experienced users know how to
c compute good approximations to the solution of this subproblem,
c independently of Quacan. In this case, the user sets the input
c parameter iguess equal to 1. In fact, in that case, Quacan is going
c to take the approximation ``computed by the user'' as initial appro-
c ximation for its own job. When you set iguess = 1, you must code
c your way to compute the initial approximation to the solution of
c the bound-constrained quadratic subproblem in the subroutine guess.
c (In fact, since the name of the subroutine is a parameter, you can
c give other name to the subroutine, as is the case with Phi and Mapro.
c Consequently, you must also declare external the name of this subroutine
c in your main code.)
c The calling sequence of guess must be:
c
c call guess(n, x, grad, ro, ra, lambda, mu, lint, uint, d)
c
c where you must consider that n is the number of variables, x is
c current point in the procedure commanded by BOX, grad is the gra-
c dient of the objective function (Augmented Lagrangian) at x, lint
c is the vector of lower bounds for d, uint is the vector of upper
c bounds of d, ro is the penalty parameter associated to linear
c constraints, ra is the vector of penalty parameters associated
c to nonlinear equality constraints, lambda is the vector
c of multipliers associated to linear constraints, mu is the
c vector of multipliers associated to nonlinear equality constraints,
c Finally, d (the output) is the approximate solution to
c the bound-constrained quadratic subproblem computed by guess.
c All real parameters in guess must be double precision. Probably,
c for computing your approximation you could need additional information
c on the current point. In this case, use common statements to make
c the communication between Phi and Guess. Probably, all the relevant
c information necessary for Guess has already been computed in Phi,
c with modo=2.
c Even if you are setting iguess=0, you must include a subroutine
c called guess in your main stack but, in this case, it will not be
c called by BOX, so, it can consist only of the statements return
c and end in order of not to disturb the compilation. It is enough,
c in the case ``iguess=0'', that you include the statements
c subroutine guess
c return
c end
c
c
c-------------------------------------------------------------------------
c
c PARAMETERS THAT DEFINE THE LINEAR AUGMENTED LAGRANGIAN TERM
C
c
C lambda^T (A x - b) + (ro/2) || A x - b ||^2
c
c
c
c m, lambda, ro, ispar, nelem, inda, a, b are the input
c parameters that define the augmented Lagrangian term.
c
c m : number of rows of the matrix A. Set m=0 if there is no
c Lagrangian term at all
c
c
c lambda : vector of m double precision positions given in the
c definition of the objective function L.
c
c ro : real*8 number given in the definition of L
c
c ispar : input integer parameter that can take the value 0 or 1
c If ispar = 0, this means that the m x n matrix A is stored
c in the array a columnwise, that is a(1,1), a(2,1), a(3,1),
c ... a(m,1), a(1,2), a(2,2),...,a(m,2)...
c You can use ispar=0 when there is not advantage on taking
c into account sparsity of A. In other words, when you judge
c that A is dense.
c If ispar=1, this means that you give only the nonzero elements
c of A, in the way explained below.
c
c nelem : this input parameter is used only when ispar = 1. It
c is the number of ``nonzero'' elements that you are reporting
c for the matrix A.
c
c inda : this input parameter is used only when ispar = 1. In
c this case, it is an integer nelem x 2 array where you
c must store, in any order, the pairs (i , j) such that the
c corresponding element a_{ij} of the matrix A is different
c from zero.
c
c a : If ispar=0, this array contains the entries of the matrix
c A, stored columnwise (see ispar above). If ispar =1, this
c array contains the nonzero elements of A in such a way
c that a(k) = a_{ij} when the row k of inda is (i, j).
c
c b : real*8 vector of m positions, mentioned in the definition
c of L(x) above.
c
c ----------------------------------------------------------------------------
c A SET OF PARAMETERS USED ONLY IN THE SUBROUTINE BOX
c ----------------------------------------------------------------------------
C
C accuracy : accuracy demanded for quacan at its first call
c within a box-iteration. The user must give this parameter belonging
C to (0,1). We recommend to set accuracy = 0.1. If the
c objective function f is quadratic we recommend accuracy=1.d-8
c so that most the work will be executed by quacan. A small
c accuracy is also desirable if the problem is badly scaled.
c Moreover, if the number of variables is small we prefer accuracy
c very small since in this case the quadratic solver is not very
c costly.
c
C acumenor : accuracy demanded for quacan at subsequent calls
c (not the first)
c within a box-iteration. The user must give this parameter belonging
C to (0,1). We recommend to set acumenor=0.5.
c If the number of variables is small we prefer acumenor
c very small since in this case the quadratic solver is not very
c costly.
c
c relarm, absarm: two input parameters that define when the trial
c point is accepted. Namely, if Pred < 0 is the value of the approxi-
c mating quadratic at its approximate minimizer and x+ is the corres-
c ponding trial point, x+ is accepted as new point when
c
c L(x+) \leq L(x) + max ( relarm * Pred, -absarm)
c
c We recommend relarm = 0.1, absarm = 1.d-6
c
c
c
C epsd : input parameter representing a tolerance for the trust
C region radius used in BOX. See usage in the description of
c (istop = 4) below.
C When BOX is unable to produce a decrease even with a trust
C region of size epsd, this possibly means that we are very close
C to a solution, which was not detected by the convergence criterion
C based on the projected gradient because, perhaps, that stopping
c parameter was given excessively small.
c Roughly speaking, BOX returns because of the
C (istop = 4) criterion, when the quadratic model restricted
C to a ball of radius epsd is unable to produce a decrease of the
C function. This is typical when a local minimizer occurs, but we
C must be careful. We recommend to set epsd = 1.d-8. However, warning!
c this is a dimensional parameter, you must take into account
c the scaling (unit measures) of your problem. Roughly speaking,
c epsd shoud be a NEGLIGIBLE DISTANCE in the x-space.
C
C nafmax : Maximum allowed number of function evaluations at each
c call of BOX. See description of (istop = 2) below.
C
C itmax : Maximum allowed number of iterations at each call of BOX.
C See description of (istop = 3) below.
C
C par : This parameter is used to define the size of the new trust
C region radius and must be given by the user belonging to (1,10).
C We recommend to set par = 4. See usage in the description of
C deltamin below.
c
C delta0 : This input parameter is the initial trust region radius at
C the beginning of iteration zero. Warning! Delta0 is dimensional.
c It should represent the radius of a ball centered on the initial
c point where you reasonably expect to be the solution.
C
C deltamin : This input parameter allows the subroutine to define
C the trust region radius. In fact, at the beginning of an
C iteration, the trust region radius is defined by BOX to be
c not less than deltamin. Due to this choice, the
C trust region radius can inherit more or less information from the
C previous iteration, according to the size of deltamin. deltamin can
C be a fraction of the diameter of the problem (1), provided that
C there are no artificial bounds. We recommend deltamin=delta0/100.
C
C istop : This output parameter tells what happened in the last call
C subroutine BOX (that is to say, the last outer iteration),
c according to the following conventions:
C istop = 0 -> Convergence :
C ||projected gradient||_\infty \leq eps.
C istop = 1 -> Convergence :
C L(x) \leq ftol
C istop = 2 -> It is achieved the maximum allowed number of
C function evaluations (nafmax).
C istop = 3 -> It is achieved the maximum allowed number of
C iterations (Itmax).
C istop = 4 -> Strangled point, i.e. trust region radius less
C than epsd.
c
C istop = 5 -> Data error. A lower bound is greater than an
C upper bound.
c
c istop = 6 -> Some input parameter is out of its prescribed
c bounds
c
C istop = 7 -> The progress (L(x_k) - L(x_{k-1})) in the last
C iteration is less than bont * (the maximum
C progress obtained in previous iterations)
c during kbont consecutive iterations.
c
C istop = 8 -> The order of the norm of the continuous
c projected gradient did not change
c during kauts consecutive
c iterations. Probably, we are asking for an
c exagerately small norm of continuous projected
c gradient for declaring convergence.
c
C bont, kbont : This parameters play the same role as dont and kdont do
c for quacan. We also suggest to set bont = 0.01 together with kbont = 5.
C
c kauts : If the order of the norm of the current
c continuous projected gradient did not change during
c kauts consecutive iterations the execution
c stops with istop = 8. Recommended: kaustky = 10. In any
c case kauts must be greater than or equal to 1.
c
c You must also set the following input quacan parameters:
c
c dont : Positive number, less than 1, used for the second convergence
c criterion by quacan, according to the description above. See
c also the comments of the subroutine quacan.
c
c eta: On input, number belonging to (0, 1) .
c If eta is close to 1, faces are fully exploited before abandoning
c them. If eta is close to 0, abandoning faces is frequent. We
c recommend eta = 0.9.
c
c factor: On input, real number greater than 1 (we recommend factor=5)
c intended to improve the current trial guess through extrapo-
c lation procedures. It is used in the subroutine Double. See
c its comments.
c mitqu: maximum number of iterations allowed on each call to quacan.
c maxmv: maximum number of matrix-vector products allowed on each call
c to quacan.
c ipqua: printing parameter of quacan. If ipqua < 0, nothing is printed.
c If ipqua = 0, only a final summary is printed
c of each call of quacan.
c Otherwise, information on quacan
c iterations are printed every
c ipqua iterations.
c
c
c lint, uint, s, gq, gc, xn, gn, d: auxiliar double
c precision vector with n positions.
c bd: auxiliar vector that must be dimensioned with n+1 positions
c when m = 0 an with at least 2n+m real*8 positions when m > 0
c
c ingc, infree: integer auxiliar vectors of n positions.
c
c ----------------------------------------------------------------------------
c END OF ``A SET OF PARAMETERS USED ONLY IN THE SUBROUTINE BOX''
c ----------------------------------------------------------------------------
c
c For consults relative to the use of this subroutine please contact:
c
c J. M. Martinez
c IMECC - UNICAMP
c C.P. 6065
c 13081-970 Campinas SP
c Brazil
c E-mail: martinez@ime.unicamp.br
c
subroutine box9903 (n, m, nonlc, roinic, ra, rainic, indivi,
* facro, facra, ilag, iguess, epslin, epsnon, ratfea,
* maxout, konout, konint, nefint, itqint, mvpint, ierla, iprint,
* x, l , u, fx, hres, hresan, res, grad,
* accuracy, acumenor, relarm, absarm, epsinic, epsdiv, epsfin,
* epsd, nafmax, itmax, par, delta0, deltamin, istop,
* bont, kbont, kauts, dont, kdont, eta, factor,
* mitqu, maxmv, phi, mapro, guess, hess, inhess, iprbox,
* ipqua, lint, uint, s, gq, gc, bd, xn, gn, d, ingc, infree,
* mu, lambda, ispar, nelem, inda, a, b)
implicit double precision(a-h,o-z)
external phi, mapro, guess
double precision l(n), lint(n), lambda(*), mu(*)
dimension x(n), u(n), grad(n), hess(*), uint(n)
dimension s(n), gq(n), gc(n), bd(*), xn(n), gn(n)
dimension d(n), ingc(n), infree(n), inda(*)
dimension a(*), b(*)
dimension res(*), hres(*)
dimension inhess(*)
dimension rainic(*), ra(*)
dimension hresan(*)
do i = 1, n
if(x(i).lt.l(i)) x(i) = l(i)
if(x(i).gt.u(i)) x(i) = u(i)
end do
c
if(iprint.ge.0)then
n10 = min0(n, 10)
m10 = min0(m, 10)
nlc10 = min0(nonlc, 10)
endif
c
c We set the initial vectors of multipliers equal to zero
c
c We set the initial penalty parameters equal to the ones
c indicated by the user
c
if(nonlc.gt.0)then
do i = 1, nonlc
ra(i) = rainic(i)
mu(i) = 0.d0
end do
endif
if(m.gt.0) then
do i =1,m
lambda(i) = 0.d0
end do
ro = roinic
endif
c
c konout is the counter of outer iterations
c
konout = 0
c
c
c Initialization of counters of inner iterations, function evaluations,
c quacan-iterations and matrix-vector products
c
konint = 0
nefint = 0
itqint = 0
mvpint = 0
ftol = - 1.d30
c
c Evaluations at the initial point
c
call phi (n, nonlc, x, mu, ra, flagra,
* fx, grad, hres, hess, inhess, 1 )
if(m.gt.0) call pena(n, x, flagra, grad, res,
* m, ilag, lambda, ro, ispar, nelem, inda, a, b, bd, 1)
if(nonlc.gt.0) then
hnor = 0.d0
do i = 1, nonlc
hnor = dmax1(hnor, dabs(hres(i)))
hresan(i) = dabs(hres(i))
end do
hnoran = hnor
endif
if(m.gt.0) then
anor = 0.d0
do i = 1, m
anor = dmax1(anor, dabs(res(i)))
end do
anoran = anor
endif
if(iprint.ge.0) then
write(10,*)
write(*,*)
write(10,*)' OUTPUT OF BOX9903'
write(10,*)
c
write(*,*)' OUTPUT OF BOX9903'
write(*,*)
write(*,*)
call imprimir (n10, m10, nlc10, nonlc,
* konout, x, istop, konint, nefint, eps,
* itqint, mvpint, anor, hnor,
* ra, ro, mu, lambda, flagra, fx)
endif
c
c
c At the first call to box, we set eps = epsinic
c
eps = epsinic
10 call box (n, x, l , u, iguess, flagra, fx, hres, res,
* grad, accuracy, acumenor, relarm, absarm, eps, epsd, ftol,
* nafmax, itmax, kont, naf, iqua, mvp, par,
* delta0, deltamin, istop, bont, kbont, kauts,
* dont, kdont, eta, factor,
* mitqu, maxmv, phi, mapro, guess, hess, inhess, iprbox, ipqua,
* lint, uint, s, gq, gc, bd, xn, gn, d, ingc, infree,
* nonlc, mu, ra,
* m, lambda, ro, ispar, nelem, inda, a, b)
c
c Increment counters of outer iterations,
c inner (BOX)-iterations, evaluations, quacan-iterations
c and matrix-vector products
c
konout = konout + 1
konint = konint + kont
nefint = nefint + naf
itqint = itqint + iqua
mvpint = mvpint + mvp
c
hnor=0.d0
if(nonlc.gt.0)then
do i=1,nonlc
hnor=dmax1(hnor, dabs(hres(i)))
end do
endif
anor=0.d0
if(m.gt.0)then
do i =1,m
anor = dmax1(anor, dabs(res(i)))
end do
endif
c If ilag = 1, update the Lagrange multiplier estimators
c
if(ilag.eq.1)then
if(nonlc.gt.0)then
do i=1,nonlc
mu(i)=mu(i) + ra(i) * hres(i)
end do
endif
c
c
if(m.gt.0)then
do i=1,m
lambda(i) = lambda(i) + ro * res(i)
end do
endif
endif
c Printing session
if(iprint.gt.0. and. mod(konout, iprint).eq.0)then
call imprimir (n10, m10, nlc10, nonlc,
* konout, x, istop, konint, nefint, eps,
* itqint, mvpint, anor, hnor,
* ra, ro, mu, lambda, flagra, fx)
endif
c Terminate if feasibility was achieved
if(eps.le.epsfin) then
if((hnor.le.epsnon.or.nonlc.le.0).and.
* (anor.le.epslin.or.m.eq.0))then
ierla = 0
if(iprint.lt.0)return
write(10,*)
write(10,*)' Feasibility was achieved.'
write(*,*)
write(*,*)' Feasibility was achieved.'
if(iprint.gt.0. and. mod(konout, iprint).eq.0)return
call imprimir (n10, m10, nlc10, nonlc,
* konout, x, istop, konint, nefint, eps,
* itqint, mvpint, anor, hnor,
* ra, ro, mu, lambda, flagra, fx)
return
endif
endif
c Terminate if number of allowed outer iterations is exhausted
if(konout.ge.maxout) then
ierla = 1
if(iprint.lt.0)return
write(10,*)
write(10,*)' Maximum number of outer iterations achieved.'
write(*,*)
write(*,*)' Maximum number of outer iterations achieved.'
if(iprint.gt.0. and. mod(konout, iprint).eq.0)return
call imprimir (n10, m10, nlc10, nonlc,
* konout, x, istop, konint, nefint, eps,
* itqint, mvpint, anor, hnor,
* ra, ro, mu, lambda, flagra, fx)
return
endif
c
c Updating penalty parameters corresponding to
c the nonlinear equality constraints
c
if(nonlc.gt.0)then
if(indivi.eq.1)then
do i = 1, nonlc
c Update, separately, penalty parameters corresponding to nonlinear
c constraints
c Increase the penalty parameter
c only if feasibility was not greatly improved
if(dabs(hres(i)) .gt. ratfea * hresan(i)) then
ra(i) = facra * ra(i)
endif
c Save the vector of constraints at the present iteration, to be used
c in the following one
hresan(i) = dabs(hres(i))
end do
else
c Case indivi .ne. 1. All the nonlinear penalty parameters are
c updated together
if(hnor.gt. ratfea * hnoran) then
c At subsequent iterations update the penalty parameter if feasibility
c was not greatly improved
do i=1,nonlc
ra(i) = facra * ra(i)
end do
endif
c Save the norm of the vector of constraints to be used at subsequent
c iterations
hnoran = hnor
endif
endif
if(m.gt.0)then
c Update the penalty parameter ro if we are at the first iteration or
c if linear feasibility was not greatly improved
if(anor .gt. ratfea * anoran) then
ro = facro * ro
endif
c Save the norm of the vector of linear constraints, to be used at
c subsequent iterations
anoran = anor
endif
c
c Compute the convergence parameter that will be used by Box at
c the next outer iteration, as an interpolation between epsinic and
c epsfin
c
eps = dmax1(eps/epsdiv, epsfin)
c Perform a new outer iteration
goto 10
end