subroutine imprimir(n10, m10, nlc10, nonlc, * konout, x, istop, konint, nefint, eps, * itqint, mvpint, anor, hnor, * ra, ro, mu, lambda, flagra, fx) c Printing subroutine of box9903 implicit double precision (a-h, o-z) dimension x(*) double precision lambda(*), mu(*), ra(*) write(10,*) write(10,*) ' Outer Augmented Lagrangian iteration :', konout write(10,*) ' (First 10) components of the current point:' write(10,*) (x(i), i=1, n10) write(*,*) write(*,*) ' Outer Augmented Lagrangian iteration :', konout write(*,*) ' (First 10) components of the current point:' write(*,*) (x(i), i=1, n10) if(m10.gt.0)then if(konout.ne.0) then write(10,*) ' (First 10) components of lambda:' write(10,*) (lambda(i), i=1, m10) write(10,*) ' Penalty parameter used for linear constraints:', * ro write(*,*) ' (First 10) components of lambda:' write(*,*) (lambda(i), i=1, m10) write(*,*) ' Penalty parameter used for linear constraints:', * ro endif endif if(nlc10.gt.0)then if(konout.ne.0) then write(10,*) ' (First 10) nonlinear equality multipliers:' write(10,*) (mu(i), i=1, nlc10) write(10,*) ' Penalty parameters used for nonlinear equality', * ' constraints (first 10):' write(10,*) (ra(i), i=1, nlc10) write(*,*) ' (First 10) nonlinear equality multipliers:' write(*,*) (mu(i), i=1, nlc10) write(*,*) ' Penalty parameters used for nonlinear equality', * ' constraints (first 10):' write(*,*) (ra(i), i=1, nlc10) ramax = 0.d0 ramin = ra(1) do i = 1, nonlc ramin = dmin1(ramin, ra(i)) ramax = dmax1(ramax, ra(i)) end do write(10,*)' Maximum nonlinear equality penalty parameter', * ' used:',ramax write(10,*)' Minimum nonlinear equality penalty parameter', * ' used:',ramin write(*,*)' Maximum nonlinear equality penalty parameter', * ' used:',ramax write(*,*)' Minimum nonlinear equality penalty parameter', * ' used:',ramin endif endif if(konout.ne.0) then write(10,*)' Diagnostic of BOX at this iteration:', istop write(10,*)' Convergence epsilon used in Box:', eps write(10,*)' Total number of BOX-iterations:', konint write(10,*)' Total number of function evaluations:', nefint write(10,*)' Total number of quacan-iterations:', itqint write(10,*)' Total number of matrix-vector products:', mvpint write(*,*)' Diagnostic of BOX at this iteration:', istop write(*,*)' Convergence epsilon used in Box:', eps write(*,*)' Total number of BOX-iterations:', konint write(*,*)' Total number of function evaluations:', nefint write(*,*)' Total number of quacan-iterations:', itqint write(*,*)' Total number of matrix-vector products:', mvpint endif if(m10.gt.0) then write(*,*)' Sup norm of A x - b:', anor write(10,*)' Sup norm of A x - b:', anor endif if(nlc10.gt.0) then write(*,*)' Maximum violation of nonlinear equality ', * ' constraints :', hnor write(10,*)' Maximum violation of nonlinear equality ', * ' constraints :', hnor endif write(10,*)' Value of the objective function:', fx write(10,*)' Value of the augmented Lagrangian function:',flagra write(*,*)' Value of the objective function:', fx write(*,*)' Value of the augmented Lagrangian function:',flagra return end C-------------------------------------------------------------------- C February 1997 C subroutine box (n, x, l , u, iguess, flagra, fx, hres, * res, grad, accuracy, acumenor, relarm, absarm, * eps, epsd, ftol, nafmax, itmax, kont, naf, iquatot, mvptot, * par, delta0, deltamin, istop, bont, kbont, kauts, * dont, kdont, eta, factor, mitqu, maxmv, phi, mapro, guess, hess, * inhess, ipr, 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 This subroutine solves the problem C C Minimize L(x) C s.t. l <= x <= u (1) 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 c + lambda^T (A x - b) + (ro/2) ||A x - b||_2^2 c c c f:\R^n \to \R, h:\R^n \to \R^{nonlc}, c are differentiable, A is an m x n matrix, b \in \R^m, c mu \in \R^{nonlc}, ra \in \R^nonlc}, c lambda \in \R^m and ro is a real number. c c c (That is, L(x) has the form of an augmented Lagrangian of f, for c nonlinear constraints of the form h(x) = 0, and c linear constraints of the form A x = b, where A is a matrix of c m rows and n columns.) C C We assume that l and u are finite n-dimensional bounds. C C The method used to solve (1) is 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 This algorithm has been modified several times after 1992. c C The main features of the method are the following: C C a) At each iteration k, the following quadratic approximation C of L(x) is considered: C C q(x) = L(x^k) + grad(x^k)^T (x - x^k) + (1/2)(x - x^k)^T B_k (x - x^k) C C where grad(x) is the gradient of L at x, and B_k is an L-Hessian C approximation (not necessarily positive semidefinite). A sequence C of subproblems of the form C C Minimize q(x) s.t. box_k (2) C C is solved approximately until the solution of (2) satisfies a C decrease condition. Then this solution is chosen as x^{k+1}. C The size of box_k depends on the natural bounds of the C problem , as well as on a trust-region size. The trust region C size is decreased each time the sufficient decrease condition C is not satisfied. C C b) The approximate solution of (2) uses the subroutine quacan, C which is a bound constrained not necessarily positive c semidefinite quadratic solver. C Given an initial approximation, quacan guarantees C to produce an approximation to the solution of (2) where the C two-norm of the projected gradient of q , is as small as C desired. In fact, the degree of accuracy desirable in the solution C of (2) is one of the main parameters of BOX, and the C best choice is currently under discussion. The stopping criterion C in (2) is C C ||Projected gradient of the quadratic||_2 c <= accuracy * ||Projected gradient of L at x^k||_2 (3) C C where accuracy < 1 represents the degree of accuracy required. c At the first call of quacan at each iteration the precision used c is accuracy. At the subsequent calls of quacan within a box-iteration c we used the precision acumenor instead of accuracy. The idea c is that accuracy should be less than or equal to accumenor, since c a great precision in quacan is only justified at the first call, c when we have the chance to take (nearly) full-Newton steps. c c In april 1992 it was suggested by numerical c experiments that acc = 0.1 was a good choice. In november 1995 c it was observed that, for badly scaled problems, values of acc c close to 0 are better. The reason is that Newton's method is invariant c to scaling and the criterion (3) makes the iteration close to c the Newton iteration when acc is close to 0. c The parameter acumenor, as different from accuracy, was introduced c in September 1996. We suggest accuracy = 1.e-3, acumenor=0.1. c However, this depends on the cost of quacan calls. If the number c of variables is small, calls to quacan are cheap. In this c case we recommend accuracy = acumenor = 1.e-8. c In some experiments in september 1996 the best accuracy was 0.1 c associated to acumenor = 0.5. C In november 1995, we incorporated a second stopping criterion c for quacan. Namely, quacan stops also when the progress obtained c within the current iteration is less than dont*(the best progress c obtained by quacan iterations before the current one) during c kdont consecutive quacan iterations. Of course c this stopping criterion may be inhibited setting dont = 0. We c recommend, preliminary, to set dont = 0.01 and kdont = 5. Perhaps c kdont should also be dependent on n, to give a chance to the con- c vergence of the conjugate gradient method, that is within quacan. C In the description of the parameters of the subroutine we shall C have the opportunity of describing further characteristics of C the method. C C C REMARK! C For running this code in a PC, you will perhaps need C to set the card $large at the beginning of the main program. C The aim of such a command is to expand C memory. C C C Parameters: c C iguess: Sometimes, due to the characteristics of your problem, c you can compute a guess of the solution of the quadratic c subproblem that must be solved by Quacan, at each iteration c of BOX. If this is the case, set iguess = 1, otherwise, set c iguess = 0. If you set iguess = 1 you need to write the sub- c routine guess, that computes your especial approximation to c the solution of the quadratic subproblem. See the comments on c the subroutine guess. C C ipr : Control for printing in BOX. If ipr< 0 nothing is printed c in BOX. If ipr=0, only a C final summary is printed. If ipr > 0, information every ipr C iterations is printed. C C n : Dimension of problem C C x : On entry, initial approximation to the solution of (1), provided C by the user. Remember that l <= x <= u. If these bounds do c not hold for your x, BOX projects your initial x on the box. c On exit, x is the final approximation C to the solution of (1), obtained by BOX. C C C l , u : On entry, these two double precision vectors are C the bounds of problem (1). That is, they are the l and the u of (1) C respectively. C C flagra : On exit, value of the objective function L(x). c c fx: On exit, value of f(x). c C grad : On exit, g is the gradient of L at the final point. c c nonlc: number of components of h. (nonlc means ``nonlinear constraints'') c C mu: input vector with nonlc components, that contributes to c define the objective function. c c c ro: real*8 parameter c that contribute to define the objective function. c c ra: on input, double precision vector with at least nonlc positions c that contributes to the definition of the objective c function. If nonlc=0, you must dimension ra with 1 position. c c c mu, lambda: real*8 nonlc- and m- vectors (respectively) c that contribute to define the objective function. c c c c res: On exit, when m is not 0, residual vector A x - b. c c hres: On exit, when nonlc is not 0, nonlinear residual h(x). c c--------------------------------------------------------------------------- c c SUBROUTINES phi AND MAPRO, THAT DEFINE Phi(x), c c c ITS GRADIENT, AND ITS HESSIAN APPROXIMATION 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 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, c 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, 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 c inhess are arrays that can contain the information relative c to the Hessian approximation. C You can choose the name for the subroutine phi. For this C reason, phi is a parameter of BOX. 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 vectors of nonlinear residuals hres and gres, 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 observing 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 and gres. 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 them. 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 C 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 C must be coded in such a way that the structure given for this C matrix within phi is compatible with formulae for the product. c Moreover, if nind < n (both 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 c SUBROUTINE GUESS 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, c lambda is the vector of multipliers associated to linear constraints, c mu is the vector of multipliers associated to nonlinear equality c 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 defini- c tion of the objective function L. c 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 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.001. 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.1. 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 eps : Input parameter given by the user for deciding convergence. c Convergence is declared if the euclidean norm of the c ``continuous projected gradient'' is less than or equal to eps. c The continuous projected gradient at x is defined as the c difference between the projection on the box of x - grad(x) and c x ( Proj [ x - grad(x)] - x ) c C epsd : Input parameter given by the user which gives tolerance for trust C region radius. See usage in the description of (istop = 4) below. C When the method 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, epsg was given C excessively small. 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 ftol : Input parameter given by the user which gives tolerance for C objective function L. c See usage in the description of (istop = 1) below. C C C nafmax : Maximum allowed number of function evaluations. See C description of (istop = 2) below. C C itmax : Maximum allowed number of iterations. See description of C (istop = 3) below. C c C kont : Number of iterations performed by BOX. C C naf : Number of function evaluations performed in BOX. c c iquatot: Total number of iterations performed by the subroutine c quacan c c mvptot: Total number of matrix-vector products done within the c subroutine quacan 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 this C subroutine, 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 whenever BOX is used in a C iteratively way (e.g. in an Augmented Lagrangian or SQP c context), together with kbont = 5 . C Otherwise, we recommend to inhibit this parameter, using C bont = 0.0. 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.1. 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 implicit double precision (a-h, o-z) double precision l(n), lint(n), lambda(*), mu(*) C dimension x(n), u(n), grad(n), hess(*), inhess(*) dimension uint(n), s(n) dimension gq(n), gc(n), bd(*), xn(n), gn(n) dimension d(n), ingc(n), infree(n) dimension a(*), b(*), inda(nelem,2) dimension res(*), hres(*), ra(*) C external phi, mapro c open (unit = 10,file = 'box.sai') c We set the quacan parameter ichan equal to 1. In this way, at each c quacan call, quacan will test if the initial approximation given c is better than the null initial approximation. If it is not, the c initial approximation (increment) given to quacan will be automatically c replaced by the null vector. ichan = 1 c nrepl is a counter of the number of times that the quadratic model c was replaced by the easy-quadratic model due to insufficient decrease of c the quadratic model. c The easy-quadratic model of the function is defined here as c qeasy(d) = g^T d + ||d||^2 / 2, therefore, its solution is the c projection of -g on the trust-region box c reg = 0 nkauts = 0 nrepl = 0 if(ipr.ge.0)then write(10,*) write(10,*)' OUTPUT OF BOX' write(10,*) n10=amin0(n,10) endif c c Check whether input parameters are correct c if(n.le.0.or.accuracy.lt.0.d0.or.acumenor.lt.0.d0 * .or.eps.lt.0.d0.or.epsd.lt.0.d0. * or.nafmax.lt.0.or.itmax.lt.0.or.par.lt.1.d0.or.delta0.le.0.d0. * or.deltamin.lt.0.d0.or.bont.lt.0.d0.or.kbont.le.0.or.dont.lt. * 0.d0.or.kdont.le.0.or.eta.le.0.d0.or.eta.gt.1.d0.or.factor.lt. * 1.d0.or.mitqu.lt.0.or.maxmv.lt.0.or. * (m.gt.0.and.ispar.lt.0).or.(m.gt.0.and.ispar.gt.1) * .or.(m.gt.0.and.(ispar.eq.1.and.nelem.le.0)).or.m.lt.0.or. * nonlc.lt.0.or.kauts.lt.1) * then ier=6 if(ipr.ge.0)then write(*,*)' Some input parameter of BOX is incorrectly posed' write(10,*)' Some input parameter of BOX is incorrectly posed' endif return endif c do 1 i=1,n if(l(i).lt.u(i)) goto1 istop=6 if(ipr.ge.0)then write(10,*)' Lower bound',i,' not smaller than' write(10,*)' upper bound', i endif return 1 continue c c Project the initial point on the feasible box, if necessary do 2 i=1,n if(x(i).lt.l(i)) x(i)=l(i) if(x(i).gt.u(i)) x(i)=u(i) 2 continue c c If lambda = 0, ila is set to 0. If not, ila is set to 1. ila=0 if(m.gt.0)then do 17 i=1,m if(lambda(i).eq.0.d0)goto17 ila=1 goto18 17 continue endif C 18 kont=0 mvptot=0 irep = 0 iquatot=0 nopro=0 progress = 0.d0 C Calculates value of objective function, the gradient, and c the Hessian information at 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, ila, lambda, ro, ispar, nelem, inda, a, b, bd, 1) naf = 1 C c Here begins the loop c C C Test convergence C 3 call phi (n, nonlc, x, mu, ra, flagra, * fx, grad, hres, hess, inhess, 2 ) if(m.gt.0) call pena(n, x, flagra, grad, res, * m, ila, lambda, ro, ispar, nelem, inda, a, b, bd, 2) If (flagra .le. ftol) then istop = 1 if(ipr.ge.0)then write(*,*) write(10,*) write(*,*) ' Convergence of BOX by small function value' write(10,*)' Convergence of BOX by small function value' write(10,*)' Warning: Do Not consider the values printed' write(10,*)' below for the norms of projected gradients' write(*,*)' Warning: Do Not consider the values printed' write(*,*)' below for the norms of projected gradients' call printibox(kont,n10,x,grad,flagra,gpnor,gpcont, * naf,iquatot,mvptot,delta, pred, ared, nrepl, irep) endif return endif c c Compute the 2-norm and the sup-norm of the projected gradient gpnor2=0.d0 gpcont=0.d0 do 4 i=1,n if((x(i).gt.l(i).and.x(i).lt.u(i)) .or. * (x(i).eq.l(i).and.grad(i).lt.0.d0) .or. * (x(i).eq.u(i).and.grad(i).gt.0.d0)) then z=x(i)-grad(i) if(z.gt.u(i)) z=u(i) if(z.lt.l(i)) z=l(i) gpnor2 = gpnor2 + grad(i)**2 gpcont = gpcont + (z-x(i))**2 endif 4 continue gpnor=dsqrt(gpnor2) gpcont = dsqrt(gpcont) c c Test whether the continuous projected gradient is small enough c to declare convergence c if(gpcont.le.eps) then istop=0 if(ipr.ge.0)then write(*,*) write(10,*) write(*,*) ' Convergence of BOX by small projected gradient' write(10,*)' Convergence of BOX by small projected gradient' call printibox(kont, n10, x, grad, flagra, gpnor, gpcont, * naf, iquatot, mvptot, delta, pred, ared, nrepl, irep) endif return endif c c Test whether the number of iterations is exhausted c if(kont.ge.itmax) then istop=3 if(ipr.ge.0)then write(*,*) write(10,*) write(*,*) ' Maximum number of BOX iterations exhausted' write(10,*)' Maximum number of BOX iterations exhausted' call printibox(kont,n10,x,grad,flagra,gpnor,gpcont, * naf,iquatot,mvptot,delta, pred, ared, nrepl, irep) endif return endif c c Test whether the number of functional evaluations is exhausted c if(naf.ge.nafmax) then istop=2 if(ipr.ge.0)then write(*,*) write(10,*) write(*,*) naf,' evaluations performed, exceeding maximum' write(10,*)naf,' evaluations performed, exceeding maximum' call printibox(kont,n10,x,grad,flagra,gpnor,gpcont, * naf,iquatot,mvptot,delta, pred, ared, nrepl, irep) endif return endif c c Test whether we performed many iterations without good progress c if(kont.ne.0) then prog=fant-flagra progress=dmax1(prog,progress) if(prog.le.bont*progress)then nopro=nopro+1 if(nopro.ge.kbont) then istop=7 if(ipr.ge.0)then write(*,*) write(10,*) write(*,*) ' BOX finished by lack of progress' write(10,*) ' BOX finished by lack of progress' call printibox(kont,n10,x,grad,flagra,gpnor,gpcont, * naf,iquatot,mvptot,delta, pred, ared, nrepl, irep) endif return endif else nopro=0 endif endif fant=flagra c c Test whether we have performed kauts iterations without good reduction c of the 2-norm of the projected gradient c if(kont.ne.0)then logp = dlog10(gpcont) if(logp.ge.logpant) then nkauts = nkauts + 1 else nkauts = 0 endif if(nkauts.ge.kauts) then istop = 8 if(ipr.ge.0)then write(*,*) write(10,*) write(*,*) ' BOX finished by no enough improvement' write(*,*) ' of the continuous projected gradient' write(10,*) ' BOX finished by no enough improvement' write(10,*) ' of the continuous projected gradient' call printibox(kont,n10,x,grad,flagra,gpnor,gpcont, * naf,iquatot,mvptot,delta, pred, ared, nrepl, irep) endif return endif logpant = logp else logpant = dlog10(gpcont) endif c c A new iteration begins here if(ipr.gt.0.and.(kont.eq.0.or.mod(kont,ipr).eq.0))then call printibox(kont,n10,x,grad,flagra,gpnor,gpcont, * naf,iquatot,mvptot,delta, pred, ared, nrepl, irep) endif c Set the convergence criterion for quacan at the first call within c the current iteration epsq=accuracy*gpnor c Update delta if(kont.eq.0)then delta=delta0 else c c If the actual reduction of the merit function was c greater than half the predicted reduction we multiply c the trust region radius by the expansion parameter par c if(ared.ge.0.5d0*pred) delta = par * delta delta=dmin1(1.d30, dmax1(delta, deltamin)) endif c Compute the bounds for quacan, lint and uint line=0 do 14 i=1,n lint(i)= dmax1(-delta, l(i)-x(i)) 14 uint(i)= dmin1(delta, u(i)-x(i)) 12 if(line.eq.1)then c c There was at least one Armijo failure at this iteration c (line=1). We update the bounds for quacan c and we use as c initial point for quacan a fraction of the previous solution c do 15 i=1,n lint(i)=dmax1(-delta, lint(i)) 15 uint(i)=dmin1(delta, uint(i)) c c do 13 i=1,n 13 s(i)=s(i)*delta/snor0 c Set the convergence criterion for quacan at not-the-first calls within c the current iteration epsq=acumenor*gpnor else c c The initial point for quacan at its first call within a particular c iteration is the null increment. However, if iguess=1 we call c the user-provided subroutine guess to compute the inital approxima- c tion to the solution of the bound-constrained quadratic subproblem. c if(iguess.eq.1) then call guess(n, x, grad, ro, ra, lambda, mu, lint, uint, s) do 26 i=1,n if(s(i).lt.lint(i)) s(i) = lint(i) if(s(i).gt.uint(i)) s(i) = uint(i) 26 continue else do 19 i=1,n 19 s(i)=0.d0 endif endif