c
c Solve, approximately, the quadratic subproblem
c
call quacan (mapro, grad, n, s, fq, gq, lint, uint, eta, factor,
* kon, mvp, ipqua, mitqu, maxmv, dont, kdont, epsq, ier,
* gc, bd, xn, gn, d, ingc, infree, hess, inhess,
* m, ro, reg, ispar, nelem, inda, a, ichan, irep0)
c
iquatot=iquatot+kon
mvptot=mvptot+mvp
irep = irep + irep0
c
c Now we test whether the quadratic model produced enough decrease,
c that is, we test whether fq is small enough.
c For this purpose, we consider the ``easy-quadratic model'', and we
c compute d, the minimizer of the easy-quadratic model in the trust region
c box. If the quadratic model did not produce a decrease at least
c equal to some fraction of the decrease produced by the easy-quadratic
c model, the solution of the quadratic model is replaced by the
c one of the easy-quadratic model
c
flimo = 0.d0
dtd = 0.d0
do 24 i = 1, n
d(i)=-grad(i)
if(d(i).lt.lint(i)) d(i) = lint(i)
if(d(i).gt.uint(i)) d(i) = uint(i)
dtd = dtd + d(i)**2
24 flimo = flimo + d(i) * grad(i)
flimo = flimo + dtd/2.d0
if(fq .gt. 1.d-6 * flimo) then
nrepl = nrepl + 1
do 25 i=1,n
25 s(i)=d(i)
fq = flimo
endif
c
c The vector s contains our trial increment. Now we test whether it
c can be accepted.
c xn will contain the trial point
c
do 9 i=1,n
xn(i)=x(i)+s(i)
c
c These are safeguards to ensure that the trial point do not
c escape from the box by rounding errors, which should be
c catastrophic.
c The quantity 1.d-14 is meant to be something a litter bit larger
c than the machine precision. You must change it if the machine
c precision changes. Here it is assumed to be 1.d-16
c
if(xn(i).gt.u(i) - 1.d-14 * dabs(u(i)))xn(i)=u(i)
if(xn(i).lt.l(i) + 1.d-14 * dabs(l(i)))xn(i)=l(i)
9 continue
c
c Compute the functional value at the trial point
call phi (n, nonlc, xn, mu, ra, fn, fxn, gn, hres,
* hess, inhess, 1)
if(m.gt.0) call pena (n, xn, fn, gn, res,
* m, ila, lambda, ro, ispar, nelem, inda, a, b, bd, 1)
naf = naf + 1
if(ipr.eq.1)then
write(10, *)' Trust region radius:', delta
write(*, *)' Trust region radius:', delta
write(10,*)' Predicted reduction:', -fq
write(10,*)' Actual reduction:', flagra - fn
write(*,*)' Predicted reduction:', -fq
write(10,*)' Reason for termination of quacan, ier =',ier
write(*,*)' Reason for termination of quacan, ier =',ier
write(10,*)' Number of iterations at last call to quacan:',kon
write(*,*)' Number of iterations at last call to quacan:',kon
write(*,*)' Actual reduction:', flagra - fn
write(10,*)' Objective function value at trial point:',fn
write(*,*)' Objective function value at trial point:',fn
endif
c
c Test if the trial point can be accepted
c
desce = dmax1 (relarm * fq, - absarm)
if ( fn. le . flagra + desce ) then
c
c The trial point is accepted and is going to be the new iterate
c
do 10 i=1,n
10 x(i)=xn(i)
pred = -fq
ared = flagra - fn
flagra=fn
fx=fxn
kont=kont+1
goto3
c
endif
c
c The trial point was not good enough. Reduce delta
c
c
c First, verify if delta is very small (strangled point)
c
if(delta.le.epsd)then
istop=4
if(ipr.ge.0)then
write(*,*)
write(10,*)
write(*,*)' BOX stops by excessively small trust region'
write(*,*)' (strangled point)'
write(10,*)' BOX stops by excessively small trust region'
write(10,*)' (strangled point)'
call printibox(kont,n10,x,grad,flagra,gpnor,gpcont,
* naf,iquatot,mvptot,delta, pred, ared, nrepl, irep)
endif
return
endif
c Updating delta after a failure of the sufficient descent
c condition
c
snor0=0.d0
do 11 i=1,n
11 snor0= dmax1(snor0, dabs(s(i)))
c
c contrac is a contraction factor for the trust region radius.
c If fn was less than flagra - 0.1 fq
c contrac is set to be equal to 2. (Therefore, the trust-region
c radius is not excessively reduced.)
c If fn was greater than flagra - 0.1 fq we set contrac equal to 10.
c That is, is fn was much larger than flagra we reduce the trust-region
c radius rather drastically.
c
if(fn.le.flagra - 0.1d0 * fq) then
contrac = 2.d0
else
contrac = 10.d0
endif
delta = snor0 / contrac
line=1
goto12
c
end
subroutine printibox(kon,n10,x,grad,f,gpnor,gpcont,
* naf,iquatot,mvp,del, pred, ared, nrepl, irep)
implicit double precision(a-h,o-z)
dimension x(n10), grad(n10)
write(10,*)
write(10,*)' BOX iteration ', kon
write(10,*)' First (at most 10) components of current point:'
write(10,*)(x(i), i=1,n10)
write(10,*)' First (at most 10) components of the gradient:'
write(10,*)(grad(i), i=1,n10)
write(10,*)' Objective function value:',f
write(*,*)
write(*,*)' BOX iteration ', kon
write(*,*)' First (at most 10) components of current point:'
write(*,*)(x(i), i=1,n10)
write(*,*)' First (at most 10) components of the gradient:'
write(*,*)(grad(i), i=1,n10)
write(*,*)' Objective function value:',f
if(kon.gt.0) then
write(10,*)' Predicted reduction:', pred
write(10,*)' Actual reduction:', ared
write(*,*)' Predicted reduction:', pred
write(*,*)' Actual reduction:', ared
endif
write(10,*)' 2-Norm of the continuous projected gradient :',
* gpcont
write(10,*)' 2-Norm of the (non-continuous) projected gradient:'
* , gpnor
write(10,*)' Functional evaluations:', naf
write(10,*)' Total number of quacan iterations:', iquatot
write(10,*)' Number of matrix-vector products:', mvp
write(10,*)' Number of times in which the solution of the'
write(10,*)' quadratic model was replaced by the solution of'
write(10,*)' the easy-quadratic model:', nrepl
write(10,*)' Number of times in which the guessed increment'
write(10,*)' or, perhaps, the backtracked increment'
write(10,*)' was worse than the null increment in terms of the'
write(10,*)' quadratic model:', irep
write(*,*)' 2-Norm of the continuous projected gradient :',
* gpcont
write(*,*)' 2-Norm of the (non-continuous) projected gradient:'
* , gpnor
write(*,*)' Functional evaluations:', naf
write(*,*)' Total number of quacan iterations:', iquatot
write(*,*)' Number of matrix-vector products:', mvp
write(*,*)' Number of times in which the solution of the'
write(*,*)' quadratic model was replaced by the solution of'
write(*,*)' the easy-quadratic model:', nrepl
write(*,*)' Number of times in which the guessed increment'
write(*,*)' or, perhaps, the backtracked increment'
write(*,*)' was worse than the null increment in terms of the'
write(*,*)' quadratic model:', irep
if(kon.ne.0)then
write(10,*)' The current point was obtained within a trust'
write(10,*)' region of radius', del
write(*,*)' The current point was obtained within a trust'
write(*,*)' region of radius', del
endif
write(10,*)
write(*,*)
return
end
C
C----------------------------------------------------------------------
subroutine quacan (mapro, b, n, x, f, g, l, u, eta, factor,
* kon, mvp, ipri, maxit, maxef, dont, kdont, eps, ier,
* gc, bd, xn, gn, d, ingc, infree, hess, inhess,
* mrowsa, ro, reg, ispar, nelem, inda, amat, ichan, irep0)
c
c Initiated in November 1995
c
c This subroutine aims to minimize a quadratic function
c
c f(x) = (1/2) x^T ( B + ro A^T A + reg I ) x + b^T x
c
c (f:\R^n \to \R) on the box given by l \leq x \leq u.
c It is especially suited for n large.
c The matrix A has mrowsa rows and n columns and can be stored dense
c or sparse according to the parameter ispar.
c Quacan replaces a previous version, avoiding
c the cumbersome driving of indices that characterizes active set
c strategies.
c The main purpose was to develop an ``easy code'', easy to manipu-
c late and modify for testing different ideas.
c The method implemented in Quacan is described in the report BFGMR
c (paper by Bielschowky, Friedlander, Gomes, Martinez and Raydan).
c It is based on a mild active set strategy that uses conjugate
c gradients inside the faces, projected searches and chopped
c gradients to leave the faces. The criterion to leave a face is
c
c ||g_c|| \geq eta ||g_p|| (1)
c
c where g_c is the chopped gradient, g_p is the projected gradient,
c defined as in BFGMR and || . || is the euclidian norm.
c
c Parameters:
c Mapro: subroutine that must be declared external in the calling
c routine or program. The calling sequence of mapro is
c call mapro(n, nind, ind, u, v)
c This subroutine must be coded by the user, taking into account that
c n is the number of variables of the problem and that v must be
c the product B u. Moreover, you must assume, when you code mapro,
c that only nind components of u are nonnull and that ind is the set
c of indices of those components. In other words, you must write
c mapro in such a way that v is the vector whose i-th entry is
c \Sum_{j=1}^{nind} B_{i ind(j)} u_ind(j)
c However, observe that you must assume that, in u, the whole vector
c is present, with its n components, even the zeroes. So, if you
c decide to code mapro without taking into account the presence
c of ind and nind, you can do it. A final advertence concerning
c your coding of mapro. You must test, in the code, if nind = n.
c In this case, do not use ind, because it is a void vector and,
c of course, has no utility. In other words, if nind=n you must
c not assume that ind is the vector (1, 2, ..., n). It is not. A
c final observation: probably, if nind is close to n, it is not
c worthwhile to use ind, due to the cost of accessing the correct
c indices. If you want, you can test, within your mapro, if
c (say) nind > n/2, and, in this case use a straightforward scalar
c product for the components of v.
c
c----------------------------------------------------------------------
c INPUT PARAMETERS THAT DEFINE RO AND THE MATRIX A
c
c mrowsa, ro, ispar, nelem, inda, amat
c
c
c mrowsa : number of rows of the matrix A. Set mrowsa = 0 if there is no
c Lagrangian term at all
c
c ro : real*8 number given in the definition of the problem
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. Set nelem = 1 if ispar is not 1.
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 amat : 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 de inda is (i, j).
c
c
c ichan : Input integer parameter. If you set ichan = 1, quacan
c tests if the initial approximation is better than the
c null initial approximation. If it is not, quacan changes
c the given initial approximation by the vector 0. In
c other words, if the quadratic value at the given initial
c approximation is greater than 0, the initial approxima-
c tion is replaced by the null vector. If you set ichan
c different from 1, the initial approximation given is
c not changed, no matter its functional value.
c
c irep0: Output parameter that can take the values 0 or 1. If
c irep0 = 1, this means that, being ichan=1, the initial
c point of quacan was really changed to the null initial
c point.
c
c-------------------------------------------------------------------------
c PARAMETER RELATED TO THE REGULARIZATION TERM
C
C reg: real*8 regularization parameter.
c
c
c-------------------------------------------------------------------------
c
c n: number of variables.
c x (dimension n): on entry, initial approximation, feasible
c to the box, on exit, best point obtained.
c f: on exit, value of f(x).
c g: when convergence ocurs, this vector (dimension n) contains the
c final values of the gradient at x.
c eta: On input, number belonging to (0, 1) according to (1) above.
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.
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 eps: Small positive number for declaring convergence when the eucli-
c dian norm of the projected gradient is less than
c or equal to eps.
c maxit: maximum number of iterations allowed.
c maxef: maximum number of matrix-vector products allowed.
c kon: on exit, number of iterations performed.
c mvp: on exit, matrix-vector products (calls to mapro)
c dont: ``lack of enough progress'' measure. The algorithm stops by
c ``lack of enough progress'' when
c f(x_k) - f(x_{k+1}) \leq dont * max \{f(x_j)-f(x_{j+1}, j= than component',i,' of u'
write(10,*)' Component',i,' of l is >= than component',i,' of u'
endif
return
endif
8 continue
c
c Test if the initial point is within to the box
do 1 i = 1, n
if(x(i).lt.l(i).or.x(i).gt.u(i))then
ier=4
if(ipri.lt.0) return
write(*,*)
write(*,*)' The initial point is not feasible'
write(*,*)
write(10,*)
write(10,*)' The initial point is not feasible'
write(10,*)
return
endif
1 continue
c
indica=0
c indica=1 means that the iteration that begins now is of conjugate
c gradient type. indica=0 means that it is of chopped gradient type
c or internal gradient type.
kon=0
mvp=0
if(ipri.ge.0) n10=min0(n,10)
noprog=0
c
c Progress is the maximum progress obtained up to now
progress=0.d0
c
c Evaluate the gradient at the initial point
call grad (n, mapro, b, x, g, mvp, hess, inhess,
* mrowsa, ro, reg, ispar, nelem, inda, amat,
* bd(n+1), bd(n+mrowsa+1))
c Evaluate the objective function at the initial point
call fuqu(n, x, g, b, f)
c
if(ichan.eq.1) then
c Test if the initial point given is better than the null vector
c If it is not, it will be changed.
if(f. gt. 0.d0) then
c Replace the given initial point by the null vector
irep0 = 1
do 23 i = 1, n
23 x(i) = 0.d0
c
c Evaluate the gradient at the initial point
call grad (n, mapro, b, x, g, mvp, hess, inhess,
* mrowsa, ro, reg, ispar, nelem, inda, amat,
* bd(n+1), bd(n+mrowsa+1))
c Evaluate the objective function at the initial point
call fuqu(n, x, g, b, f)
c
endif
endif
c
c Here begins the iteration loop
c
c nfree is the number of components of the internal gradient
c (number of free variables)
c infree is the vector of indices of components of the internal
c gradient (free variables).
c ncgc is the number of components of the chopped gradient
c ingc is the vector of indices of the components of the chopped
c gradient.
c
10 nfree=0
ncgc=0
do 3 i=1,n
if(x(i).gt.l(i).and.x(i).lt.u(i))then
nfree=nfree+1
infree(nfree)=i
endif
if(x(i).eq.l(i).and.g(i).lt.0.d0)then
ncgc=ncgc+1
ingc(ncgc)=i
endif
if(x(i).eq.u(i).and.g(i).gt.0.d0)then
ncgc=ncgc+1
ingc(ncgc)=i
endif
3 continue
c
c
c Compute the squared 2-norm of the internal gradient
c However, if indica=1, this is not necessary since it has
c been already computed at the end of the previous iteration
if(indica.eq.0)then
ginor2=0.d0
if(nfree.gt.0)then
do 4 i=1,nfree
j=infree(i)
4 ginor2=ginor2+g(j)**2
endif
endif
c Compute the 2-norm of the chopped gradient
gcnor2=0.d0
if(ncgc.gt.0)then
do 22 i=1,ncgc
j=ingc(i)
22 gcnor2=gcnor2+g(j)**2
endif
c Compute the 2-norm of the projected gradient
gcnor=dsqrt(gcnor2)
gpnor2=ginor2+gcnor2
gpnor=dsqrt(gpnor2)
c
c Test convergence criterion
c
if(gpnor.le.eps)then
ier=0
if(ipri.ge.0)then
write(*,*)
write(10,*)
write(*,*)' Convergence detected at iteration ', kon
write(10,*)' Convergence detected at iteration ', kon
call printi(kon,n,n10,x,g,f,gpnor,gcnor,nfree,ncgc,mvp,ichop,
* internal, iconju)
write(10,*)' Use Quacan and you will be happy for ever!!'
write(*,*)' Use Quacan and you will be happy for ever!!'
c
write(*,*)
write(10,*)
endif
return
endif
c
c Test maximum number of iterations
c
if(kon.ge.maxit)then
ier=1
if(ipri.ge.0)then
write(*,*)
write(10,*)
write(*,*)' Maximum number of iterations ', kon,' reached'
write(10,*)' Maximum number of iterations ', kon,' reached'
call printi(kon,n,n10,x,g,f,gpnor,gcnor,nfree,ncgc,mvp,
* ichop,internal, iconju)
endif
return
endif
c
c Test maximum number of matrix-vector products
c
if(mvp.ge.maxef)then
ier=2
if(ipri.ge.0)then
write(*,*)
write(10,*)
write(*,*)' Maximum number of matrix-vector products ',
* maxef,' reached'
write(10,*)' Maximum number of matrix-vector products ',
* maxef,' reached'
call printi(kon,n,n10,x,g,f,gpnor,gcnor,nfree,ncgc,mvp,
* ichop,internal, iconju)
endif
return
endif
if(kon.ne.0)then
prog=dmax1(0.d0, fant-f)
progress=dmax1(prog, progress)
if(prog.le.dont*progress)then
noprog=noprog+1
if(noprog.ge.kdont)then
ier=3
if(ipri.ge.0)then
write(*,*)
write(10,*)
write(*,*)' Quacan stops because of lack of progress'
write(10,*)' Quacan stops because of lack of progress'
call printi(kon,n,n10,x,g,f,gpnor,gcnor,nfree,ncgc,mvp,
* ichop,internal, iconju)
endif
return
endif
else
noprog=0
endif
endif
fant=f
c
c Printing session
c
if(ipri.gt.0)then
if(mod(kon,ipri).eq.0.or.kon.eq.0)then
call printi(kon,n,n10,x,g,f,gpnor,gcnor,nfree,ncgc,mvp,
* ichop,internal, iconju)
endif
endif
c
c
c A new iteration is going to be performed
kon=kon+1
c
c Test if the current face must be leaved or not
c
if(gcnor.ge.eta*gpnor)then
c ichop, internal and iconju are to indicate the printer which
c type of iteration has been performed
if(ipri.ge.0)then
ichop=1
internal=0
iconju=0
endif
c
c The current face must be abandoned
c
c Use the chopped gradient to leave it
c
c
c gc is the (negative) chopped gradient vector
c
do 5 i=1,n
5 gc(i)=0.d0
do 6 i=1,ncgc
j=ingc(i)
gc(j)= -g(j)
6 continue
c
c Compute the minimum breakpoint alfabreak along the chopped
c gradient direction. (Observe that all the components of gc
c that appear dividing are different from zero.)
c
j=ingc(1)
if(x(j).eq.l(j)) then
alfabreak = (u(j) - l(j))/gc(j)
jbreak=j
else
alfabreak = (l(j) - u(j))/gc(j)
jbreak=j
endif
if(ncgc.gt.1)then
do 9 i = 2, ncgc
j=ingc(i)
if(x(j).eq.l(j)) then
alfa = (u(j) - l(j))/gc(j)
if(alfa.lt.alfabreak)then
alfabreak=alfa
jbreak=j
endif
else
alfa = (l(j) - u(j))/gc(j)
if(alfa.lt.alfabreak)then
alfabreak=alfa
jbreak=j
endif
endif
9 continue
endif
c
c
c
c Compute the product of the Hessian times the chopped gradient
call mapro(n, ncgc, ingc, gc, bd, hess, inhess)
if(mrowsa.gt.0)call roata(n, ncgc, ingc, gc, bd,
* mrowsa, ro, reg, ispar, nelem, inda, amat, bd(n+1),
* bd(n+mrowsa+1))
mvp=mvp+1
c
c Compute the denominator for the linesearch
c
den=0.d0
do 7 i=1,ncgc
c In the following line an overflow appeared on 28/12/95, running gedan
7 den=den + gc(ingc(i)) * bd(ingc(i))
c
c If den .le. 0, the next iteration will be the minimum breakpoint
c along the chopped gradient direction, or, perhaps, something
c better (after the double-stepping procedure)
c
if(den.le.0.d0) then
call double (n, b, x, f, g, gc, ingc, ncgc, l, u, mapro, factor,
* alfabreak, jbreak, xn, gn, bd, mvp, hess, inhess,
* mrowsa, ro, reg, ispar, nelem, inda, amat)
c
c The new iteration has been computed and stored in x
c The new gradient is stored in g
c The new objective function value is at f
c The iteration stops here
c The next iteration will not be conjugate
indica=0
goto10
endif
c Compute the numerator of the linesearch coefficient
c
alfa= gcnor2/den
if(alfa.ge.alfabreak) then
call double (n, b, x, f, g, gc, ingc, ncgc, l, u, mapro, factor,
* alfabreak, jbreak, xn, gn, bd, mvp, hess, inhess,
* mrowsa, ro, reg, ispar, nelem, inda, amat)
c
c The new iteration has been computed and stored in x
c The new gradient is stored in g
c The new objective function value is at f
c The iteration stops here
c The next iteration will not be conjugate (indica=0)
indica=0
goto10
endif
c
c If the minimizer along the direction of the chopped gradient is
c before the minimum breakpoint, we adopt this minimizer as new
c iterate
c
do 12 i=1,ncgc
j=ingc(i)
12 x(j)=x(j)+alfa*gc(j)
c
c Compute the gradient at the new point x
do 13 i=1,n
13 g(i)=g(i)+alfa*bd(i)
c
c Compute the new objective function value
call fuqu(n, x, g, b, f)
c The iteration stops here
c
c The next iteration will not be conjugate (indica=0)
indica=0
goto10
c Here finishes the ``if'' that begins with ``if(gcnor.ge.eta*gpnor)''
c
else
c Now, we consider the case where the face must not be abandoned
c that is, gcnor < eta*gpnor
c In this case, we perform conjugate gradient iterations
c However, we need an indicator for telling us if the iteration
c must be ``conjugate gradient'' or merely gradient. A conjugate
c gradient iteration will take place only if the indicator
c indica is equal to 1. Otherwise, this will be a gradient iteration
c along the (negative) internal gradient. Observe that we set indica=0 at
c the steps of the algorithm that corresponds to chopped gradient
c iterations.
c Put the search direction on vector d
if(indica.eq.0)then
c Gradient search direction
do 15 i=1,n
15 d(i)=0.d0
do 14 j=1,nfree
i=infree(j)
14 d(i)=-g(i)
endif
if(indica.eq.1)then
c Conjugate search direction
do 19 j=1,nfree
i=infree(j)
19 d(i)= -g(i) + beta*d(i)
endif
c
c Compute the minimum breakpoint along d
alfabreak=-1.d0
do 16 i=1,nfree
j=infree(i)
if(d(j).eq.0.d0)goto16
alfa = (u(j)-x(j))/d(j)
if(alfa.gt.0.d0)then
if(alfabreak.eq.-1.d0.or.alfa.lt.alfabreak)then
alfabreak=alfa
jbreak=j
goto16
endif
else
alfa = (l(j)-x(j))/d(j)
if(alfa.le.0.d0)goto16
if(alfabreak.eq.-1.d0.or.alfa.lt.alfabreak)then
alfabreak=alfa
jbreak=j
endif
endif
16 continue
c
c Compute the denominator of the linesearch
call mapro(n, nfree, infree, d, bd, hess, inhess)
if(mrowsa.gt.0)call roata(n, nfree, infree, d, bd,
* mrowsa, ro, reg, ispar, nelem, inda, amat, bd(n+1),
* bd(n+mrowsa+1))
mvp=mvp+1
den=0.d0
do 17 i=1,nfree
j=infree(i)
17 den=den+bd(j)*d(j)
c
c If the denominator of the linesearch is .le. 0, go to the breakpoint
c
if(den.le.0.d0)then
call double (n, b, x, f, g, d, infree, nfree,
* l, u, mapro, factor, alfabreak, jbreak, xn, gn, bd, mvp,
* hess, inhess,
* mrowsa, ro, reg, ispar, nelem, inda, amat )
c
c The new iteration has been computed and stored in x
c The new gradient is stored in g
c The new objective function value is at f
c The iteration stops here
indica=0
c ichop, internal and iconju are to indicate the printer which
c type of iteration has been performed
if(ipri.ge.0)then
ichop=0
internal=1
iconju=0
endif
goto10
endif
c Compute the linesearch coefficient
alfa= ginor2/den
if(alfa.ge.alfabreak) then
call double (n, b, x, f, g, d, infree, nfree, l, u,
* mapro, factor, alfabreak, jbreak, xn, gn, bd, mvp,
* hess, inhess, mrowsa, ro, reg, ispar, nelem, inda, amat)
c
c The new iteration has been computed and stored in x
c The new gradient is stored in g
c The new objective function value is at f
c The iteration stops here
c ichop, internal and iconju are to indicate the printer which
c type of iteration has been performed
if(ipri.ge.0)then
ichop=0
if(indica.eq.0)then
internal=1
iconju=0
else
iconju=1
internal=0
endif
endif
c The next iteration will not be conjugate (indica=0)
indica=0
goto10
else
c
c The new iterate is in the same face. The next iteration will
c be conjugate
do 18 i=1,nfree
j=infree(i)
18 x(j)=x(j)+alfa*d(j)
do 20 i=1,n
20 g(i)=g(i)+alfa*bd(i)
gant2=ginor2
ginor2=0.d0
do 21 j=1,nfree
i=infree(j)
21 ginor2=ginor2+g(i)**2
beta=ginor2/gant2
call fuqu(n, x, g, b, f)
if(ipri.ge.0)then
ichop=0
if(indica.eq.0)then
internal=1
iconju=0
else
iconju=1
internal=0
endif
endif
indica=1
goto10
endif
endif
end
c End of subroutine Quacan
subroutine grad (n, mapro, b, x, g, mvp, hess, inhess,
* mrowsa, ro, reg, ispar, nelem, inda, amat, u, v)
implicit double precision (a-h, o-z)
dimension x(n), b(n), g(n), nada(1)
dimension hess(*), inhess(*)
dimension inda(nelem,2), amat(*)
c This subroutine is called by quacan to compute the gradient
c of the quadratic f at the point x.
c So, the output vector g is
c g = (B + ro A^T A + reg I) x + b
c The product B x is computed using the user subroutine mapro,
c which is described in the comments of Quacan.
c mvp increments in 1 the number of matrix-vector products.
c
c
call mapro (n, n, nada, x, g, hess, inhess)
c
c
c
c Add to B x the term ro A^T A x
c
if(mrowsa.ne.0)then
call roata (n, n, nada, x, g,
* mrowsa, ro, reg, ispar, nelem, inda, amat, u, v)
endif
mvp=mvp+1
do 1 i=1,n
1 g(i)=g(i)+b(i)
return
end
c End of subroutine grad
subroutine printi(kon,n,n10,x,g,f,gpnor,gcnor,nfree,ncgc,mvp,
* ichop,internal, iconju)
implicit double precision(a-h,o-z)
dimension x(n10), g(n10)
write(10,*)
write(10,*)' Quacan-iteration ', kon
write(*,*)' Quacan-iteration ', kon
if(kon.gt.0)then
if(ichop.eq.1)then
write(10,*)' This point comes from a chopped gradient search'
write(*,*)' This point comes from a chopped gradient search'
endif
if(internal.eq.1)then
write(10,*)' This point comes from an internal gradient search'
write(*,*)' This point comes from an internal gradient search'
endif
if(iconju.eq.1)then
write(10,*)' This point comes from a conj. gradient search'
write(*,*)' This point comes from a conj. gradient search'
endif
endif
write(10,*)' First (at most 10) components of current point:'
write(10,*)(x(i), i=1,n10)
write(10,*)' First (at most 10) components of the gradient:'
write(10,*)(g(i), i=1,n10)
write(10,*)' Value of the quadratic objective function:',f
write(10,*)' 2-Norm of the projected gradient:', gpnor
write(10,*)' 2-Norm of the chopped gradient:', gcnor
write(10,*)' Number of free variables:', nfree
nbound=n-nfree
write(10,*)' Number of variables on their bounds:',nbound
write(10,*)' Number of components of chopped gradient:',ncgc
write(10,*)' Number of matrix-vector products:', mvp
write(10,*)
write(*,*)
write(*,*)' First (at most 10) components of current point:'
write(*,*)(x(i), i=1,n10)
write(*,*)' Value of the quadratic objective function:',f
write(*,*)' Norm of the projected gradient:', gpnor
write(*,*)' Norm of the chopped gradient:', gcnor
write(*,*)' Number of free variables:', nfree
write(*,*)' Number of variables on their bounds:',nbound
write(*,*)' Number of components of chopped gradient:',ncgc
write(*,*)' Number of matrix-vector products:', mvp
write(*,*)
return
end
subroutine fuqu(n, x, g, b, f)
implicit double precision (a-h,o-z)
dimension x(n), g(n), b(n)
c This subroutine evaluates
c f(x) = (1/2) x^T B x + b^T x
c using g(x) = B x + b
c
c We use the formula
c
c f(x) = (1/2) x^T [B x + b] + b^T x / 2 = (1/2) x^T g(x) + b^T x / 2
c = x^T (g(x) + b)/2
c
f=0.d0
do 7 i=1,n
7 f=f+ x(i) * (g(i) + b(i))
f=f/2.d0
return
end
subroutine double (n, b, x, f, g, d, infree, nfree, l, u,
* mapro, factor, alfabreak, jbreak, xn, gn, bd, mvp,
* hess, inhess, mrowsa, ro, reg, ispar, nelem, inda, amat)
implicit double precision (a-h, o-z)
double precision l(n)
external mapro
dimension x(n), d(n), u(n), g(n), b(n)
dimension infree(n)
dimension xn(n), gn(n), bd(*)
dimension hess(*), inhess(*)
dimension inda(nelem,2), amat(*)
c
c Compute the trial point corresponding to the minimum break
c We use jbreak to put variables that, by rounding errors,
c are free, on the corresponding bounds. For example, the va-
c riable alfabreak certainly must be on a bound, but the variables
c that, perhaps, have a distance to a bound that is less than or
c equal to the distance of x(jbreak) to its bound, will be also
c put on the corresponding bound.
c
c
c Moreover, in dic 4 1995 we decided to put a variable on its bound
c if the difference between this bound and the variable is less
c than 1.d-14 (relative error). This modification is present in
c the version box9512, but not in box9511d. We assume a machine
c precision around 1.d-16 If the machine precision changes, you
c must change 1.d-14 consequently.
do 7 i=1,n
7 xn(i)=x(i)
do 1 i=1,nfree
j=infree(i)
xn(j)=x(j) + alfabreak*d(j)
1 continue
tole = dmin1(u(jbreak)-xn(jbreak), xn(jbreak)-l(jbreak))
tole = dmax1(0.d0, tole)
tole = 1.5d0*tole
c The following two lines should not be necessary
c if we neglected rounding errors. See comment above.
do 10 i=1,nfree
j=infree(i)
tolo=dmax1(tole, 1.d-14 * dabs(u(j)))
if(xn(j).ge.u(j)-tolo) xn(j)=u(j)
tolo=dmax1(tole, 1.d-14 * dabs(l(j)))
if(xn(j).le.l(j)+tolo) xn(j)=l(j)
10 continue
c
c Compute the gradient at this xn
c
do 2 i=1,n
2 gn(i)=g(i)+alfabreak*bd(i)
c
c Compute the functional value at this xn
c
call fuqu(n, xn, gn, b, fn)
c
c Copy xn in x, since this new point is necessarily better than
c the old one
c
do 3 i=1,nfree
j=infree(i)
3 x(j)=xn(j)
alfa=alfabreak
do 4 i=1,n
4 g(i)=gn(i)
f=fn
c
c Initiate de double-step procedure, with the aim of obtaining
c a point even better
c
9 alfa=factor*alfa
do 8 i=1,nfree
j=infree(i)
xn(j)=x(j)+alfa*d(j)
if(xn(j).gt.u(j))xn(j)=u(j)
if(xn(j).lt.l(j))xn(j)=l(j)
8 continue
call grad(n, mapro, b, xn, gn, mvp, hess, inhess,
* mrowsa, ro, reg, ispar, nelem, inda, amat,
* bd(n+1), bd(n+mrowsa+1))
call fuqu(n, xn, gn, b, fn)
c If the doubling step did not improve, return
if(fn.ge.f)return
f=fn
do 5 i=1,nfree
j=infree(i)
5 x(j)=xn(j)
do 6 i=1,n
6 g(i)=gn(i)
goto9
end
subroutine roata (n, nfree, infree, d, bd,
* m, ro, reg, ispar, nelem, inda, a, u, v)
implicit double precision (a-h, o-z)
c
c This subroutine adds the term ro A^T A d + reg I d to the vector bd
c See the comments on the use of mapro, and the comments of the
c calling subroutines for more details on storage of A, etc
c
dimension d(n), bd(*), inda(nelem,2), a(*), u(*), v(*)
dimension infree(*)
if(ispar.eq.0.and.nfree.eq.n)then
c
c Compute u = A d
c
do 1 i=1,m
u(i)=0.d0
do 1 j=1,n
k=(j-1)*m+i
c
c a(k) is [A]_ij
c
1 u(i)=u(i)+a(k)*d(j)
c
c Compute A^T u, pre-multiply by ro and add to bd
c
6 do 2 j=1,n
z=0.d0
do 3 i=1,m
k=(j-1)*m+i
c
c a(k) is [A]_ij
c
3 z=z+a(k)*u(i)
2 bd(j)=bd(j) + ro*z
if(reg.ne.0.d0) then
do 14 i=1,n
14 bd(i)=bd(i)+reg*d(i)
endif
return
endif
if(ispar.eq.0.and.nfree.lt.n)then
c
c Compute u = A d
c
do 5 i=1,m
u(i)=0.d0
do 5 k=1,nfree
j=infree(k)
kk=(j-1)*m+i
c
c a(kk) is [A]_ij
c
5 u(i)=u(i)+ a(kk)*d(j)
c
c Compute A^T u, pre-multiply by ro and add to bd
c
do 7 j=1,n
z=0.d0
do 8 i=1,m
k=(j-1)*m+i
c
c a(k) is [A]_ij
c
8 z=z+a(k)*u(i)
7 bd(j)=bd(j) + ro*z
if(reg.ne.0.d0) then
do 15 i=1,n
15 bd(i)=bd(i)+reg*d(i)
endif
return
endif
c
if(ispar.eq.1)then
c
c Compute u = A d
c
do 13 i=1,m
13 u(i)=0.d0
do 9 k=1,nelem
i=inda(k,1)
j=inda(k,2)
c
c a(k) is [A]_ij
c
9 u(i)=u(i)+a(k)*d(j)
c
c Compute A^T u, pre-multiply by ro and add to bd
c
do 10 i=1,n
10 v(i)=0.d0
do 11 k=1,nelem
i=inda(k,1)
j=inda(k,2)
c
c a(k) is [A]_ij
c
11 v(j)=v(j) + a(k)*u(i)
do 12 j=1,n
12 bd(j)=bd(j)+ro*v(j)
if(reg.ne.0.d0) then
do 16 i=1,n
16 bd(i)=bd(i)+reg*d(i)
endif
return
endif
c
c if(ispar.eq.1.and.nfree.lt.n) then ????
c I do not know how to take advantage of sparsity both in A and
c d for obtaining the product A d
c
end
subroutine pena (n, x, f, g, res,
* m, ila, lambda, ro, ispar, nelem, inda, a, b, aux, modo)
implicit double precision (a-h,o-z)
double precision lambda(*)
dimension x(n), g(n), inda(nelem, 2), a(*), b(*), res(*)
dimension aux(*)
c
c This subroutine adds to f the term
c
c lambda^T (A x - b) + (ro/2) || A x - b ||^2
c
c if modo = 1 and adds to g the term
c
c A^T [ lambda + ro (A x - b) ]
c
c if modo = 2 .
c
c We store A x - b in the vector res and we use the fact that, in
c BOX a call of pena with modo=2 is always immediately
c preceeded by a call with
c modo=1 so that we can use the values of res computed in modo=1
c for the computations for modo=2 .
c
c See the comments of BOX for the storage of A
c
c
c Compute the residual vector res
c
if(modo.eq.1)then
if(m.gt.0)then
if(ispar.eq.0)then
do 1 i=1,m
res(i) = - b(i)
do 1 j=1,n
k=(j-1)*m+i
c
c a(k) is [A]_ij
c
1 res(i)=res(i)+a(k)*x(j)
else
do 2 i = 1, m
2 res(i) = - b(i)
do 3 k = 1, nelem
i = inda (k, 1)
j = inda (k, 2)
c
c a(k) is [A]_ij
c
3 res(i) = res(i) + a(k) * x(j)
endif
pes=0.d0
if(ila.eq.1)then
do 4 i=1,m
4 pes=pes+lambda(i)*res(i)
endif
z=0.d0
do 5 i=1,m
5 z=z+res(i)**2
f = f + pes + ro * z / 2.d0
endif
return
else
c
c Now, modo = 2
c
if(m.ne.0)then
do 7 i=1,m
7 aux(n+i) = lambda(i) + ro*res(i)
c
do 6 j=1,n
6 aux(j)=0.d0
if(ispar.eq.0)then
do 10 j=1,n
do 10 i=1,m
k=(j-1)*m+i
10 aux(j)=aux(j)+ a(k)*aux(n+i)
else
do 8 k=1,nelem
i=inda(k,1)
j=inda(k,2)
8 aux(j)=aux(j)+a(k)*aux(n+i)
endif
do 9 j=1,n
9 g(j)=g(j)+aux(j)
endif
return
endif
end
c
C----------------------------------------------------------------------
C
integer function mult( p, q)
C
Integer p, q, p0, p1, q0, q1
C
p1 = p/10000
p0 = mod(p,10000)
q1 = q/10000
q0 = mod(q,10000)
mult = mod( mod( p0*q1+p1*q0,10000)*10000+p0*q0,100000000)
return
end
C
C----------------------------------------------------------------------
C
real*8 function rnd(sem)
C
integer sem, mult
C
sem = mod( mult( sem, 3141581) + 1, 100000000)
rnd = sem/100000000.0d0
return
end
C
C----------------------------------------------------------------------