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