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