subroutine ivmhess(n,m,x,y,arow,acol,aval,anze,b,epsdsn, + epsden,epsnlsn,maxoutit,maxinnit,maxfc,p, + extrap,hess,output,kprint,eta,theta,tau,tol, + lmin,lmax,gamma,sterel,steabs,epsmac,infty, + tau1,tau2,f,dsupn,deucn2,amax,outit,innit,fcnt, + gcnt,hcnt,modh,flag) C Subroutine IVM implements the Inexact Variable Metric method C to find the local minimizers of a given function with convex C constraints (linear inequality constraints in this implementation), C described in: C C R. Andreani, E. G. Birgin, J. M. Martínez and J. Yuan, "Spectral C projected gradient and variable metric methods for optimization C with linear inequalities", IMA Journal of Numerical Analysis 25, C pp. 221-252, 2005. C C The user must supply the external subroutines evalf, evalg, C and evalH to evaluate the objective function, its gradient and C its hessian, respectively. Depending on parameter hess, evalH C subroutine may be an empty subroutine. C C This version 22 SET 2003 by E.G.Birgin, J.M.Martinez and J.Yuan. C Final revision ... by E.G.Birgin, J.M.Martinez and J. Yuan. C C ON ENTRY: C C n integer, C number of variables, C C m integer, C number of linear inequality constraints, C C x double precision x(n), C initial estimation of the solution, C C y double precision y(m), C initial estimation of the Lagrange multipliers, C C lda integer C leading dimension of matrix A (see below), C C A double precision A(m,n), C linear inequalities coefficients matrix, C C b double precision b(m), C linear inequalities right hand side, C C epsdsn double precision, C tolerance for the stopping criterion || amax d ||_sup < epsdsn, C where d is the IVM search direction and amax is C the maximum step which preserves feasibility, C C epsden double precision, C tolerance for the stopping criterion || amax d ||_2 < epsden, C where d is the IVM search direction and amax is C the maximum step which preserves feasibility, C C epsnlsn double precision, C tolerance for the stopping criterion of the inner C solver based on the sup-norm of the projected gradient C of the Lagrangian function, C C maxoutit integer, C maximum number of outer iterations, C C maxinnit integer, C maximum number of inner iterations, C C maxfc integer, C maximum number of function evaluations, C C p integer, C number of previous function values to be considered C in the nonmonotone line search, C C RECOMMENDED = 10 C C extrap logical, C if extrap = TRUE then, when solving the subproblem (using C an active set strategy), extrapolations are tried to increase C the dimension of the active set (it may be recomended for high C dimensional subproblems which comes from original problems with C many constraints); else extrapolations are not used at all. C C RECOMMENDED = .true. C C hess logical, C if hess = TRUE then the true Hessian is used in the C quadratic subproblems (the Identity is used if the function C is linear); else the spectral approximation is used. C C RECOMMENDED = try both! C C output integer, C < 0 means no print at all; C = 0 just initial and final information; C = 1 information at each outer iteration; C = 2 initial and final information for the inner C jobs: inner subproblems and line searches; C = 3 information of the inner iterations and C the line search iterations; C = 4 detailed information of the inner iterations. C C kprint integer, C maximum number of components to be printed when an C array must be showed, C C RECOMMENDED = 5 C C eta double precision, C eta in (0,1] controls the inexactness for solving the subproblem. C eta = 1 means "solve the subproblem exactly". C C RECOMMENDED: 0.9 C C theta double precision, C theta in (eta,1) is used to warranty that the iterates are C interior enough. C C RECOMMENDED: 0.95 C C tau double precision, C tau > 0 is used to compute the search direction of a C subproblem of type q(y) = 0.5 y^t B y + b^t y as C (B + sigma I) dk = - \nabla q(yk). See formulae (44) and C the end of Section 7 for details. C C RECOMMENDED = 10^-6 C C tol double precision, C tol in (0,1) is used, i the algorithm to solve the subproblems C (which uses an active set strategy) to decide if the current C face must be abandoned or not. Near 0 means to abandon the C face as fast as possible. Near 1 means try to find the solution C inside the current face. C C RECOMMENDED: 0.9 C C lmin double precision, C lmin > 0, lower bound for the safeguarded spectral C steplenght, C C RECOMMENDED: 10^-3 C C lmax double precision, C lmax >= lmin, upper bound for the safeguarded spectral C steplenght, C C RECOMMENDED: 10^3 C C gamma double precision, C gamma > 0, constant for the Armijo-like stooping criterion C of the line search, C C RECOMMENDED: 10^-4 C C sterel double precision, C epsrel > 0, a relative small number used for finite difference C approximation of derivatives, C C RECOMMENDED: 10^-7 C C steabs double precision, C epsabs > 0, an absolute small number used for finite difference C approximation of derivatives, C C RECOMMENDED: 10^-10 C C epsmac double precision, C epsmac > 0, machine epsilon, i.e., the smallest positive value C such that 1 + epsmac > 1. C C RECOMMENDED: In IEEE double precision floating point C arithmetic, the machine epsilon is something near 10^-16. C But the recomendation is to fix its value running the C following subroutine: C C subroutine mcheps(epsmac) C double precision epsmac C epsmac = 1.0d0 C 20 if ( ( 1.0d0 + epsmac ) .ne. 1.0d0 ) then C epsmac = epsmac / 2.0d0 C go to 20 C end if C epsmac = 2.0d0 * epsmac C end C C infty double precision, C a big number, C C RECOMMENDED = 10^20 C C The next two parameters are parameters for subroutine modchl C which computes the "Eskow-Schnabel" modified Cholesky factorization C of a given symetric matrix, taken from: C C ALGORITHM 695, COLLECTED ALGORITHMS FROM ACM. C PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 17, NO. 3, SEPTEMBER, 1991, PP. 306-312. C C tau1 double precision, C tau1 > 0, tolerance used for determining when to switch phase 2 C C RECOMMENDED: 10^-6 (approx. epsmac ^ 1/3, where epsmac C is the machine epsilon (approx. 10^-16). C C tau2 double precision, C tolerance used for determining the maximum condition C number of the final 2X2 submatrix. C C RECOMMENDED: 10^-6 (approx. epsmac ^ 1/3, where epsmac C tau2 > 0, is the machine epsilon (approx. 10^-16). C C ON RETURN: C C x double precision x(n), C final estimation of to the local minimizer, C C y double precision y(m), C final estimation of the Lagrange multipliers, C C f double precision, C function value at the approximation to the local C minimizer, C C dsupn double precision, C || d ||_sup at the final iteration, C where d is the IVM search direction, C C deucn2 double precision, C || d ||_2^2 at the final iteration, C where d is the IVM search direction, C C amax double precision, C maximum step to preserve feasibility in the last C computed direction, C C outit integer, C number of outer iterations, C C innit integer, C number of inner iterations, C C fcnt integer, C number of function evaluations, C C gcnt integer, C number of gradient evaluations, C C hcnt integer, C number of Hessian evaluations, C C flag integer, C termination parameter: C 0= convergence with sup-norm of the IVM search direction, C 1= convergence with euclidian-norm of the IVM search direction, C 6= too many outer iterations, C 7= too many inner iterations, C 8= too many functional evaluations, C -1= error in evalf subroutine, C -2= error in evalg subroutine, C -3= error in evalH subroutine, C -4= error in cholcol subroutine, C -5= error in sscol subroutine. C PARAMETERS integer mmax parameter (mmax = 100 000) integer nmax parameter (nmax = 4 000) integer pmax parameter (pmax = 100) C SCALAR ARGUMENTS logical hess,extrap integer anze,fcnt,flag,gcnt,hcnt,innit,kprint,m,maxfc, + maxinnit,maxoutit,modh,n,outit,output,p double precision amax,dsupn,deucn2,epsdsn,epsden,epsnlsn,f,eta, + theta,tau,tol,lmin,lmax,gamma,tau1,tau2,sterel,steabs, + epsmac,infty C ARRAY ARGUMENTS integer arow(anze),acol(anze) double precision aval(anze),b(m),x(n),y(m) C LOCAL SCALARS integer i,j,k,inform,nind,nprint,mprint,adyind double precision fbest,fnew,gtd,redc,L,Q,nLcsupn,nLpsupn,num, + den,alpha,atemp,fmax,adymax,adymaxi,maxelem,lambda,sts, + sty,gtg,tsmall,Ltmp,Qtmp,Lnew,Qnew,adytmp,adynew,amaxtmp, + amaxnew,redctmp,redcnew,nLpi,w1max,w2max,nLmax,dymax, + atamax,Hsafemax,sigma C LOCAL ARRAYS integer ind(mmax),pe(nmax),arini(mmax+1) double precision C(nmax,nmax),H(nmax,nmax),Hsafe(nmax,nmax), + E(nmax),xnew(nmax),xtmp(nmax),xbest(nmax),g(nmax), + gnew(nmax),gtmp(nmax),d(nmax),dnew(nmax),dtmp(nmax), + ynew(mmax),ytmp(mmax),ybest(mmax),nL(mmax),nLnew(mmax), + nLtmp(mmax),nLc(mmax),dy(mmax),lastfv(0:pmax-1), + bbar(mmax),w1(nmax),w2(nmax),w3(nmax),w4(nmax),w5(mmax) C EXTERNAL SUBROUTINES external modchl,solve,cholcol,sscol,evallnlhess C INTRINSIC FUNCTIONS intrinsic abs,max,min,mod,sqrt C =============================================================== C INEXACT VARIABLE METRIC METHOD C =============================================================== C =============================================================== C INITALIZATION C =============================================================== C PRINT PRESENTATION INFORMATION if ( output .ge. 0 ) then write(*, fmt=9010) write(10,fmt=9010) end if C CONSTANTS nprint = min(kprint,n) mprint = min(kprint,m) C COUNTERS outit = 0 innit = 0 fcnt = 0 gcnt = 0 hcnt = 0 modh = 0 C SET WHERE EACH LINE OF SPARSE MATRIX A STARTS do i = 1,m arini(i) = 0 end do arini(m+1) = anze + 1 do i = 1,anze if ( arini(arow(i)) .eq. 0 ) then arini(arow(i)) = i end if end do C INITIALIZE LAST FUNCTIONAL VALUES VECTOR do i = 0,p - 1 lastfv(i) = - infty end do C EVALUATE THE FUNCTION call evalf(n,x,f,inform) fcnt = fcnt + 1 if ( inform .ne. 0 ) then ! ERROR IN EVALF SUBROUTINE, STOP flag = -1 go to 900 end if C EVALUATE THE GRADIENT call evalg(n,x,g,inform) gcnt = gcnt + 1 if ( inform .ne. 0 ) then ! ERROR IN EVALG SUBROUTINE, STOP flag = -2 go to 900 end if C STORE FUNCTION VALUE FOR THE NONMONOTONE LINE SEARCH lastfv(0) = f C INITIALIZE BEST SOLUTION, BEST FUNCTIONAL VALUE AND ITS C ASSOCIATED LAGRANGE MULTIPLIERS fbest = f do i = 1,n xbest(i) = x(i) end do do i = 1,m ybest(i) = y(i) end do C DEFINE INITIAL SPECTRAL STEPLENGHT tsmall = 0.0d0 do i = 1,n tsmall = max( tsmall, abs( x(i) ) ) end do tsmall = max( sterel * tsmall, steabs ) do i = 1,n xtmp(i) = x(i) - tsmall * g(i) end do call evalg(n,xtmp,gtmp,inform) gcnt = gcnt + 1 gtg = 0.0d0 do i = 1,n gtg = gtg + g(i) ** 2 end do sts = tsmall ** 2 * gtg sty = 0.0d0 do i = 1,n sty = sty + g(i) * ( gtmp(i) - g(i) ) end do sty = - tsmall * sty if ( sty .le. 0.0d0 ) then lambda = lmax else lambda = min( lmax, max( lmin, sts / sty ) ) end if C =============================================================== C MAIN OUTER LOOP C =============================================================== 100 continue outit = outit + 1 C =============================================================== C STEP 1: COMPUTE THE SEARCH DIRECTION C =============================================================== C --------------------------------------------------------------- C Inner subproblem C --------------------------------------------------------------- C To compute the search direction we want to C C Minimize (approx) Qk(d) = 0.5 d^t Hk d + gk^t d C subject to A(xk + d) <= b. C Equivalently (via duality), we want to C C Minimize (approx) -Lk(y) = 0.5 (gk + A^t y) Hk^-1 (gk + A^t y) + C y^t (b - A xk) C subject to y >= 0. C The primal solution is then obtained by: C C d = - Hk^-1 (gk + A^t y) C --------------------------------------------------------------- C Initialization C --------------------------------------------------------------- if ( output .ge. 2 ) then write(*, fmt=9200) write(10,fmt=9200) end if C Compute bbar = b - A xk do i = 1,m bbar(i) = b(i) end do do i = 1,anze bbar(arow(i)) = bbar(arow(i)) - aval(i) * x(acol(i)) end do c Just to avoid bbar(i) = -10^-16 due to rounding errors do i = 1,m bbar(i) = max( 0.0d0, bbar(i) ) end do C Evaluate the Hessian, if ( hess ) then call evalH(n,x,nmax,H,inform) hcnt = hcnt + 1 if ( inform .ne. 0 ) then ! Error in evalh subroutine, stop flag = -3 go to 900 end if C Verify if the function is a linear function maxelem = 0.0d0 do j = 1,n do i = 1,n maxelem = max( maxelem, abs( H(i,j) ) ) end do end do if ( maxelem .le. epsmac ) then do j = 1,n do i = 1,n H(i,j) = 0.0d0 end do H(j,j) = 1.0d0 end do if ( output .ge. 2 ) then write(*, fmt=9201) write(10,fmt=9201) end if else if ( output .ge. 2 ) then write(*, fmt=9202) write(10,fmt=9202) end if end if C or use the spectral approximation else do j = 1,n do i = 1,n H(i,j) = 0.0d0 end do H(j,j) = 1.0d0 / lambda end do if ( output .ge. 2 ) then write(*, fmt=9203) lambda write(10,fmt=9203) lambda end if end if C Save the Hessian do j = 1,n do i = 1,n Hsafe(i,j) = H(i,j) end do end do C Modified Cholesky factorization of Hk (Eskow-Schanbel) call modchl(nmax,n,H,w1,tau1,tau2,pe,E) C COUNT THE NUMBER OF TIMES THE HESSIAN IS ACTUALLY MODIFIED if ( E(n) .gt. 0.0d0 ) then modh = modh + 1 if ( hess ) then if (output .ge. 2 ) then write(*, 9204) E(n) write(10,9204) E(n) end if end if end if C Add E to the saved Hessian do i = 1,n Hsafe(pe(i),pe(i)) = Hsafe(pe(i),pe(i)) + E(i) end do C --------------------------------------------------------------- C Given y, this subroutine computes: C C (a) The corresponding primal point: d = - Hk^-1 (gk + A^t y); C C (b) The maximum step to mantain feasibility of xk + d in the C primal feasible set A( xk + alpha d ) <= b; C C (c) The reduction constant to warranty that xk + redc * d is C an (enough) interior point; C C (d) Lk(d,y) and Qk(redc d); C C (e) The gradient of -L as a function of y. C --------------------------------------------------------------- call evallnlhess(m,n,nmax,y,g,arow,acol,aval,anze,bbar,H,pe,d, + amax,redc,L,Q,nL,w1,w5,w3,w4,theta,infty) C --------------------------------------------------------------- C Inner main loop C --------------------------------------------------------------- 200 continue C Compute chopped antigradient and projected antigradient norms; C and store in nind the number of free variables and in array nind C its identifiers. nind = 0 nLcsupn = 0.0d0 nLpsupn = 0.0d0 do i = 1,m if ( y(i) .gt. 0.0d0 ) then nind = nind + 1 ind(i) = nind nLc(i) = 0.0d0 nLpi = - nL(i) else ind(i) = 0 nLc(i) = max( 0.0d0, - nL(i) ) nLpi = max( 0.0d0, - nL(i) ) end if nLcsupn = max( nLcsupn, abs( nLc(i) ) ) nLpsupn = max( nLpsupn, abs( nLpi ) ) end do C Print inner iteration information if ( output .ge. 3 ) then if ( output .ge. 4 ) then write(*, fmt=9210) innit,mprint,m,(y(i),i=1,mprint) write(*, fmt=9220) nind,m write(*, fmt=9240) -L,nLcsupn,nLpsupn write(*, fmt=9250) nprint,n,(d(i),i=1,nprint) write(*, fmt=9260) amax,redc,Q write(10,fmt=9210) innit,mprint,m,(y(i),i=1,mprint) write(10,fmt=9220) nind,m write(10,fmt=9240) -L,nLcsupn,nLpsupn write(10,fmt=9250) nprint,n,(d(i),i=1,nprint) write(10,fmt=9260) amax,redc,Q else write(*, fmt=9205) innit,-L,nLpsupn,nind write(10,fmt=9205) innit,-L,nLpsupn,nind end if end if C --------------------------------------------------------------- C Inner stooping criteria C --------------------------------------------------------------- C Testing precognized inner stooping criterion: C Qk(redc d) <= eta L(d,y) if ( Q .le. eta * L ) then ! Precongnized stooping criterion satisfied, ! stop the inner algorithm if ( output .ge. 2 ) then write(*, fmt=9270) write(10,fmt=9270) end if go to 300 end if C Testing projected-gradient norm inner stooping criterion if ( nLpsupn .le. epsnlsn ) then ! Stationary point of the subproblem, ! stop the inner algorithm if ( output .ge. 2 ) then write(*, fmt=9280) write(10,fmt=9280) end if go to 300 end if C Testing number of inner iterations if ( innit .ge. maxinnit ) then ! Maximum number of inner iterations exceeded, ! stop the algorithm flag = 7 go to 900 end if C --------------------------------------------------------------- C Doing an inner iteration C --------------------------------------------------------------- innit = innit + 1 C --------------------------------------------------------------- C Testing whether we leave the current face or not and compute C the corresponding direction C --------------------------------------------------------------- if ( nLpsupn * tol .le. nLcsupn ) then C ----------------------------------------------------------- C Leaving the face in the chopped gradient direcion C ----------------------------------------------------------- if ( output .ge. 4 ) then write(*, fmt=9290) write(10,fmt=9290) end if do i = 1,m dy(i) = nLc(i) end do else C ----------------------------------------------------------- C Interior to the face iteration (all the operations of this C job are done into the reduced space) C ----------------------------------------------------------- if ( output .ge. 4 ) then write(*, fmt=9300) write(10,fmt=9300) end if C Compute direction dy = - (A Hk^-1 A^t + sigma I)^-1 nL C in two steps: C C First, solving (sigma Hk + A^t A) r = - A^t nL; and C second dy = - ( 1 / sigma ) ( nL + A r ). C Compute C = sigma Hk + A^t A C Compute A^t A do j = 1,n do i = j,n C(i,j) = 0.0d0 end do end do do k = 1,m if ( ind(k) .ne. 0 ) then do j = arini(k),arini(k+1)-1 do i = j,arini(k+1)-1 C(acol(i),acol(j)) = + C(acol(i),acol(j)) + aval(i) * aval(j) end do end do end if end do atamax = 0.0d0 do j = 1,n do i = j,n atamax = max( atamax, abs( C(i,j) ) ) end do end do C Compute sigma Hsafemax = 0.0d0 do j = 1,n do i = j,n Hsafemax = max( Hsafemax, abs( Hsafe(i,j) ) ) end do end do sigma = tau * atamax / Hsafemax C Add sigma Hk do j = 1,n do i = j,n C(i,j) = C(i,j) + sigma * Hsafe(i,j) C(j,i) = C(i,j) end do end do C Cholesky factorization of C call cholcol(nmax,n,C,epsmac,inform) if ( inform .ne. 0 ) then flag = -4 go to 900 end if C Solve system C w2 = - A^t nL C Compute w1 = - A^t nL do j = 1,n w1(j) = 0.0d0 end do do i = 1,anze if ( ind(arow(i)) .ne. 0 ) then w1(acol(i)) = w1(acol(i)) - aval(i) * nL(arow(i)) end if end do C Solve system C w2 = w1 call sscol(nmax,n,C,w1,w2,epsmac,inform) if ( inform .ne. 0 ) then flag = -5 go to 900 end if C Compute dy = - ( 1 / sigma ) ( nL + A w2 ) do i = 1,m if (ind(i) .ne. 0 ) then dy(i) = nL(i) else dy(i) = 0.0d0 end if end do do i = 1,anze if ( ind(arow(i)) .ne. 0 ) then dy(arow(i)) = dy(arow(i)) + aval(i) * w2(acol(i)) end if end do do i = 1,m dy(i) = - dy(i) / sigma end do end if if ( output .ge. 4 ) then write(*, fmt=9310) mprint,(dy(i),i=1,mprint) write(10,fmt=9310) mprint,(dy(i),i=1,mprint) end if C --------------------------------------------------------------- C Compute maximum step C --------------------------------------------------------------- adymax = infty adyind = 0 do i = 1,m if ( dy(i) .lt. 0.0d0 ) then adymaxi = - y(i) / dy(i) if ( adymaxi .le. adymax ) then adymax = adymaxi adyind = i end if end if end do if ( output .ge. 4 ) then if ( adyind .ne. 0 ) then write(*, fmt=9320) adymax,adyind,y(adyind) write(10,fmt=9320) adymax,adyind,y(adyind) else write(*, fmt=9321) adymax write(10,fmt=9321) adymax end if end if C --------------------------------------------------------------- C Compute the one-dimensional minimizer of the quadratic -Lk(d,y) C along the dy direction (as the primal is feasible the dual is C bounded and, so, this minimizer exists). C --------------------------------------------------------------- C alphady = - nL^t dy / ( dy^t A Hk^-1 A^t dy ) C Compute w2 = Hk^-1 A^t dy C w1 = A^t dy do j = 1,n w1(j) = 0.0d0 end do do i = 1,anze w1(acol(i)) = w1(acol(i)) + aval(i) * dy(arow(i)) end do C Solve system Hk w2 = w1 call solve(nmax,n,H,pe,w2,w1,w3,w4) C Compute the one-dimensional minimizer C num = nL^t dy nLmax = 0.0d0 dymax = 0.0d0 do i = 1,m nLmax = max( nLmax, abs( nL(i) ) ) dymax = max( dymax, abs( dy(i) ) ) end do if ( nLmax .eq. 0.0d0 .or. dymax .eq. 0.0d0 ) then num = 0.0d0 else num = 0.0d0 do i = 1,m num = num + ( nL(i) / nLmax ) * ( dy(i) / dymax ) end do num = num * nLmax * dymax end if C den = v^t r w1max = 0.0d0 w2max = 0.0d0 do i = 1,n w1max = max( w1max, abs( w1(i) ) ) w2max = max( w2max, abs( w2(i) ) ) end do if ( w1max .eq. 0.0d0 .or. w2max .eq. 0.0d0 ) then den = 0.0d0 else den = 0.0d0 do i = 1,n den = den + ( w1(i) / w1max ) * ( w2(i) / w2max ) end do den = den * w1max * w2max end if C One-dimensional minimizer adynew = - num / den if ( output .ge. 4 ) then write(*, fmt=9340) adynew write(10,fmt=9340) adynew c write(*, fmt=9341) adynew,num,den c write(10,fmt=9341) adynew,num,den end if adynew = min( adynew, adymax ) C --------------------------------------------------------------- C Compute the first candidate to be the new iterate (without C considering the extrapolation) C --------------------------------------------------------------- do i = 1,m ynew(i) = y(i) + adynew * dy(i) end do if ( adynew .eq. adymax .and. adyind .ne. 0 ) then ynew(adyind) = 0.0d0 end if do i = 1,m if ( ynew(i) .le. epsmac ) then ynew(i) = 0.0d0 end if end do C --------------------------------------------------------------- C Given y, this subroutine computes: C C (a) The corresponding primal point: d = - Hk^-1 (gk + A^t y); C C (b) The maximum step to mantain feasibility of xk + d in the C primal feasible set A( xk + alpha d ) <= b; C C (c) The reduction constant to warranty that xk + redc * d is C an (enough) interior point; C C (d) Lk(d,y) and Qk(redc d); C C (e) The gradient of -L as a function of y. C --------------------------------------------------------------- call evallnlhess(m,n,nmax,ynew,g,arow,acol,aval,anze,bbar,H, + pe,dnew,amaxnew,redcnew,Lnew,Qnew,nLnew,w1,w5,w3, + w4,theta,infty) C --------------------------------------------------------------- C Extrapolation C --------------------------------------------------------------- if ( extrap ) then 250 continue adytmp = 2.0d0 * adynew do i = 1,m ytmp(i) = max( 0.0d0, y(i) + adytmp * dy(i) ) end do call evallnlhess(m,n,nmax,ytmp,g,arow,acol,aval,anze,bbar, + H,pe,dtmp,amaxtmp,redctmp,Ltmp,Qtmp,nLtmp,w1, + w5,w3,w4,theta,infty) if ( -Ltmp .le. -Lnew ) then adynew = adytmp do i = 1,m ynew(i) = ytmp(i) nLnew(i) = nLtmp(i) end do do i = 1,n dnew(i) = dtmp(i) end do amaxnew = amaxtmp redcnew = redctmp Lnew = Ltmp Qnew = Qtmp go to 250 end if end if C --------------------------------------------------------------- C Set next iterate C --------------------------------------------------------------- do i = 1,m y(i) = ynew(i) nL(i) = nLnew(i) end do do i = 1,n d(i) = dnew(i) end do amax = amaxnew redc = redcnew L = Lnew Q = Qnew C --------------------------------------------------------------- C End of the inner iteration C --------------------------------------------------------------- go to 200 C --------------------------------------------------------------- C End of the inner main loop C --------------------------------------------------------------- C =============================================================== C PRINT OUTER ITERATION INFORMATION AND TEST STOPPING CRITERIA C =============================================================== 300 continue C --------------------------------------------------------------- C Compute d norms C --------------------------------------------------------------- dsupn = 0.0d0 deucn2 = 0.0d0 do i = 1,n dsupn = max( dsupn, abs( d(i) ) ) deucn2 = deucn2 + d(i) ** 2 end do C --------------------------------------------------------------- C Print outer iteration information C --------------------------------------------------------------- if ( output .ge. 1 ) then write(*, fmt=9020) outit,f,dsupn,amax write(10,fmt=9020) outit,f,dsupn,amax end if C --------------------------------------------------------------- C Stooping criterion C --------------------------------------------------------------- if ( abs( redc ) * dsupn .le. epsdsn ) then ! IVM SEARCH DIRECTION SUP-NORM STOPPING CRITERION SATISFIED, STOP flag = 0 go to 900 end if if ( redc ** 2 * deucn2 .le. epsden ** 2 ) then ! IVM SEARCH DIRECTION 2-NORM STOPPING CRITERION SATISFIED, STOP flag = 1 go to 900 end if if ( outit .gt. maxoutit ) then ! MAXIMUM NUMBER OF OUTER ITERATIONS EXCEEDED, STOP flag = 6 go to 900 end if C =============================================================== C STEP 2: NONMONOTONE LINE SEARCH C =============================================================== C --------------------------------------------------------------- C Line search initialization C --------------------------------------------------------------- if ( output .ge. 2 ) then write(*, fmt=9500) write(10,fmt=9500) end if C Compute g^t d gtd = 0.0d0 do i = 1,n gtd = gtd + g(i) * d(i) end do C Compute the maximum functional value of the last p iterations fmax = lastfv(0) do i = 1,p - 1 fmax = max(fmax,lastfv(i)) end do C Compute first trial alpha = redc do i = 1,n xnew(i) = x(i) + alpha * d(i) end do call evalf(n,xnew,fnew,inform) fcnt = fcnt + 1 if (inform .ne. 0) then ! Error in evalf subroutine, stop flag = -1 go to 900 end if C --------------------------------------------------------------- C Line search main loop C --------------------------------------------------------------- 400 continue C Print iteration information if ( output .ge. 3 ) then write(*, fmt=9510) alpha,fnew write(10,fmt=9510) alpha,fnew end if C --------------------------------------------------------------- C Line search stooping criteria C --------------------------------------------------------------- if ( fnew .le. fmax + gamma * alpha * gtd ) then ! Nonmonotone Armijo-like stopping criterion satisfied, ! stop the line search if ( output .ge. 2 ) then write(*, fmt=9520) write(10,fmt=9520) end if go to 500 end if if ( fcnt .ge. maxfc ) then ! Maximum number of functional evaluations exceeded, ! stop the algorithm flag = 8 go to 900 end if C --------------------------------------------------------------- C Line search iteration C --------------------------------------------------------------- C Safeguarded quadratic interpolation if ( alpha .le. 0.1d0 ) then alpha = 0.5d0 * alpha else atemp = ( - gtd * alpha ** 2 ) / + ( 2.0d0 * ( fnew - f - alpha * gtd ) ) if ( atemp .lt. 0.1d0 .or. atemp .gt. 0.9d0 * alpha ) then atemp = 0.5d0 * alpha end if alpha = atemp end if C Compute trial point do i = 1,n xnew(i) = x(i) + alpha * d(i) end do C Evaluate the function call evalf(n,xnew,fnew,inform) fcnt = fcnt + 1 if (inform .ne. 0) then ! Error in evalf subroutine, stop flag = -1 go to 900 end if C --------------------------------------------------------------- C End of the line search iteration C --------------------------------------------------------------- go to 400 C --------------------------------------------------------------- C End of the line search main loop C --------------------------------------------------------------- C =============================================================== C PREPARING FOR THE NEXT OUTER ITERATION C =============================================================== 500 continue C EVALUATE THE GRADIENT AT THE NEW ITERATE call evalg(n,xnew,gnew,inform) gcnt = gcnt + 1 if ( inform .ne. 0 ) then ! ERROR IN EVALG SUBROUTINE, STOP flag = -2 go to 900 end if C COMPUTE THE SPECTRAL STEPLENGHT sts = 0.0d0 sty = 0.0d0 do i = 1,n sts = sts + ( xnew(i) - x(i) ) ** 2 sty = sty + ( xnew(i) - x(i) ) * ( gnew(i) - g(i) ) end do if ( sty .le. 0.0d0 ) then lambda = lmax else lambda = min( lmax, max( lmin, sts / sty ) ) end if C SET THE NEW CURRENT POINT AND ITS GRADIENT do i = 1,n x(i) = xnew(i) g(i) = gnew(i) end do C SET NEW FUNCTION VALUE AND SAVE IT FOR THE NONMONOTONE C LINE SEARCH f = fnew lastfv(mod(outit,p)) = f C COMPARE THE NEW FUNCTION VALUE AGAINST THE BEST FUNCTION C VALUE AND, IF SMALLER, UPDATE THE BEST FUNCTION C VALUE, THE CORRESPONDING BEST POINT AND ITS ASSOCIATED C LAGRANGE MULTIPLIERS if ( f .lt. fbest ) then fbest = f do i = 1,n xbest(i) = x(i) end do do i = 1,m ybest(i) = y(i) end do end if C =============================================================== C END OF THE OUTER ITERATION C =============================================================== go to 100 C =============================================================== C END OF THE MAIN LOOP C =============================================================== C =============================================================== C RETURN C =============================================================== 900 continue C SET X, Y AND F WITH THE BEST SOLUTION, ITS ASSOCIATED LAGRANGE C MULTIPLIERS AND ITS FUNCTIONAL VALUE f = fbest do i = 1,n x(i) = xbest(i) end do do i = 1,m y(i) = ybest(i) end do C WRITE OUTPUT STATUS if ( output .ge. 0 ) then if ( flag .eq. 0 ) then write(*, fmt=9600) flag,epsdsn write(10,fmt=9600) flag,epsdsn else if ( flag .eq. 1 ) then write(*, fmt=9610) flag,epsden write(10,fmt=9610) flag,epsden else if ( flag .eq. 6 ) then write(*, fmt=9620) flag,maxoutit write(10,fmt=9620) flag,maxoutit else if ( flag .eq. 7 ) then write(*, fmt=9630) flag,maxinnit write(10,fmt=9630) flag,maxinnit else if ( flag .eq. 8 ) then write(*, fmt=9640) flag,maxfc write(10,fmt=9640) flag,maxfc else if ( flag .eq. -1 ) then write(*, fmt=9650) flag write(10,fmt=9650) flag else if ( flag .eq. -2 ) then write(*, fmt=9660) flag write(10,fmt=9660) flag else if ( flag .eq. -3 ) then write(*, fmt=9670) flag write(10,fmt=9670) flag else if ( flag .eq. -4 ) then write(*, fmt=9680) flag write(10,fmt=9680) flag else if ( flag .eq. -5 ) then write(*, fmt=9690) flag write(10,fmt=9690) flag end if end if C WRITE STATISTICS if ( output .ge. 0 ) then write(*, fmt=9700) f,dsupn,sqrt(deucn2),amax,flag write(*, fmt=9710) outit,fcnt,gcnt,hcnt,modh,innit write(10,fmt=9700) f,dsupn,sqrt(deucn2),amax,flag write(10,fmt=9710) outit,fcnt,gcnt,hcnt,modh,innit end if return C NON-EXECUTABLE STATEMENTS 9010 format(/,1X,'INEXACT VARIABLE METRIC METHOD') 9020 format(/,1X,'OUTIT= ',I10,' F= ',1P,D14.7,' DSUPN= ',1P,D14.7, + ' AMAX= ',1P,D14.7) 9200 format(/,6X,'INNER SOLVER (COMPUTING THE SEARCH DIRECTION FOR ', + 'THE NEXT ITERATION)') 9201 format( 6X,'(The function is linear, using Gk = I)') 9202 format( 6X,'(Using Gk = the true Hessian)') 9203 format( 6X,'(Using Gk = (1 / lambda) I, where lambda = ', + 1P,D14.7,' is the spectral', + /,6X,'steplenght)') 9204 format( 6X,'(',1P,D14.7,' was added to the diagonal)') 9205 format( 6X,'INNIT = ',I6,' -L =',1P,D14.7,' NLPSUPN = ', + 1P,D14.7,' FACDIM = ',I6) 9210 format(/,6X,'INNER ITERATION: ',I6, + /,6X,'Current y (first ',I6,' components over a total of ', + I6,'):', + /,1(6X,6(1X,1PD11.4))) 9220 format( 6X,'Number of free variables at this point: ',I6, + ' (over a total of ',I6, ')') 9230 format( 6X,'Indices of the free variables (first ',I6,'): ', + /,1(6X,10(1X,I6))) 9240 format( 6X,'Functional value = ',1PD14.7, + /,6X,'Chopped gradient sup-norm = ',1PD14.7, + /,6X,'Projected gradient sup-norm = ',1PD14.7) 9250 format( 6X,'Associated primal point (first ',I6,' components ', + 'over a total of ',I6,'):', + /,1(6X,6(1X,1PD11.4))) 9260 format( 6X,'Maximum primal step = ',1PD14.7, + /,6X,'Including constant of reduction = ',1PD14.7, + /,6X,'Q(redc d) = ',1PD14.7) 9270 format( 6X,'Precognized stooping criterion satisfied') 9280 format( 6X,'Projected gradient norm stooping criterion ', + 'satisfied') 9290 format( 6X,'Leaving the face') 9300 format( 6X,'Remaining into the face') 9310 format( 6X,'dy direction (first ',I6, ' components): ', + /,1(6X,6(1X,1PD11.4))) 9320 format( 6X,'Maximum step: ',1PD14.7,' Variable index: ',I6, + ' value: ',1PD14.7) 9321 format( 6X,'There is no maximum step. Maximum step: ',1PD14.7) 9330 format( 6X,'The one-dimensional minimizer does not exist') 9340 format( 6X,'One-dimensional minimizer: ',1PD14.7) 9341 format( 6X,'One-dimensional minimizer: ',1PD14.7,1PD14.7,1PD14.7) 9500 format(/,6X,'LINE SEARCH') 9510 format( 6X,'Step = ',1P,D14.7,' F = ',1P,D14.7) 9520 format( 6X,'Armijo-like stooping criterion satisfied') 9600 format(/,1X,'IVM terminated with flag ',I2,': ', + 'convergence with the sup-norm of the search ', + /,1X,'direction smaller than ',1P,D14.7, + ' (parameter epsdsn)') 9610 format(/,1X,'IVM terminated with flag ',I2,': ', + 'convergence with the euclidian-norm of the ', + 'search', + /,1X,'direction smaller than ',1P,D14.7, + ' (parameter epsden)') 9620 format(/,1X,'IVM terminated with flag ',I2,': ', + 'too many outer iterations ', + /,1X,'(more than parameter maxoutit = ',I6,')') 9630 format(/,1X,'IVM terminated with flag ',I2,': ', + 'too many inner iterations ', + /,1X,'(more than parameter maxinnit = ',I6,')') 9640 format(/,1X,'IVM terminated with flag ',I2,': ', + 'too many functional evaluations ', + /,1X,'(more than parameter maxfc = ',I6,')') 9650 format(/,1X,'IVM terminated with flag ',I2,': ', + 'Error in evalf subroutine') 9660 format(/,1X,'IVM terminated with flag ',I2,': ', + 'Error in evalg subroutine') 9670 format(/,1X,'IVM terminated with flag ',I2,': ', + 'Error in evalH subroutine') 9680 format(/,1X,'IVM terminated with flag ',I2,': ', + 'Error in cholcol subroutine') 9690 format(/,1X,'IVM terminated with flag ',I2,': ', + 'Error in sscol subroutine') 9700 format(/,1X,'F = ', 1P,D17.10, + /,1X,'DSUPNORM = ',1X,1P,D16.10, + /,1X,'DEUCNORM = ',1X,1P,D16.10, + /,1X,'AMAX = ',1X,1P,D16.10, + /,1X,'FLAG = ',15X,I2) 9710 format(/,1X,'OUTIT = ',1X,I16, + /,1X,'FCNT = ',1X,I16, + /,1X,'GCNT = ',1X,I16, + /,1X,'HCNT = ',1X,I16, + /,1X,'MODH = ',1X,I16, + /,1X,'INNIT = ',1X,I16) end c ****************************************************** c ****************************************************** subroutine cholcol(lda,n,A,epsmac,flag) integer lda,n,flag double precision A(lda,n),epsmac integer i,j,k flag = 0 do k = 1,n if ( A(k,k) .le. epsmac ) then flag = -1 return end if A(k,k) = sqrt( A(k,k) ) do i = k+1,n A(i,k) = A(i,k) / A(k,k) end do do j = k+1,n do i = j,n A(i,j) = A(i,j) - A(i,k) * A(j,k) end do end do end do return end c ****************************************************** c ****************************************************** subroutine sscol(lda,n,A,b,x,epsmac,flag) integer lda,n,flag double precision A(lda,n),b(n),x(n),epsmac integer i,j flag = 0 do i = 1,n x(i) = b(i) end do C Solving G y = b do j = 1,n if ( A(j,j) .le. epsmac ) then flag = -1 return end if x(j) = x(j) / A(j,j) do i = j+1,n x(i) = x(i) - x(j) * A(i,j) end do end do C Solving G^t x = y do i = n,1,-1 do j = i+1,n x(i) = x(i) - A(j,i) * x(j) end do x(i) = x(i) / A(i,i) if ( A(i,i) .le. epsmac ) then flag = -1 return end if end do return end c ****************************************************** c ****************************************************** subroutine evallnlhess(m,n,ldh,y,g,arow,acol,aval,anze,bbar,H,pe, + d,amax,redc,L,Q,nL,w1,w2,w3,w4,theta,infty) C SCALAR ARGUMENTS integer anze,m,n,ldh double precision amax,redc,L,Q,theta,infty C ARRAY ARGUMENTS integer arow(anze),acol(anze),pe(n) double precision y(m),g(n),aval(anze),H(ldh,n),d(n), + bbar(m),nL(m),w1(n),w2(m),w3(n),w4(n) C LOCAL SCALARS integer i,j double precision tmp,dmax,gmax,w1max,ymax,w2mbbarmax C --------------------------------------------------------------- C Computing corresponding primal point d = - Hk^-1 ( gk + A^t y ) C --------------------------------------------------------------- C w1 = - ( gk + A^t y ) do j = 1,n w1(j) = - g(j) end do do i = 1,anze w1(acol(i)) = w1(acol(i)) - aval(i) * y(arow(i)) end do C Solve system Hk d = w1 call solve(ldh,n,H,pe,d,w1,w3,w4) C --------------------------------------------------------------- C Computing maximum step to mantain feasibility of xk + d C in the primal feasible set A( xk + alpha d ) <= b C --------------------------------------------------------------- C w2 = A d do i = 1,m w2(i) = 0.0d0 end do do i = 1,anze w2(arow(i)) = w2(arow(i)) + aval(i) * d(acol(i)) end do C amax = min_{i | [A d]_i > 0} { [bbar]_i / [A d]_i } amax = infty do i = 1,m if ( w2(i) .gt. 0.0d0 ) then amax = min( amax, bbar(i) / w2(i) ) end if end do C --------------------------------------------------------------- C Computing the reduction constant to warranty that xk + redc * d C is an (enough) interior point C --------------------------------------------------------------- redc = min( 1.0d0, theta * amax ) C --------------------------------------------------------------- C Computing Lk(d,y) and Qk(redc d) C --------------------------------------------------------------- dmax = 0.0d0 gmax = 0.0d0 w1max = 0.0d0 do i = 1,n dmax = max( dmax, abs( d(i) ) ) gmax = max( gmax, abs( d(i) ) ) w1max = max( w1max, abs( w1(i) ) ) end do ymax = 0.0d0 w2mbbarmax = 0.0d0 do i = 1,m ymax = max( ymax, abs( y(i) ) ) w2mbbarmax = max( w2mbbarmax, abs( w2(i) - bbar(i) ) ) end do dmax = 1.0d0 gmax = 1.0d0 w1max = 1.0d0 ymax = 1.0d0 w2mbbarmax = 1.0d0 C Computing 0.5 d^t Hk d = - d^t (gk + A^t y) = d^t w1 if ( dmax .eq. 0.0d0 .or. w1max .eq. 0.0d0 ) then tmp = 0.0d0 else tmp = 0.0d0 do i = 1,n tmp = tmp + ( d(i) / dmax ) * ( w1(i) / w1max ) end do tmp = tmp * dmax * w1max end if L = 0.5d0 * tmp Q = 0.5d0 * redc ** 2 * tmp C Computing gk^t d if ( gmax .eq. 0.0d0 .or. dmax .eq. 0.0d0 ) then tmp = 0.0d0 else tmp = 0.0d0 do i = 1,n tmp = tmp + ( g(i) /gmax ) * ( d(i) / dmax ) end do tmp = tmp * gmax * dmax end if L = L + tmp Q = Q + redc * tmp C Computing y^t ( A xk - b + A d ) = y^t ( w2 - bbar ) if ( ymax .eq. 0.0d0 .or. w2mbbarmax .eq. 0.0d0 ) then tmp = 0.0d0 else tmp = 0.0d0 do i = 1,m tmp = tmp + ( y(i) / ymax ) * ( w2(i) - bbar(i) ) / w2mbbarmax end do tmp = tmp * ymax * w2mbbarmax end if L = L + tmp C --------------------------------------------------------------- C Computing - \nabla L(d,y) = A Hk^-1 (gk + A^t y) + bbar = bbar - A d C --------------------------------------------------------------- C nL is the antigradient of L. As we are minimizing -L, nL is the C positive gradient of the function being minimized. do i = 1,m nL(i) = bbar(i) - w2(i) end do end