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----------------------------------------------------------------------