subroutine ivmispg(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, + 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 + tau 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 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,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 COMMON SCALARS double precision lambda C COMMON ARRAYS double precision g(nmax),bbar(mmax) C LOCAL SCALARS integer i,inform,nprint,mprint double precision fbest,fnew,gtd,redc,L,Q,nLpsupn,alpha,atemp,fmax, + sts,sty,gtg,tsmall,nLpi C LOCAL ARRAYS double precision xnew(nmax),xtmp(nmax),xbest(nmax),gnew(nmax), + gtmp(nmax),d(nmax),ybest(mmax),nL(mmax),lastfv(0:pmax-1), + w1(nmax),w3(nmax),w4(nmax),w5(mmax) C GENCAN VARIABLES logical Gnearlyq integer Gmaxitnfp,Gmaxitngp,Gmaxit,Gmaxfc,Gcgmaxit,Gcgscre, + Gmaxitnqmp,Ggtype,Ghtvtype,Gtrtype,Giprint,Giter,Gfcnt, + Ggcnt,Gcgcnt,Gspgiter,Gspgfcnt,Gtniter,Gtnfcnt,Gtnstpcnt, + Gtnintcnt,Gtnexgcnt,Gtnexbcnt,Gtnintfe,Gtnexgfe,Gtnexbfe, + Gflag,Gw2(mmax),Gmininterp,Gncomp double precision Gepsgpen,Gepsgpsn,Gepsnfp,Gfmin,Gdelta0,Gcggpnf, + Gcgepsi,Gcgepsf,Gepsnqmp,Gf,Gg(mmax),Ggpeucn2,Ggpsupn, + Gw11(mmax),Gw12(mmax),Gw13(mmax),Gw14(5*mmax),Geta, + Gdelmin,Glammax,Glammin,Gtheta,Ggamma,Gbeta,Gsigma1, + Gsigma2,Gsterel,Gsteabs,Gepsrel,Gepsabs,Ginfty,Gnint, + Gnext,Gl(mmax),Gu(mmax) C COMMON BLOCKS common /ldata2/ g,bbar,lambda C EXTERNAL SUBROUTINES external evallnlispg,gencan 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 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 C ---> In the case of the ISPG, Hk = (1 / lambda) I 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 C ---> In the case of the ISPG, C C Minimize (approx) -Lk(y) = 0.5 lambda || (gk + A^t y) ||_2^2 + 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 ---> In the case of the ISPG, C C d = - lambda (gk + A^t y) C --------------------------------------------------------------- C Initialization C --------------------------------------------------------------- if ( output .ge. 2 ) then write(*, fmt=9200) write(*, fmt=9203) lambda write(10,fmt=9200) write(10,fmt=9203) lambda 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 --------------------------------------------------------------- C Inner main loop C --------------------------------------------------------------- 200 continue 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 evallnlispg(m,n,y,g,arow,acol,aval,anze,bbar,lambda,d,amax, + redc,L,Q,nL,w1,w5,w3,w4,theta,infty) 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. nLpsupn = 0.0d0 do i = 1,m if ( y(i) .gt. 0.0d0 ) then nLpi = - nL(i) else nLpi = max( 0.0d0, - nL(i) ) end if 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=9240) -L,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=9240) -L,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 write(10,fmt=9205) innit,-L,nLpsupn 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 --------------------------------------------------------------- c innit = innit + 1 C GENCAN ITERATION, STARTING FROM y TO MINIMZE C C -L(y) = 0.5 lambda || (gk + A^t y) ||_2^2 + bbar^t y C C SUBJECT TO y >= 0. do i = 1,m Gl(i) = 0.0d0 Gu(i) = infty end do Gepsgpen = 0.0d0 Gepsgpsn = epsnlsn Gepsnfp = 0.0d0 Gmaxitnfp = 10000 Gmaxitngp = 10000 Gfmin = -1.0d+99 Gmaxit = 100 Gmaxfc = 10000 Gdelta0 = -1.0d0 Gcgmaxit = -1 Gcgscre = 2 Gcggpnf = Gepsgpsn Gcgepsi = 1.0d-1 Gcgepsf = 1.0d-5 Gepsnqmp = 1.0d-4 Gmaxitnqmp = 5 Gnearlyq = .true. Ggtype = 0 Ghtvtype = 0 Gtrtype = 0 Giprint = -1 Geta = 0.9d0 Gdelmin = 0.1d0 Glammax = 1.0d+03 Glammin = 1.0d-03 Gtheta = 1.0d-06 Ggamma = 1.0d-04 Gbeta = 0.5d0 Gsigma1 = 0.1d0 Gsigma2 = 0.9d0 Gnint = 2 Gnext = 2 Gmininterp = 4 Gncomp = 5 Gsterel = 1.0d-07 Gsteabs = 1.0d-10 Gepsrel = 1.0d-10 Gepsabs = 1.0d-20 Ginfty = 1.0d+20 call gencan(m,y,Gl,Gu,Gepsgpen,Gepsgpsn,Gepsnfp,Gmaxitnfp, +Gmaxitngp,Gfmin,Gmaxit,Gmaxfc,Gdelta0,Gcgmaxit,Gcgscre,Gcggpnf, +Gcgepsi,Gcgepsf,Gepsnqmp,Gmaxitnqmp,Gnearlyq,Ggtype,Ghtvtype, +Gtrtype,Giprint,Gf,Gg,Ggpeucn2,Ggpsupn,Giter,Gfcnt,Ggcnt,Gcgcnt, +Gspgiter,Gspgfcnt,Gtniter,Gtnfcnt,Gtnstpcnt,Gtnintcnt,Gtnexgcnt, +Gtnexbcnt,Gtnintfe,Gtnexgfe,Gtnexbfe,Gflag,Gw11,Gw12, +Gw13,Gw2(1),Gw14,Geta,Gdelmin,Glammax,Glammin,Gtheta, +Ggamma,Gbeta,Gsigma1,Gsigma2,Gnint,Gnext,Gmininterp,Gncomp, +Gsterel,Gsteabs,Gepsrel,Gepsabs,Ginfty) do i = 1,m if ( y(i) .le. epsmac ) then y(i) = 0.0d0 end if end do innit = innit + Giter 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 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,innit write(10,fmt=9700) f,dsupn,sqrt(deucn2),amax,flag write(10,fmt=9710) outit,fcnt,gcnt,hcnt,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)') 9205 format( 6X,'INNIT = ',I6,' -L =',1P,D14.7,' NLPSUPN = ', + 1P,D14.7) 9210 format(/,6X,'INNER ITERATION: ',I6, + /,6X,'Current y (first ',I6,' components over a total of ', + I6,'):', + /,1(6X,6(1X,1PD11.4))) 9240 format( 6X,'Functional value = ',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) 9330 format( 6X,'The one-dimensional minimizer does not exist') 9340 format( 6X,'One-dimensional minimizer: ',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') 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,'INNIT = ',1X,I16) end c ****************************************************** c ****************************************************** subroutine evallnlispg(m,n,y,g,arow,acol,aval,anze,bbar,lambda,d, + amax,redc,L,Q,nL,w1,w2,w3,w4,theta,infty) C SCALAR ARGUMENTS integer m,n,anze double precision lambda,amax,redc,L,Q,theta,infty C ARRAY ARGUMENTS integer arow(anze),acol(anze) double precision y(m),g(n),aval(anze),d(n),bbar(m),nL(m),w1(n), + w2(m),w3(n),w4(n) C LOCAL SCALARS integer i,j double precision tmp C --------------------------------------------------------------- C Computing corresponding primal point d = - lambda ( 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 d = - lambda w1 do i = 1,n d(i) = lambda * w1(i) end do 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 if ( amax .lt. 0.0d0 ) then amax = infty do i = 1,m if ( w2(i) .gt. 0.0d0 ) then amax = min( amax, bbar(i) / w2(i) ) write(*,*) amax,bbar(i),w2(i) end if end do stop end if 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 --------------------------------------------------------------- C Computing 0.5 d^t Hk d = - d^t (gk + A^t y) = d^t w1 tmp = 0.0d0 do i = 1,n tmp = tmp + d(i) * w1(i) end do L = 0.5d0 * tmp Q = 0.5d0 * redc ** 2 * tmp C Computing gk^t d tmp = 0.0d0 do i = 1,n tmp = tmp + g(i) * d(i) end do L = L + tmp Q = Q + redc * tmp C Computing y^t ( A xk - b + A d ) = y^t ( w2 - bbar ) tmp = 0.0d0 do i = 1,m tmp = tmp + y(i) * ( w2(i) - bbar(i) ) end do L = L + tmp C --------------------------------------------------------------- C Computing - \nabla L(d,y) = A lambda (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 c **************************************************** c **************************************************** c Last update: March 26th, 2003 c See report of last modifications at the end of this file. subroutine gencan(n,x,l,u,epsgpen,epsgpsn,epsnfp,maxitnfp, *maxitngp,fmin,maxit,maxfc,udelta0,ucgmaxit,cgscre,cggpnf,cgepsi, *cgepsf,epsnqmp,maxitnqmp,nearlyq,gtype,htvtype,trtype,iprint, *f,g,gpeucn2,gpsupn,iter,fcnt,gcnt,cgcnt,spgiter,spgfcnt, *tniter,tnfcnt,tnstpcnt,tnintcnt,tnexgcnt,tnexbcnt,tnintfe,tnexgfe, *tnexbfe,flag,s,y,d,ind,w,eta,delmin,lammax,lammin,theta,gamma, *beta,sigma1,sigma2,nint,next,mininterp,ncomp,sterel,steabs,epsrel, *epsabs,infty) logical nearlyq integer n,maxitnfp,maxitngp,maxit,maxfc,maxitnqmp,ucgmaxit, *cgscre,gtype,htvtype,trtype,iprint,iter,fcnt,gcnt,cgcnt, *spgiter,spgfcnt,tniter,tnfcnt,tnstpcnt,tnintcnt,tnexgcnt,tnexbcnt, *tnintfe,tnexgfe,tnexbfe,flag,ind(n),mininterp,ncomp double precision x(n),l(n),u(n),epsgpen,epsgpsn,epsnfp,fmin, *udelta0,cggpnf,cgepsi,cgepsf,epsnqmp,f,g(n),gpeucn2,gpsupn,s(n), *y(n),d(n),w(5*n),eta,delmin,lammax,lammin,theta,gamma,beta,sigma1, *sigma2,nint,next,sterel,steabs,epsrel,epsabs,infty c Solves the box-constrained minimization problem c c Minimize f(x) c subject to l \leq x \leq u c c using a method described in c c E. G. Birgin and J. M. Martinez, "Large-scale active-set c box-constrained optimization method with spectral projected c gradients", Computational Optimization and Applications, c 23, pp. 101-125, 2002. c c Subroutines evalf and evalg must be supplied by the user to c evaluate the function f and its gradient, respectively. The c calling sequences are c c call evalf(n, x, f, inform) c c call evalg(n, x, g, inform) c c where x is the point where the function (the gradient) must be c evaluated, n is the number of variables and f (g) is the c functional value (the gradient vector). The real parameters c x, f, g must be double precision. c c A subroutine evalhd to compute the Hessian times vector products c is optional. If this subroutine is not provided an incremental c quotients version will be used instead. The calling sequence of c this subroutine should be c c call evalhd(nind,ind,n,x,u,hu,inform) c c where x is the point where the approx-Hessian is being considered, c n is the number of variables, u is the vector which should be c multiplied by the approx-Hessian H and hu is the vector where the c product should be placed. The information about the matrix H must c be passed to evalhd by means of common declarations. The necessary c computations must be done in evalg. The real parameters x, u, hu c must be double precision. c c This subroutine must be coded by the user, taking into account c that n is the number of variables of the problem and that hu must c be the product H u. Moreover, you must assume, when you code evalhd, 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 evalhd in such a way that hu is the vector whose i-th entry is c c hu(i) = \Sum_{j=1}^{nind} H_{i,ind(j)} u_ind(j) c c Moreover, the only components of hu that you need to compute are c those which corresponds to the indices ind(1),...,ind(nind). c However, observe that you must assume that, in u, the whole c vector is present, with its n components, even the zeroes. So, c if you decide to code evalhd without taking into account the c presence of ind and nind, you can do it. A final observation: c probably, if nind is close to n, it is not worthwhile to use ind, c due to the cost of accessing the correct indices. If you want, c you can test, within your evalhd, if (say) nind > n/2, and, in c this case use a straightforward scalar product for the components c of hu. c c Example: Suppose that the matrix H is full. The main steps of c evalhd could be: c c do i= 1, nind c indi= ind(i) c hu(indi)= 0.0d0 c do j= 1, nind c indj= ind(j) c hu(indi)= hu(indi) + H(indi,indj) * u(indj) c end do c end do c c c On Entry c c n integer c the space dimension c c x double precision x(n) c initial estimate to the solution c c l double precision l(n) c lower bounds c c u double precision u(n) c upper bounds c c epsgpen double precision c small positive number for declaring convergence when the c euclidian norm of the projected gradient is less than c or equal to epsgpen c c RECOMMENDED: epsgpen = 1.0d-5 c c epsgpsn double precision c small positive number for declaring convergence when the c infinite norm of the projected gradient is less than c or equal to epsgpsn c c RECOMMENDED: epsgpsn = 1.0d-5 c c epsnfp double precision c ``lack of enough progress'' measure. The algorithm stops by c ``lack of enough progress'' when f(x_k) - f(x_{k+1}) <= c epsnfp * max { f(x_j)-f(x_{j+1}, j= 0, gencan iterations information is printed. c If iprint >= 1, line searches and conjugate gradient information c is printed. c c RECOMMENDED: iprint = 1 c c s,y,d,w double precision s(n),y(n),d(n),w(5*n) c working vectors c c ind integer ind(n) c working vector c c eta double precision c constant for deciding abandon the current face or not c We abandon the current face if the norm of the internal c gradient (here, internal components of the continuous c projected gradient) is smaller than (1-eta) times the c norm of the continuous projected gradient. Using eta=0.9 c is a rather conservative strategy in the sense that c internal iterations are preferred over SPG iterations. c c RECOMMENDED: eta = 0.9 c c delmin double precision c minimum ``trust region'' to compute the Truncated Newton c direction c c RECOMMENDED: delmin = 0.1 c c lammin, lammax double precision c The spectral steplength, called lambda, is projected c inside the box [lammin,lammax] c c RECOMMENDED: lammin = 10^{-10} and lammax = 10^{10} c c theta double precision c constant for the angle condition, i.e., at iteration k c we need a direction d_k such that c <= -theta ||g||_2 ||d_k||_2, c where g_k is \nabla f(x_k) c c RECOMMENDED: theta = 10^{-6} c c gamma double precision c constant for the Armijo crtierion c f(x + alpha d) <= f(x) + gamma * alpha * <\nabla f(x),d> c c RECOMMENDED: gamma = 10^{-4} c c beta double precision c constant for the beta condition c .ge. beta * c if (x_k + d_k) satisfies the Armijo condition but does not c satisfy the beta condition then the point is accepted, but c if it satisfied the Armijo condition and also satisfies the c beta condition then we know that there is the possibility c for a succesful extrapolation c c RECOMMENDED: beta = 0.5 c c sigma1, sigma2 double precision c constant for the safeguarded interpolation c if alpha_new \notin [sigma1, sigma*alpha] then we take c alpha_new = alpha / nint c c RECOMMENDED: sigma1 = 0.1 and sigma2 = 0.9 c c nint double precision c constant for the interpolation. See the description of c sigma1 and sigma2 above. Sometimes we take as a new trial c step the previous one divided by nint c c RECOMMENDED: nint = 2.0 c c next double precision c constant for the extrapolation c when extrapolating we try alpha_new = alpha * next c c RECOMMENDED: next = 2.0 c c mininterp integer c constant for testing if, after having made at least mininterp c interpolations, the steplength is too small. In that case c failure of the line search is declared (may be the direction c is not a descent direction due to an error in the gradient c calculations) c c RECOMMENDED: mininterp = 4 c c (use mininterp > maxfc for inhibit this stopping criterion) c c ncomp integer c this constant is just for printing. In a detailed printing c option, ncomp component of the current point will be printed c c RECOMMENDED: ncomp = 5 c c sterel, steabs double precision c this constants mean a ``relative small number'' and ``an c absolute small number'' for the increments in finite c difference approximations of derivatives c c RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} c c epsrel, epsabs, infty double precision c this constants mean a ``relative small number'', ``an c absolute small number'', and ``infinite or a very big c number''. Basically, a quantity A is considered negligeble c with respect to another quantity B if c |A| < max ( epsrel * |B|, epsabs ) c c RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20} and c infty = 10^{+20} c c c On Return c c x double precision x(n) c final estimation to the solution c c f double precision c function value at the final estimation c c g double precision g(n) c gradient at the final estimation c c gpeucn2 double precision c squared 2-norm of the continuous projected c gradient g_p at the final estimation (||g_p||_2^2) c c gpsupn double precision c ||g_p||_inf at the final estimation c c iter integer c number of iterations c c fcnt integer c number of function evaluations c c gcnt integer c number of gradient evaluations c c cgcnt integer c number of conjugate gradient iterations c c spgiter integer c number of SPG iterations c c spgfcnt integer c number of function evaluations in SPG-directions line searches c c tniter integer c number of Truncated Newton iterations c c tnfcnt integer c number of function evaluations in TN-directions line searches c c tnintcnt integer c number of times a backtracking in a TN-drection was needed c c tnexgcnt integer c number of times an extrapolation in a TN-direction was c succesful in decreass the function value c c tnexbcnt integer c number of times an extrapolation was aborted in the first c extrapolated point by augment of the function value c c flag integer c This output parameter tells what happened in this c subroutine, according to the following conventions: c c 0= convergence with small euclidian-norm of the c projected gradient (smaller than epsgpen); c c 1= convergence with small infinite-norm of the c projected gradient (smaller than epsgpsn); c c 2= the algorithm stopped by ``lack of enough progress'', c that means that f(x_k) - f(x_{k+1}) <= c epsnfp * max { f(x_j)-f(x_{j+1}, j maxfc for inhibit this criterion; c c 7= it was achieved the maximum allowed number of c iterations (maxit); c c 8= it was achieved the maximum allowed number of c function evaluations (maxfc); c c 9= error in evalf subroutine; c c 10= error in evalg subroutine; c c 11= error in evalhd subroutine. character*3 ittype integer i,itnfp,itngp,cgflag,lsflag,nind,nprint,cgmaxit,cgiter, *fcntprev,tnintprev,tnexgprev,tnexbprev,rbdtype,rbdind,inform double precision lambda,gieucn2,gpi,fprev,sts,sty,epsgpen2,delta, *currprog,bestprog,gexp,gexpprev,cgeps,xnorm,ometa2,amax,amaxx, *acgeps,bcgeps,kappa,gpeucn20,gpsupn0 c ======================================================= c Initialization c ======================================================= c Set some initial values: c counters, iter= 0 fcnt= 0 gcnt= 0 cgcnt= 0 spgiter= 0 spgfcnt= 0 tniter= 0 tnfcnt= 0 tnstpcnt= 0 tnintcnt= 0 tnexgcnt= 0 tnexbcnt= 0 tnintfe= 0 tnexgfe= 0 tnexbfe= 0 c just for printing, nprint= min0(n,ncomp) c for testing convergence, epsgpen2= epsgpen**2 c for testing wether abandon the current face or not, c (ometa2 means '(one minus eta) squared') ometa2= (1.0d0-eta)**2 c for testing progress in f, fprev= infty bestprog= 0.0d0 itnfp= 0 c for testing progress in the order of the gradient norm, and gexpprev= infty itngp= 0 c Print problem information if(iprint.ge.0) then write(*,977) n write(*,978) nprint,(l(i),i=1,nprint) write(*,979) nprint,(u(i),i=1,nprint) write(*,980) nprint,(x(i),i=1,nprint) write(10,977) n write(10,978) nprint,(l(i),i=1,nprint) write(10,979) nprint,(u(i),i=1,nprint) write(10,980) nprint,(x(i),i=1,nprint) end if c Project initial guess. If the initial guess is infeasible, c projection puts it in the box. do i= 1, n x(i)= dmax1(l(i),dmin1(x(i),u(i))) end do c Compute x euclidian-norm xnorm= 0.0d0 do i= 1, n xnorm= xnorm + x(i)**2 end do xnorm= dsqrt(xnorm) c Compute function and gradient at the initial point call evalq(n,x,f,inform) fcnt= fcnt + 1 if (inform.ne.0) then flag = 9 if (iprint.ge.0) then write(*,999) flag write(10,999) flag end if go to 500 end if if (gtype.eq.0) then call evalnq(n,x,g,inform) else if (gtype.eq.1) then call evalgdiff(n,x,g,sterel,steabs,inform) end if gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.0) then write(*,1000) flag write(10,1000) flag end if go to 500 end if c Compute continuous project gradient infinite- and euclidian- c norm, internal gradient euclidian-norm, and store in nind the c number of free variables and in array p their identifiers. nind= 0 gpsupn= 0.0d0 gpeucn2= 0.0d0 gieucn2= 0.0d0 do i= 1, n gpi= dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i) gpsupn= dmax1(gpsupn, dabs(gpi)) gpeucn2= gpeucn2 + gpi**2 if (x(i).gt.l(i).and.x(i).lt.u(i)) then gieucn2= gieucn2 + gpi**2 nind= nind + 1 ind(nind)= i end if end do c Compute a linear relation between gpeucn2 and cgeps2, i.e., c scalars a and b such that c c a * log10(||g_P(x_0)||_2^2) + b = log10(cgeps_0^2) and c c a * log10(||g_P(x_f)||_2^2) + b = log10(cgeps_f^2), c c where cgeps_0 and cgeps_f are provided. Note that if c cgeps_0 is equal to cgeps_f then cgeps will be allways c equal to cgeps_0 and cgeps_f. c We introduce now a linear relation between gpsupn and cgeps also. if (cgscre.eq.1) then acgeps= 2.0d0*dlog10(cgepsf/cgepsi)/dlog10(cggpnf**2/gpeucn2) bcgeps= 2.0d0*dlog10(cgepsi)-acgeps*dlog10(gpeucn2) else ! if (cgscre.eq.2) then acgeps= dlog10(cgepsf/cgepsi)/dlog10(cggpnf/gpsupn) bcgeps= dlog10(cgepsi)-acgeps*dlog10(gpsupn) end if c And it will be used for the linear relation of cgmaxit gpsupn0= gpsupn gpeucn20= gpeucn2 c Print initial information if(iprint.ge.0) then write(*,981) iter write(*,985) nprint,(x(i),i=1,nprint) write(*,986) nprint,(g(i),i=1,nprint) write(*,987) nprint, * (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint) write(*,988) min0(nprint,nind),nind, * (ind(i),i=1,min0(nprint,nind)) write(*,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n, * spgiter,tniter,fcnt,gcnt,cgcnt write(10,981) iter write(10,985) nprint,(x(i),i=1,nprint) write(10,986) nprint,(g(i),i=1,nprint) write(10,987) nprint, * (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint) write(10,988) min0(nprint,nind),nind, * (ind(i),i=1,min0(nprint,nind)) write(10,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n, * spgiter,tniter,fcnt,gcnt,cgcnt end if c Test whether the initial functional value is very small if (f.le.fmin) then flag= 4 if (iprint.ge.0) then write(*,994) flag,fmin write(10,994) flag,fmin end if go to 500 end if c Test whether the number of functional evaluations is c exhausted (i.e., maxfc = 1) if (fcnt.ge.maxfc) then flag= 8 if (iprint.ge.0) then write(*,998) flag,maxfc write(10,998) flag,maxfc end if go to 500 end if c ======================================================= c Main loop c ======================================================= 100 continue c ======================================================= c Test stopping criteria c ======================================================= c Test whether the continuous projected gradient euclidian-norm c is small enough to declare convergence if (gpeucn2.le.epsgpen2) then flag= 0 if (iprint.ge.0) then write(*,990) flag,epsgpen write(10,990) flag,epsgpen end if go to 500 end if c Test whether the continuous projected gradient infinite-norm c is small enough to declare convergence if (gpsupn.le.epsgpsn) then flag= 1 if (iprint.ge.0) then write(*,991) flag,epsgpsn write(10,991) flag,epsgpsn end if go to 500 end if c Test whether we performed many iterations without good progress c of the functional value currprog= fprev - f bestprog= dmax1(currprog, bestprog) if (currprog.le.epsnfp*bestprog) then itnfp= itnfp + 1 if (itnfp.ge.maxitnfp) then flag= 2 if (iprint.ge.0) then write(*,992) flag,epsnfp,maxitnfp write(10,992) flag,epsnfp,maxitnfp end if go to 500 endif else itnfp= 0 endif c Test whether we have performed many iterations without good c reduction of the euclidian-norm of the projected gradient gexp= dlog10(gpeucn2) if (gexp.ge.gexpprev) then itngp= itngp + 1 if(itngp.ge.maxitngp) then flag= 3 if (iprint.ge.0) then write(*,993) flag,maxitngp write(10,993) flag,maxitngp end if go to 500 endif else itngp= 0 endif gexpprev= gexp c Test whether the number of iterations is exhausted if (iter.ge.maxit) then flag= 7 if (iprint.ge.0) then write(*,997) flag,maxit write(10,997) flag,maxit end if go to 500 end if c Note that the stopping criteria related to small functional c values ( f <= fmin ) and number of functional evaluations c are detected inside the line searches and do not need to c be tested here c ======================================================= c The stopping criteria were not satisfied, a new c iteration will be made c ======================================================= iter= iter + 1 c ======================================================= c Save current values, f, x and g c ======================================================= fprev= f do i= 1, n s(i)= x(i) y(i)= g(i) end do c ======================================================= c Compute new iterate c ======================================================= c We abandon the current face if the norm of the internal gradient c (here, internal components of the continuous projected gradient) c is smaller than (1-eta) times the norm of the continuous projected c gradient. Using eta=0.9 is a rather conservative strategy in the c sense that internal iterations are preferred over SPG iterations. c Replace eta=0.9 by other tolerance in (0,1) if you find it convenient. if (gieucn2.le.ometa2*gpeucn2) then c =================================================== c Some constraints should be abandoned. Compute c the new iterate using an SPG iteration c =================================================== ittype= 'SPG' spgiter= spgiter + 1 c Compute spectral steplength if (iter.eq.1.or.sty.le.0.0d0) then lambda= dmax1(1.0d0,xnorm)/dsqrt(gpeucn2) else lambda= sts/sty end if lambda= dmin1(lammax,dmax1(lammin,lambda)) c Perform a line search with safeguarded quadratic c interpolation along the direction of the spectral c continuous projected gradient fcntprev= fcnt call spgls(n,x,f,g,l,u,lambda,fmin,maxfc,iprint,fcnt,lsflag, * w(1),w(n+1),gamma,sigma1,sigma2,nint,mininterp,sterel,steabs, * epsrel,epsabs,infty) spgfcnt= spgfcnt + (fcnt-fcntprev) if (lsflag.eq.9) then flag = 9 if (iprint.ge.0) then write(*,999) flag write(10,999) flag end if go to 500 end if c Compute the gradient of the new iterate if (gtype.eq.0) then call evalnq(n,x,g,inform) else if (gtype.eq.1) then call evalgdiff(n,x,g,sterel,steabs,inform) end if gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.0) then write(*,1000) flag write(10,1000) flag end if go to 500 end if else c =================================================== c The new iterate will belong to the closure of c the current face c =================================================== ittype= 'TN ' tniter= tniter + 1 c Compute trust-region radius if (iter.eq.1) then if(udelta0.lt.0.d0) then delta= dmax1(delmin,100.0d0*dmax1(1.0d0,xnorm)) else delta= udelta0 end if else delta= dmax1(delmin,10.0d0*dsqrt(sts)) end if c Shrink the point, its gradient and the bounds call Gshrink(nind,ind,n,x,g,l,u) c Compute the descent direction solving the newtonian c system by conjugate gradients c Set conjugate gradient stopping criteria. Default values are c taken if you set ucgeps < 0 and ucgmaxit < 0, respectively. c Otherwise, the parameters cgeps and cgmaxit will be the ones c set by the user. if(ucgmaxit.lt.0) then if (nearlyq) then cgmaxit= nind else if (cgscre.eq.1) then kappa= dlog10(gpeucn2/gpeucn20)/ * dlog10(epsgpen2/gpeucn20) else ! if (cgscre.eq.2) then kappa= dlog10(gpsupn/gpsupn0)/ * dlog10(epsgpsn/gpsupn0) end if kappa= dmax1(0.0d0,dmin1(1.0d0,kappa)) cgmaxit= int( * (1.0d0-kappa)*dmax1(1.0d0,10.0d0*dlog10(dfloat(nind))) * + kappa*dfloat(nind) ) end if c cgmaxit = 2 * nind else cgmaxit= ucgmaxit end if if (cgscre.eq.1) then cgeps= dsqrt(10.0d0**(acgeps*dlog10(gpeucn2)+bcgeps)) else ! if (cgscre.eq.2) then cgeps= 10.0d0**(acgeps*dlog10(gpsupn)+bcgeps) end if cgeps= dmax1(cgepsf,dmin1(cgepsi,cgeps)) c Call conjugate gradients call cg(nind,ind,n,x,g,delta,l,u,cgeps,epsnqmp,maxitnqmp, * cgmaxit,nearlyq,gtype,htvtype,trtype,iprint,ncomp,d,cgiter, * rbdtype,rbdind,cgflag,w(1),w(n+1),w(2*n+1),w(3*n+1),w(4*n+1), * theta,sterel,steabs,epsrel,epsabs,infty) cgcnt= cgcnt + cgiter if (cgflag.eq.10) then flag = 10 if (iprint.ge.0) then write(*,1000) flag write(10,1000) flag end if go to 500 else if (cgflag.eq.11) then flag = 11 if (iprint.ge.0) then write(*,1001) flag write(10,1001) flag end if go to 500 end if c Compute maximum step if (cgflag.eq.2) then amax= 1.0d0 else amax= infty do i= 1, nind if (d(i).gt.0.0d0.and.u(i).lt.infty) then amaxx= (u(i)-x(i))/d(i) if (amaxx.lt.amax) then amax= amaxx rbdind= i rbdtype= 2 end if else if (d(i).lt.0.0d0.and.l(i).gt.-infty) then amaxx= (l(i)-x(i))/d(i) if (amaxx.lt.amax) then amax= amaxx rbdind= i rbdtype= 1 end if end if end do end if c Perform the line search tnintprev= tnintcnt tnexgprev= tnexgcnt tnexbprev= tnexbcnt fcntprev= fcnt call tnls(nind,ind,n,x,f,g,d,l,u,amax,rbdtype,rbdind,fmin, * maxfc,gtype,iprint,fcnt,gcnt,tnintcnt,tnexgcnt,tnexbcnt, * lsflag,w(1),w(n+1),gamma,beta,sigma1,sigma2,nint,next, * mininterp,sterel,steabs,epsrel,epsabs,infty) if (lsflag.eq.9) then flag = 9 if (iprint.ge.0) then write(*,999) flag write(10,999) flag end if go to 500 else if (lsflag.eq.10) then flag = 10 if (iprint.ge.0) then write(*,1000) flag write(10,1000) flag end if go to 500 end if tnfcnt= tnfcnt + (fcnt-fcntprev) if (tnintcnt.gt.tnintprev) then tnintfe= tnintfe + (fcnt-fcntprev) else if (tnexgcnt.gt.tnexgprev) then tnexgfe= tnexgfe + (fcnt-fcntprev) else if (tnexbcnt.gt.tnexbprev) then tnexbfe= tnexbfe + (fcnt-fcntprev) else tnstpcnt= tnstpcnt + 1 end if c Expand the point, its gradient and the bounds call Gexpand(nind,ind,n,x,g,l,u) c If the line search (interpolation) in the Truncated Newton c direction stopped due to a very small step (lsflag=6) we c will discard this iteration and force a SPG iteration c Note that tnls subroutine was coded in such a way that in c case of lsflag=6 termination the subroutine discards all c what was made and returns with the same point it started if (lsflag.eq.6) then if (iprint.ge.0) then write(*,*) write(*,*) * ' The previous NT iteration was discarded due to', * ' a termination for very small step in the line ', * ' search. A SPG iteration will be forced now. ' write(10,*) write(10,*) * ' The previous NT iteration was discarded due to', * ' a termination for very small step in the line ', * ' search. A SPG iteration will be forced now. ' end if ittype= 'SPG' spgiter= spgiter + 1 c Compute spectral steplength if (iter.eq.1.or.sty.le.0.0d0) then lambda= dmax1(1.0d0,xnorm)/dsqrt(gpeucn2) else lambda= sts/sty end if lambda= dmin1(lammax,dmax1(lammin,lambda)) c Perform a line search with safeguarded quadratic c interpolation along the direction of the spectral c continuous projected gradient fcntprev= fcnt call spgls(n,x,f,g,l,u,lambda,fmin,maxfc,iprint,fcnt, * lsflag,w(1),w(n+1),gamma,sigma1,sigma2,nint,mininterp, * sterel,steabs,epsrel,epsabs,infty) spgfcnt= spgfcnt + (fcnt-fcntprev) if (lsflag.eq.9) then flag = 9 if (iprint.ge.0) then write(*,999) flag write(10,999) flag end if go to 500 end if c Compute the gradient in the new iterate if (gtype.eq.0) then call evalnq(n,x,g,inform) else if (gtype.eq.1) then call evalgdiff(n,x,g,sterel,steabs,inform) end if gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.0) then write(*,1000) flag write(10,1000) flag end if go to 500 end if end if end if c ======================================================= c Prepare for the next iteration c ======================================================= c This adjustment/projection is "por lo que las putas pudiera" do i= 1, n if (x(i).le.l(i)+dmax1(epsrel*dabs(l(i)),epsabs)) then x(i)= l(i) else if (x(i).ge.u(i)-dmax1(epsrel*dabs(u(i)),epsabs)) then x(i)= u(i) end if end do c Compute x euclidian-norm xnorm= 0.0d0 do i= 1, n xnorm= xnorm + x(i)**2 end do xnorm= dsqrt(xnorm) c Compute s = x_{k+1} - x_k, y = g_{k+1} - g_k, and sts= 0.0d0 sty= 0.0d0 do i= 1, n s(i)= x(i)-s(i) y(i)= g(i)-y(i) sts= sts + s(i)*s(i) sty= sty + s(i)*y(i) end do c Compute continuous project gradient infinite- and euclidian- c norm, internal gradient euclidian-norm, and store in nind the c number of free variables and in array p its identifiers. nind= 0 gpsupn= 0.0d0 gpeucn2= 0.0d0 gieucn2= 0.0d0 do i= 1, n gpi= dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i) gpsupn= dmax1(gpsupn, dabs(gpi)) gpeucn2= gpeucn2 + gpi**2 if (x(i).gt.l(i).and.x(i).lt.u(i)) then gieucn2= gieucn2 + gpi**2 nind= nind + 1 ind(nind)= i end if end do c Print information of this iteration if (iprint.ge.0) then write(*,983) iter,ittype write(*,985) nprint,(x(i),i=1,nprint) write(*,986) nprint,(g(i),i=1,nprint) write(*,987) nprint, * (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint) write(*,988) min0(nprint,nind),nind, * (ind(i),i=1,min0(nprint,nind)) write(*,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n, * spgiter,tniter,fcnt,gcnt,cgcnt write(10,983) iter,ittype write(10,985) nprint,(x(i),i=1,nprint) write(10,986) nprint,(g(i),i=1,nprint) write(10,987) nprint, * (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint) write(10,988) min0(nprint,nind),nind, * (ind(i),i=1,min0(nprint,nind)) write(10,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n, * spgiter,tniter,fcnt,gcnt,cgcnt end if c ======================================================= c Test some stopping criteria that may occur inside the c line searches c ======================================================= if (lsflag.eq.4) then flag= 4 if (iprint.ge.0) then write(*,994) flag,fmin write(10,994) flag,fmin end if go to 500 else if (lsflag.eq.6) then flag= 6 if (iprint.ge.0) then write(*,996) flag,mininterp,epsrel,epsabs write(10,996) flag,mininterp,epsrel,epsabs end if go to 500 else if (lsflag.eq.8) then flag= 8 if (iprint.ge.0) then write(*,998) flag,maxfc write(10,998) flag,maxfc end if go to 500 end if c ======================================================= c Iterate c ======================================================= go to 100 c ======================================================= c End of main loop c ======================================================= c ======================================================= c Report output status and return c ======================================================= 500 continue c Print final information if (iprint.ge.0) then write(*,982) iter write(*,985) nprint,(x(i),i=1,nprint) write(*,986) nprint,(g(i),i=1,nprint) write(*,987) nprint, * (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint) write(*,988) min0(nprint,nind),nind, * (ind(i),i=1,min0(nprint,nind)) write(*,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n, * spgiter,tniter,fcnt,gcnt,cgcnt write(10,982) iter write(10,985) nprint,(x(i),i=1,nprint) write(10,986) nprint,(g(i),i=1,nprint) write(10,987) nprint, * (dmin1(u(i),dmax1(l(i),x(i)-g(i)))-x(i),i=1,nprint) write(10,988) min0(nprint,nind),nind, * (ind(i),i=1,min0(nprint,nind)) write(10,1002) f,dsqrt(gpeucn2),dsqrt(gieucn2),gpsupn,nind,n, * spgiter,tniter,fcnt,gcnt,cgcnt end if return c Non-executable statements 977 format(/1X, 'Entry to GENCAN. Number of variables: ',I7) 978 format(/1X,'Lower bounds (first ',I6, ' components): ', */,6(1X,1PD11.4)) 979 format(/1X,'Upper bounds (first ',I6, ' components): ', */,6(1X,1PD11.4)) 980 format(/1X,'Initial point (first ',I6, ' components): ', */,6(1X,1PD11.4)) 981 format(/1X,'GENCAN iteration: ',I6, ' (Initial point)') 982 format(/1X,'GENCAN iteration: ',I6, ' (Final point)') 983 format(/,1X,'GENCAN iteration: ',I6, *' (This point was obtained using a ',A3,' iteration)') 985 format(1X,'Current point (first ',I6, ' components): ', */,6(1X,1PD11.4)) 986 format(1X,'Current gradient (first ',I6, ' components): ', */,6(1X,1PD11.4)) 987 format(1X,'Current continuous projected gradient (first ',I6, *' components): ',/,6(1X,1PD11.4)) 988 format(1X,'Current free variables (first ',I6, *', total number ',I6,'): ',/,10(1X,I6)) 990 format(/1X,'Flag of GENCAN= ',I2, *' (convergence with Euclidian-norm of the projected gradient', */,1X,'smaller than ',1PD11.4,')') 991 format(/1X,'Flag of GENCAN= ',I2, *' (convergence with sup-norm of the projected gradient', */,1X,'smaller than ',1PD11.4,')') 992 format(/1X,'Flag of GENCAN= ',I2, *' (The algorithm stopped by lack of enough progress. This means', */,1X,'that f(x_k) - f(x_{k+1}) .le. ',1PD11.4, *' * max [ f(x_j)-f(x_{j+1}, j < k ]',/,1X,'during ',I7, *' consecutive iterations') 993 format(/1X,'Flag of GENCAN= ',I2, *' (The algorithm stopped because the order of the', */,1X,'Euclidian-norm of the continuous projected gradient did', *' not change during ',/,1X,I7,' consecutive iterations.', *' Probably, an exagerately small norm of the',/,1X,'continuous', *' projected gradient is required for declaring convergence') 994 format(/1X,'Flag of GENCAN= ',I2, *' (The algorithm stopped because the functional value is', */,1X,'smaller than ',1PD11.4) 996 format(/1X,'Flag of GENCAN= ',I2, *' (Too small step in a line search. After having made at least ', */,1X,I7,' interpolations, the steplength becames small. Small', *' means that we were',/,1X,'at point x with direction d and took', *' a step alpha such that',/,1X,'alpha * |d_i| .lt. max [', *1PD11.4,' * |x_i|,',1PD11.4,' ] for all i') 997 format(/1X,'Flag of GENCAN= ',I2, *' (It was exceeded the maximum allowed number of iterations', */,1X,'(maxit=',I7,')') 998 format(/1X,'Flag of GENCAN= ',I2, *' (It was exceeded the maximum allowed number of functional', */,1X,'evaluations (maxfc=',I7,')') 999 format(/1X,'Flag of GENCAN= ',I2, *' (Error in evalf subroutine)') 1000 format(/1X,'Flag of GENCAN= ',I2, *' (Error in evalg subroutine)') 1001 format(/1X,'Flag of GENCAN= ',I2, *' (Error in evalhd subroutine)') 1002 format(1X,'Functional value: ', 1PD11.4, */,1X,'Euclidian-norm of the continuous projected gradient: ', *1PD11.4, */,1X,'Euclidian-norm of the internal projection of gp: ',1PD11.4, */,1X,'Sup-norm of the continuous projected gradient: ',1PD11.4, */,1X,'Free variables at this point: ',I7, *' (over a total of ',I7,')', */,1X,'SPG iterations: ',I7, */,1X,'TN iterations: ',I7, */,1X,'Functional evaluations: ',I7, */,1X,'Gradient evaluations: ',I7, */,1X,'Conjugate gradient iterations: ',I7) end c ****************************************************** c ****************************************************** subroutine spgls(n,x,f,g,l,u,lambda,fmin,maxfc,iprint,fcnt,flag, *xtrial,d,gamma,sigma1,sigma2,nint,mininterp,sterel,steabs,epsrel, *epsabs,infty) integer n,maxfc,iprint,fcnt,flag,mininterp double precision x(n),f,g(n),l(n),u(n),lambda,fmin,xtrial(n),d(n), *gamma,sigma1,sigma2,nint,sterel,steabs,epsrel,epsabs,infty c Safeguarded quadratic interpolation, used in SPG. c c On Entry c c n integer c the order of the x c c x double precision x(n) c current point c c f double precision c function value at the current point c c g double precision g(n) c gradient vector at the current point c c l double precision l(n) c lower bounds c c u double precision u(n) c upper bounds c c lambda double precision c spectral steplength c c fmin double precision c functional value for the stopping criteria f <= fmin c c maxfc integer c maximum number of funtion evaluations c c output logical c TRUE: print some information at each iteration, c FALSE: no print. c c xtrial, d double precision xtrial(n), d(n) c working vectors c c gamma double precision c constant for the Armijo crtierion c f(x + alpha d) <= f(x) + gamma * alpha * <\nabla f(x),d> c c RECOMMENDED: gamma = 10^{-4} c c sigma1, sigma2 double precision c constant for the safeguarded interpolation c if alpha_new \notin [sigma1, sigma*alpha] then we take c alpha_new = alpha / nint c c RECOMMENDED: sigma1 = 0.1 and sigma2 = 0.9 c c nint double precision c constant for the interpolation. See the description of c sigma1 and sigma2 above. Sometimes we take as a new trial c step the previous one divided by nint c c RECOMMENDED: nint = 2.0 c c mininterp integer c constant for testing if, after having made at least mininterp c interpolations, the steplength is soo small. In that case c failure of the line search is declared (may be the direction c is not a descent direction due to an error in the gradient c calculations) c c RECOMMENDED: mininterp = 4 c c sterel, steabs double precision c this constants mean a ``relative small number'' and ``an c absolute small number'' for the increments in finite c difference approximations of derivatives c c RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} c c epsrel, epsabs, infty double precision c this constants mean a ``relative small number'', ``an c absolute small number'', and ``infinite or a very big c number''. Basically, a quantity A is considered negligeble c with respect to another quantity B if c |A| < max ( epsrel * |B|, epsabs ) c c RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20} and c infty = 10^{+20} c c On Return c c x double precision c final estimation of the solution c c f double precision c functional value at the final estimation c c fcnt integer c number of function evaluations used in the line search c c flag integer c This output parameter tells what happened in this c subroutine, according to the following conventions: c c 0= convergence with an Armijo-like criterion c (f(xnew) <= f(x) + 1.0d-4 * alpha * ); c c 4= the algorithm stopped because the functional value c is very small (f <= fmin); c c 6= too small step in the line search. After having made at c least mininterp interpolations, the steplength becames c small. ``small steplength'' means that we are at point c x with direction d and step alpha, and, for all i, c c |alpha * d(i)| .le. max ( epsrel * |x(i)|, epsabs ). c c In that case failure of the line search is declared c (maybe the direction is not a descent direction c due to an error in the gradient calculations). Use c mininterp > maxfc for inhibit this criterion; c c 8= it was achieved the maximum allowed number of c function evaluations (maxfc); c c 9= error in evalf subroutine. c We use the Armijo parameter gamma=1.0d-4. You can replace it by c some other number in (0,1) if you find it convenient. logical samep integer i,interp,inform double precision ftrial,gtd,alpha,atemp c Initialization interp= 0 c Compute first trial point, spectral projected gradient c direction, and directional derivative . alpha= 1.0d0 gtd= 0.0d0 do i= 1, n xtrial(i)= dmin1(u(i),dmax1(l(i),x(i)-lambda*g(i))) d(i)= xtrial(i)-x(i) gtd= gtd + g(i)*d(i) end do call evalq(n,xtrial,ftrial,inform) fcnt= fcnt + 1 if (inform.ne.0) then flag = 9 if (iprint.ge.1) then write(*,1000) flag write(10,1000) flag end if go to 500 end if c Print initial information if (iprint.ge.1) then write(*,980) lambda write(*,999) alpha,ftrial,fcnt write(10,980) lambda write(10,999) alpha, ftrial, fcnt end if c Main loop 100 continue c Test Armijo stopping criterion if (ftrial.le.f+gamma*alpha*gtd) then f= ftrial do i= 1, n x(i)= xtrial(i) end do flag= 0 if (iprint.ge.1) then write(*,990) flag write(10,990) flag end if go to 500 end if c Test whether f is very small if (ftrial.le.fmin) then f= ftrial do i= 1, n x(i)= xtrial(i) end do flag= 4 if (iprint.ge.1) then write(*,994) flag write(10,994) flag end if go to 500 end if c Test whether the number of functional evaluations is exhausted if (fcnt.ge.maxfc) then if (ftrial.lt.f) then f= ftrial do i= 1, n x(i)= xtrial(i) end do end if flag= 8 if (iprint.ge.1) then write(*,998) flag write(10,998) flag end if go to 500 end if c Compute new step (safeguarded quadratic interpolation) interp= interp + 1 if (alpha.lt.sigma1) then alpha= alpha/nint else atemp= (-gtd*alpha**2)/(2.0d0*(ftrial-f-alpha*gtd)) if (atemp.lt.sigma1 .or. atemp.gt.sigma2*alpha) then alpha= alpha/nint else alpha= atemp end if end if c Compute new trial point do i= 1, n xtrial(i)= x(i) + alpha*d(i) end do call evalq(n,xtrial,ftrial,inform) fcnt= fcnt + 1 if (inform.ne.0) then flag = 9 if (iprint.ge.1) then write(*,1000) flag write(10,1000) flag end if go to 500 end if c Print information of this iteration if (iprint.ge.1) then write(*,999) alpha,ftrial,fcnt write(10,999) alpha,ftrial,fcnt end if c Test whether at least mininterp interpolations were made and c two consecutive iterates are close enough samep= .true. do i= 1, n if (dabs(alpha*d(i)).gt.dmax1(epsrel*dabs(x(i)),epsabs)) then samep= .false. end if end do if (interp.gt.mininterp.and.samep) then if (ftrial.lt.f) then f= ftrial do i= 1, n x(i)= xtrial(i) end do end if flag= 6 if (iprint.ge.1) then write(*,996) flag write(10,996) flag end if go to 500 end if c Iterate go to 100 c Return 500 continue return c Non-executable statements 980 format(/,6x,'SPG (spectral steplength ',1PD11.4,')', */,/,6x,'SPG Line search') 999 format(6x,'Alpha= ',1PD11.4,' F= ',1PD11.4,' FE= ',I5) 990 format(6x,'Flag of SPG Line search= ',I2, *' (Convergence with an Armijo-like criterion)') 994 format(6x,'Flag of SPG Line search= ',I2, *' (Small functional value, smaller than ',/,6X,'parameter fmin)') 996 format(6x,'Flag of SPG Line search= ',I2, *' (Too small step in the interpolation process)') 998 format(6x,'Flag of SPG Line search= ',I2, *' (Too many functional evaluations)') 1000 format(6x,'Flag of SPG Line search= ',I2, *' (Error in evalf subroutine)') end c ****************************************************** c ****************************************************** subroutine cg(nind,ind,n,x,g,delta,l,u,eps,epsnqmp,maxitnqmp, *maxit,nearlyq,gtype,htvtype,trtype,iprint,ncomp,s,iter,rbdtype, *rbdind,flag,w,y,r,d,sprev,theta,sterel,steabs,epsrel,epsabs,infty) logical nearlyq integer nind,ind(nind),n,maxitnqmp,maxit,gtype,htvtype,trtype, *iprint,ncomp,iter,rbdtype,rbdind,flag double precision x(n),g(n),delta,l(n),u(n),eps,epsnqmp,s(n),w(n), *y(n),r(n),d(n),sprev(n),theta,sterel,steabs,epsrel,epsabs,infty c This subrotuine implements the conjugate-gradient method for c minimizing the quadratic approximation q(s) of f(x) at x, where c c q(s) = 1/2 s^T H s + g^T s, c c H = \nabla^2 f(x), c c g = \nabla f(x), c c subject to ||s||_2 <= delta and l <= x + s <= u. c c The method returns an approximation s to the solution such c that ||H s + g||_2 <= eps * ||g||_2; or converges to the c boundary of ||s||_2 <= delta and l <= x + s <= u; or finds c a point s and a direction d such that q(s + alpha d) = q(s) c for any alpha, i.e., d^T H d = g^T d = 0. c c On Entry c c nind integer c number of free variables (this is thee dimension in c which this subroutine will work) c c ind integer ind(n) c array which contains, in the first nind positions, the c identifiers of the free variables c c n integer c dimension of the full space c c x double precision x(n) c point at which f function is being approximated by the c quadratic model c c The first nind positions of x contains the free variables c x_ind(1), x_ind(2), ..., x_ind(nind). c c g double precision g(n) c linear coeficient of the quadratic function c c This is \nabla f(x) and it also contains in the first nind c positions the components g_ind(1), g_ind(2), ..., g_ind(nind). c c IMPORTANT: the linear algebra of this subroutine lies in c a space of dimension nind. The value of the full dimension n, c the non-free variables (which are at the end of array x) and c its gradient components (which are at the and of array g) c are, at this moment, being used to approximate the hessian c times vector products by incremental quotients. c c delta double precision c trust region radius (||s||_2 <= delta) c c l double precision l(n) c lower bounds on x + s. It components are ordered in the c same way as x and g. c c u double precision u(n) c upper bounds on x + s. It components are ordered in the c same way as x, g and l. c c eps double precision c tolerance for the stopping criterion c ||H s + g||_2 < eps * ||g||_2 c c epsnqmp double precision c see below c c maxitnqmp integer c This and the previous one parameter are used for a stopping c criterion of the conjugate gradient subalgorithm. If the c progress in the quadratic model is less or equal than a c fraction of the best progress ( epsnqmp * bestprog ) during c maxitnqmp consecutive iterations then CG is stopped by not c enough progress of the quadratic model. c c RECOMMENDED: epsnqmp = 1.0d-4, maxitnqmp = 5 c c maxit integer c maximum number of iterations allowed c c nearlyq logical c if function f is (nearly) quadratic, use the option c nerlyq = TRUE. Otherwise, keep the default option. c c if, in an iteration of CG we find a direction d such c that d^T H d <= 0 then we take the following decision: c c (i) if nearlyq = TRUE then take direction d and try to c go to the boundary chosing the best point among the two c point at the boundary and the current point. c c (ii) if nearlyq = FALSE then we stop at the current point. c c RECOMMENDED: nearlyq = FALSE c c gtype integer c type of gradient calculation c gtype = 0 means user suplied evalg subroutine, c gtype = 1 means central diference approximation. c c RECOMMENDED: gtype = 0 c c (provided you have the evalg subroutine) c c htvtype integer c type of gradient calculation c htvtype = 0 means user suplied evalhd subroutine, c htvtype = 1 means incremental quotients approximation. c c RECOMMENDED: htvtype = 1 c c (you take some risk using this option but, a menos que c voce tenha uma boa evalhd subroutine, incremental c quotients is a very cheap option) c c trtype integer c type of trust-region radius c trtype = 0 means 2-norm trust-region c trtype = 1 means infinite-norm trust-region c c RECOMMENDED: trtype = 0 c c output logical c TRUE: print some information at each iteration, c FALSE: no print. c c w,y,r,d,sprev double precision w(n),y(n),r(n),d(n),sprev(n) c working vectors c c theta double precision c constant for the angle condition, i.e., at iteration k c we need a direction d_k such that c <= -theta ||g||_2 ||d_k||_2, c where g_k is \nabla f(x_k) c c RECOMMENDED: theta = 10^{-6} c c sterel, steabs double precision c this constants mean a ``relative small number'' and ``an c absolute small number'' for the increments in finite c difference approximations of derivatives c c RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} c c epsrel, epsabs, infty double precision c this constants mean a ``relative small number'', ``an c absolute small number'', and ``infinite or a very big c number''. Basically, a quantity A is considered negligeble c with respect to another quantity B if c |A| < max ( epsrel * |B|, epsabs ) c c RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20} and c infty = 10^{+20} c c On Return c c s double precision s(n) c final estimation of the solution c c iter integer c number of conjugate-gradient iterations performed c c flag integer c termination parameter: c c 0= convergence with ||H s + g||_2 <= eps * ||g||_2; c c 1= convergence to the boundary of ||s||_2 <= delta; c c 2= convergence to the boundary of l - x <= s <= u - x; c c 3= stopping with s= s_k such that c <= -theta ||g||_2 ||s_k||_2 and c > -theta ||g||_2 ||s_{k+1}||_2; c c 4= not enough progress of the quadratic model during c maxitnqmp iterations, i.e., during maxitnqmp iterations c |q-qprev| .le. max ( epsrel * |q|, epsabs ); c c 6= very similar consecutive iterates, c for two consecutive iterates x and y, for all i c |x(i)-y(i)| .le. max ( epsrel * |x(i)|, epsabs ); c c 7= stopping with d such that d^T H d = 0 and g^T d = 0; c c 8= too many iterations; c c 11= error in evalhd subroutine. character*5 rbdtypea logical goth,samep integer i,itnqmp,rbdposatype,rbdnegatype,rbdposaind,rbdnegaind, *inform double precision rnorm2,rnorm2prev,beta,alpha,amax,amax1,amax2, *dtw,dnorm2,dts,aa,bb,cc,dd,gts,dtr,gnorm2,snorm2,q,snorm2prev, *qprev,amax1n,amax2n,amaxn,amax2x,amax2nx,qamax,qamaxn,norm2s, *bestprog,currprog c ======================================================= c Initialization c ======================================================= goth= .false. gnorm2= norm2s(nind,g) iter= 0 itnqmp= 0 qprev= infty bestprog= 0.0d0 do i= 1, nind s(i)= 0.0d0 r(i)= g(i) end do q= 0.0d0 gts= 0.0d0 snorm2= 0.0d0 rnorm2= gnorm2 c ======================================================= c Print initial information c ======================================================= if (iprint.ge.1) then write(*,980) maxit,eps if (trtype.eq.0) then write(*,981) delta else if (trtype.eq.1) then write(*,982) delta else write(*,983) end if write(*,984) iter,rnorm2,dsqrt(snorm2),q write(10,980) maxit,eps if (trtype.eq.0) then write(10,981) delta else if (trtype.eq.1) then write(10,982) delta else write(10,983) end if write(10,984) iter,rnorm2,dsqrt(snorm2),q end if c ======================================================= c Main loop c ======================================================= 100 continue c ======================================================= c Test stopping criteria c ======================================================= c if ||r||_2 = ||H s + g||_2 <= eps * ||g||_2 then stop if (rnorm2.le.eps**2*gnorm2.or. *(rnorm2.le.1.0d-10.and.iter.ne.0)) then flag= 0 if (iprint.ge.1) then write(*,990) flag write(10,990) flag end if go to 500 end if c if the maximum number of iterations was achieved then stop if (iter.ge.maxit) then flag= 8 if (iprint.ge.1) then write(*,998) flag write(10,998) flag end if go to 500 end if c ======================================================= c Compute direction c ======================================================= if (iter.eq.0) then do i= 1, nind d(i)= -r(i) end do dnorm2= rnorm2 dtr= -rnorm2 else beta= rnorm2/rnorm2prev do i= 1, nind d(i)= -r(i) + beta*d(i) end do dnorm2= rnorm2 - 2.0d0*beta*(dtr+alpha*dtw) + beta**2*dnorm2 dtr= -rnorm2 + beta*(dtr+alpha*dtw) end if c Force d to be a descent direction of q(s), i.e., c <\nabla q(s), d> = = \le 0. if (dtr.gt.0.0d0) then do i= 1, nind d(i)= -d(i) end do dtr= -dtr end if c ======================================================= c Compute d^T H d c ======================================================= c w= A d call calchd(nind,ind,x,d,g,n,x,gtype,htvtype,w,y,sterel,steabs, *goth,inform) if (inform.ne.0) then flag = 11 if (iprint.ge.1) then write(*,1011) flag write(10,1011) flag end if go to 500 end if c Compute d^T w and ||w||^2 dtw= 0.0d0 do i= 1, nind dtw= dtw + d(i)*w(i) end do c ======================================================= c Compute maximum step c ======================================================= c amax1 > 0 and amax1n < 0 are the values of alpha such c that ||s + alpha * d||_2 or ||s + alpha * d||_\infty = delta dts= 0.0d0 do i= 1, nind dts= dts + d(i)*s(i) end do c 2-norm trust-region radius if (trtype.eq.0) then aa= dnorm2 bb= 2.0d0*dts cc= snorm2-delta**2 dd= dsqrt(bb**2-4.0d0*aa*cc) amax1= (-bb+dd)/(2.0d0*aa) amax1n= (-bb-dd)/(2.0d0*aa) c infinite-norm trust-region radius else if (trtype.eq.1) then amax1= infty amax1n= -infty do i= 1, nind if (d(i).gt.0.0d0) then amax1= dmin1(amax1, (delta-s(i))/d(i)) amax1n= dmax1(amax1n, (-delta-s(i))/d(i)) else if (d(i).lt.0.0d0) then amax1= dmin1(amax1, (-delta-s(i))/d(i)) amax1n= dmax1(amax1n, (delta-s(i))/d(i)) end if end do end if c amax2 > 0 and amax2n < 0 are the maximum and the minimum c values of alpha such that l - x <= s + alpha * d <= u - x, c respectively amax2= infty amax2n= -infty do i= 1, nind if (d(i).gt.0.0d0) then if (u(i).lt.infty) then amax2x= (u(i)-x(i)-s(i))/d(i) if (amax2x.lt.amax2) then amax2= amax2x rbdposaind= i rbdposatype= 2 end if end if if (l(i).gt.-infty) then amax2nx= (l(i)-x(i)-s(i))/d(i) if (amax2nx.gt.amax2n) then amax2n= amax2nx rbdnegaind= i rbdnegatype= 1 end if end if else if (d(i).lt.0.0d0) then if (l(i).gt.-infty) then amax2x= (l(i)-x(i)-s(i))/d(i) if (amax2x.lt.amax2) then amax2= amax2x rbdposaind= i rbdposatype= 1 end if end if if (u(i).lt.infty) then amax2nx= (u(i)-x(i)-s(i))/d(i) if (amax2nx.gt.amax2n) then amax2n= amax2nx rbdnegaind= i rbdnegatype= 2 end if end if end if end do c Compute amax as the minimum among amax1 and amax2, and c amaxn as the minimum among amax1n and amax2n. Moreover c change amaxn by -amaxn to have amax and amaxn as maximum c steps along d direction (and not -d in the case of amaxn) amax= dmin1(amax1, amax2) amaxn= dmax1(amax1n, amax2n) c ======================================================= c Compute the step (and the quadratic functional value at c the new point) c ======================================================= qprev= q c If d^T H d > 0 then take the conjugate gradients step if (dtw.gt.0.0d0) then alpha= dmin1(amax, rnorm2/dtw) q= q + 0.5d0*alpha**2*dtw + alpha*dtr c If d^T H d <= 0 and function f is nearly quadratic then c take the point with the minimum functional value (q) among c the current one and the ones which are at the boundary, i.e., c the best one between q(s), q(s + amax*d) and q(s + amaxn*d). else qamax= q + 0.5d0*amax**2*dtw + amax*dtr c If we are at iteration zero then take the maximum positive c step in the minus gradient direction if (iter.eq.0) then alpha= amax q= qamax c If we are not in the first iteration then if function f is c nearly quadratic and q(s + amax * d) or q(s + amaxn * d) is c smaller than q(s), go to the best point in the boundary else if (nearlyq.and.(qamax.lt.q.or.qamaxn.lt.q)) then qamaxn= q + 0.5d0*amaxn**2*dtw + amaxn*dtr if (qamax.lt.qamaxn) then alpha= amax q= qamax else alpha= amaxn q= qamaxn end if c Else, stop at the current point else flag= 7 if (iprint.ge.1) then write(*,997) flag write(10,997) flag end if go to 500 end if end if c ======================================================= c Compute new s c ======================================================= do i= 1, nind sprev(i)= s(i) s(i)= s(i) + alpha*d(i) end do snorm2prev= snorm2 snorm2= snorm2 + alpha**2*dnorm2 + 2.0d0*alpha*dts c ======================================================= c Compute the residual r = H s + g c ======================================================= rnorm2prev= rnorm2 do i= 1, nind r(i)= r(i) + alpha*w(i) end do rnorm2= norm2s(nind,r) c ======================================================= c Increment number of iterations c ======================================================= iter= iter + 1 c ======================================================= c Print information of this iteration c ======================================================= if (iprint.ge.1) then write(*,984) iter,dsqrt(rnorm2),dsqrt(snorm2),q write(10,984) iter,dsqrt(rnorm2),dsqrt(snorm2),q end if c ======================================================= c Test other stopping criteria c ======================================================= c Test angle condition gts= 0.0d0 do i= 1, nind gts= gts + g(i)*s(i) end do if (gts.gt.0.0d0 .or. gts**2.lt.theta**2*gnorm2*snorm2) then do i= 1, nind s(i)= sprev(i) end do snorm2= snorm2prev q= qprev flag= 3 if (iprint.ge.1) then write(*,993) flag write(10,993) flag end if go to 500 end if c If we are in the boundary of the box also stop if (alpha.eq.amax2.or.alpha.eq.amax2n) then if (alpha.eq.amax2) then rbdind= rbdposaind rbdtype= rbdposatype else ! if (alpha.eq.amax2n) then rbdind= rbdnegaind rbdtype= rbdnegatype end if if (rbdtype.eq.1) then rbdtypea= 'lower' else ! if (rbdtype.eq.2) then rbdtypea= 'upper' end if flag= 2 if (iprint.ge.1) then write(*,992) flag,ind(rbdind),rbdtypea write(10,992) flag,ind(rbdind),rbdtypea end if go to 500 end if c If we are in the boundary of the trust region then stop if (alpha.eq.amax1.or.alpha.eq.amax1n) then flag= 1 if (iprint.ge.1) then write(*,991) flag write(10,991) flag end if go to 500 end if c If two consecutive iterates are much close then stop samep= .true. do i= 1, nind if (dabs(alpha*d(i)).gt.dmax1(epsrel*dabs(s(i)),epsabs)) then samep= .false. end if end do if (samep) then flag= 6 if (iprint.ge.1) then write(*,996) flag write(10,996) flag end if go to 500 end if c Test whether we performed many iterations without good progress c of the quadratic model c if (dabs(q-qprev).le.dmax1(epsrel*dabs(qprev),epsabs)) then c itnqmp= itnqmp + 1 c if (itnqmp.ge.maxitnqmp) then c flag= 4 c if (iprint.ge.1) then c write(*,994) flag,itnqmp c write(10,994) flag,itnqmp c end if c go to 500 c endif c else c itnqmp= 0 c endif c Test whether we performed many iterations without good progress c of the quadratic model currprog= qprev - q bestprog= dmax1(currprog, bestprog) if (currprog.le.epsnqmp*bestprog) then itnqmp= itnqmp + 1 if (itnqmp.ge.maxitnqmp) then flag= 4 if (iprint.ge.1) then write(*,994) flag,itnqmp,epsnqmp,bestprog write(10,994) flag,itnqmp,epsnqmp,bestprog end if go to 500 endif else itnqmp= 0 endif c ======================================================= c Iterate c ======================================================= go to 100 c ======================================================= c End of main loop c ======================================================= c ======================================================= c Return c ======================================================= 500 continue c Print final information if (iprint.ge.1) then write(*,985) min0(nind,ncomp),(s(i),i=1,min0(nind,ncomp)) write(10,985) min0(nind,ncomp),(s(i),i=1,min0(nind,ncomp)) end if return c Non-executable statements 980 format(/,6x,'Conjugate gradients (maxit= ',I7,' acc= ',1PD11.4, *')') 981 format(6x,'Using Euclidian trust region (delta= ',1PD11.4, *')') 982 format(6x,'Using sup-norm trust region (delta= ',1PD11.4,')') 983 format(6x,'Unknown trust-region type') 984 format(6x,'CG iter= ',I5,' rnorm: ',1PD11.4,' snorm= ',1PD11.4, *' q= ',1PD11.4) 985 format(/,6x,'Truncated Newton direction (first ',I6, *' components): ',/,1(6x,6(1PD11.4,1x))) 990 format(6x,'Flag of CG= ',I2,' (Convergence with small residual)') 991 format(6x,'Flag of CG= ',I2, *' (Convergence to the trust region boundary)') 992 format(6x,'Flag of CG= ',I2, *' (Convergence to the boundary of the box constraints,',/,6x, *'taking step >= 1, variable ',I7,' will reaches its ',A5, *' bound)') 993 format(6x,'Flag of CG= ',I2, *' (The next CG iterate will not satisfy the angle condition)') 994 format(6x,'Flag of CG= ',I2, *' (Not enough progress in the quadratic model. This means',/,6x, *'that the progress of the last ',I7,' iterations was smaller ', *'than ',/,6x,1PD11.4,' times the best progress (',1PD11.4,')') 996 format(6x,'Flag of CG= ',I2, *' (Very near consecutive iterates)') 997 format(6x,'Flag of CG= ',I2, *' (d such that d^T H d = 0 and g^T d = 0 was found)') 998 format(6x,'Flag of CG= ',I2,' (Too many GC iterations)') 1011 format(6x,'Flag of CG= ',I2,' (Error in evalhd subroutine)') end c ****************************************************** c ****************************************************** subroutine tnls(nind,ind,n,x,f,g,d,l,u,amax,rbdtype,rbdind,fmin, *maxfc,gtype,iprint,fcnt,gcnt,intcnt,exgcnt,exbcnt,flag,xplus, *xtemp,gamma,beta,sigma1,sigma2,nint,next,mininterp,sterel,steabs, *epsrel,epsabs,infty) integer nind,ind(nind),n,rbdtype,rbdind,maxfc,gtype,iprint,fcnt, *gcnt,intcnt,exgcnt,exbcnt,flag,mininterp double precision x(n),f,g(n),d(n),l(n),u(n),amax,fmin,xplus(n), *xtemp(n),gamma,beta,sigma1,sigma2,nint,next,sterel,steabs,epsrel, *epsabs,infty c This subrotuine implements the line search described in a working c paper by E. G. Birgin and J. M. Martinez. c c On Entry c c nind integer c number of free variables (this is thee dimension in c which this subroutine will work) c c ind integer ind(n) c array which contains, in the first nind positions, the c identifiers of the free variables c c n integer c dimension of the full space c c x double precision x(n) c current point c c The first nind positions of x contains the free variables c x_ind(1), x_ind(2), ..., x_ind(nind). c c f double precision c functional value at x c c g double precision g(n) c gradient vector at x c c It also contains in the first nind positions the components c g_ind(1), g_ind(2), ..., g_ind(nind). c c IMPORTANT: the linear algebra of this subroutine lies in c a space of dimension nind. The value of the full dimension n, c the non-free variables (which are at the end of array x) and c its gradient components (which are at the and of array g) c are also used and actualized any time the gradient is being c evaluated. c c d double precision d(nind) c descent direction c c l double precision l(nind) c lower bounds on x. It components are ordered in the c same way as x and g. c c u double precision u(nind) c upper bounds on x. It components are ordered in the c same way as x, g and l. c c fmin double precision c functional value for the stopping criteria f <= fmin c c maxfc integer c maximum number of funtion evaluations c c gtype integer c type of gradient calculation c gtype = 0 means user suplied evalg subroutine, c gtype = 1 means central diference approximation. c c RECOMMENDED: gtype = 0 c c (provided you have the evalg subroutine) c c output logical c TRUE: print some information at each iteration, c FALSE: no print. c c xplus, xtemp double precision xplus(nind),xtemp(nind) c working vectors c c gamma double precision c constant for the Armijo crtierion c f(x + alpha d) <= f(x) + gamma * alpha * <\nabla f(x),d> c c RECOMMENDED: gamma = 10^{-4} c c beta double precision c constant for the beta condition c < beta * c if (x_k + d_k) satisfies the Armijo condition but does not c satisfy the beta condition then the point is accepted, but c if it satisfied the Armijo condition and also satisfies the c beta condition then we know that there is the possibility c for a succesful extrapolation c c RECOMMENDED: beta = 0.5 c c sigma1, sigma2 double precision c constant for the safeguarded interpolation c if alpha_new \notin [sigma1, sigma*alpha] then we take c alpha_new = alpha / nint c c RECOMMENDED: sigma1 = 0.1 and sigma2 = 0.9 c c nint double precision c constant for the interpolation. See the description of c sigma1 and sigma2 above. Sometimes we take as a new trial c step the previous one divided by nint c c RECOMMENDED: nint = 2.0 c c next double precision c constant for the extrapolation c when extrapolating we try alpha_new = alpha * next c c RECOMMENDED: next = 2.0 c c mininterp integer c constant for testing if, after having made at least mininterp c interpolations, the steplength is soo small. In that case c failure of the line search is declared (may be the direction c is not a descent direction due to an error in the gradient c calculations) c c RECOMMENDED: mininterp = 4 c c sterel, steabs double precision c this constants mean a ``relative small number'' and ``an c absolute small number'' for the increments in finite c difference approximations of derivatives c c RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} c c epsrel, epsabs, infty double precision c this constants mean a ``relative small number'', ``an c absolute small number'', and ``infinite or a very big c number''. Basically, a quantity A is considered negligeble c with respect to another quantity B if c |A| < max ( epsrel * |B|, epsabs ) c c RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20} and c infty = 10^{+20} c c On Return c c x double precision x(n) c new current point c c f double precision c functional value at x c c g double precision g(n) c gradient vector at x c c fcnt integer c number of funtional evaluations used in this line search c c gcnt integer c number of gradient evaluations used in this line search c c intcnt integer c number of interpolations c c exgcnt integer c number of good extrapolations c c exbcnt integer c number of bad extrapolations c c flag integer c This output parameter tells what happened in this c subroutine, according to the following conventions: c c 0= convergence with an Armijo-like criterion c (f(xnew) <= f(x) + 1.0d-4 * alpha * ); c c 4= the algorithm stopped because the functional value c is very small (f <= fmin); c c 6= soo small step in the line search. After having made at c least mininterp interpolations, the steplength becames c small. ``small steplength'' means that we are at point c x with direction d and step alpha, and, for all i, c c |alpha * d(i)| .le. max ( epsrel * |x(i)|, epsabs ). c c In that case failure of the line search is declared c (may be the direction is not a descent direction c due to an error in the gradient calculations). Use c mininterp > maxfc for inhibit this criterion; c c 8= it was achieved the maximum allowed number of c function evaluations (maxfc); c c 9= error in evalf subroutine; c c 10= error in evalg subroutine. c Armijo parameter and beta parameter. You can modify them if c necessary logical samep integer i,interp,inform double precision fplus,atemp,ftemp,gtd,gptd,alpha,fbext c ======================================================= c Compute directional derivative c ======================================================= gtd= 0.0d0 do i= 1, nind gtd= gtd + g(i)*d(i) end do c ======================================================= c Compute first trial c ======================================================= alpha= dmin1(1.0d0, amax) do i= 1, nind xplus(i)= x(i) + alpha*d(i) end do if (alpha.eq.amax) then if (rbdtype.eq.1) then xplus(rbdind)= l(rbdind) else ! if (rbdtype.eq.2) then xplus(rbdind)= u(rbdind) end if end if call calcf(nind,ind,xplus,n,x,fplus,inform) fcnt= fcnt + 1 if (inform.ne.0) then flag = 9 if (iprint.ge.1) then write(*,1000) flag write(10,1000) flag end if go to 500 end if c Print initial information if (iprint.ge.1) then write(*,980) amax write(*,999) alpha,fplus,fcnt write(10,980) amax write(10,999) alpha,fplus,fcnt end if c ======================================================= c Test Armijo and beta-condition and decide for accepting c the trial point, interpolate or extrapolate. c ======================================================= if (amax.gt.1.0d0) then c x + d belongs to the interior of the feasible set (amax > 1) if (iprint.ge.1) then write(*,*) ' x+d belongs to int of the feasible set' write(10,*) ' x+d belongs to int of the feasible set' end if c Verify Armijo if (fplus.le.f+gamma*alpha*gtd) then c Armijo condition holds if (iprint.ge.1) then write(*,*) ' Armijo condition holds' write(10,*) ' Armijo condition holds' end if call calcg(nind,ind,xplus,n,x,gtype,g,sterel,steabs, * inform) gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.1) then write(*,1001) flag write(10,1001) flag end if go to 500 end if gptd= 0.0d0 do i= 1, nind gptd= gptd + g(i)*d(i) end do c Verify directional derivative (beta condition) if (gptd.ge.beta*gtd) then c Step = 1 was ok, finish the line search if (iprint.ge.1) then write(*,*) ' The beta condition is also true' write(*,*) ' Line search is over' write(10,*) ' The beta condition is also true' write(10,*) ' Line search is over' end if f= fplus do i= 1, nind x(i)= xplus(i) end do flag= 0 if (iprint.ge.1) then write(*,990) flag write(10,990) flag end if go to 500 else c Extrapolate if (iprint.ge.1) then write(*,*) ' The beta-condition does not hold' write(*,*) ' We will extrapolate' write(10,*)' The beta-condition does not hold' write(10,*)' We will extrapolate' end if c f before extrapolation fbext= fplus go to 100 end if else c Interpolate if (iprint.ge.1) then write(*,*) ' Armijo does not hold' write(*,*) ' We will interpolate' write(10,*) ' Armijo does not hold' write(10,*) ' We will interpolate' end if go to 200 end if else c x + d does not belong to the feasible set (amax <= 1) if (iprint.ge.1) then write(*,*) ' x+d does not belong to box-interior' write(10,*) ' x+d does not belong to box-interior' end if if (fplus.lt.f) then c Extrapolate if (iprint.ge.1) then write(*,*) ' f(x+d) < f(x)' write(*,*) ' We will extrapolate' write(10,*) ' f(x+d) < f(x)' write(10,*) ' We will extrapolate' end if c f before extrapolation fbext= fplus go to 100 else c Interpolate if (iprint.ge.1) then write(*,*) ' f(x+d) >= f(x)' write(*,*) ' We will interpolate' write(10,*) ' f(x+d) >= f(x)' write(10,*) ' We will interpolate' end if go to 200 end if end if c ======================================================= c Extrapolation c ======================================================= 100 continue c Test f going to -inf 120 if (fplus.le.fmin) then c Finish the extrapolation with the current point f= fplus do i= 1, nind x(i)= xplus(i) end do call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.1) then write(*,1001) flag write(10,1001) flag end if go to 500 end if if (f.lt.fbext) then exgcnt= exgcnt + 1 else exbcnt= exbcnt + 1 end if flag= 4 if (iprint.ge.1) then write(*,994) flag write(10,994) flag end if go to 500 end if c Test maximum number of functional evaluations if (fcnt.ge.maxfc) then c Finish the extrapolation with the current point f= fplus do i= 1, nind x(i)= xplus(i) end do call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.1) then write(*,1001) flag write(10,1001) flag end if go to 500 end if if (f.lt.fbext) then exgcnt= exgcnt + 1 else exbcnt= exbcnt + 1 end if flag= 8 if (iprint.ge.1) then write(*,998) flag write(10,998) flag end if go to 500 end if c Chose new step if (alpha.lt.amax.and.next*alpha.gt.amax) then atemp= amax else atemp= next*alpha end if c Compute new trial point do i= 1, nind xtemp(i)= x(i) + atemp*d(i) end do if (atemp.eq.amax) then if (rbdtype.eq.1) then xtemp(rbdind)= l(rbdind) else ! if (rbdtype.eq.2) then xtemp(rbdind)= u(rbdind) end if end if c Project if (atemp.gt.amax) then do i= 1, nind xtemp(i)= dmax1(l(i),dmin1(xtemp(i), u(i))) end do end if c Test if this is not the same point as the previous one. c This test is performed only when alpha >= amax. if(alpha.ge.amax) then samep= .true. do i= 1, nind if (dabs(xtemp(i)-xplus(i)).gt. * dmax1(epsrel*dabs(xplus(i)),epsabs)) then samep= .false. end if end do if (samep) then c Finish the extrapolation with the current point f= fplus do i= 1, nind x(i)= xplus(i) end do call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.1) then write(*,1001) flag write(10,1001) flag end if go to 500 end if if (f.lt.fbext) then exgcnt= exgcnt + 1 else exbcnt= exbcnt + 1 end if flag= 0 if (iprint.ge.1) then write(*,990) flag write(10,990) flag end if go to 500 end if end if c Evaluate function call calcf(nind,ind,xtemp,n,x,ftemp,inform) fcnt= fcnt + 1 if (inform.ne.0) then flag = 9 if (iprint.ge.1) then write(*,1000) flag write(10,1000) flag end if go to 500 end if c Print information of this iteration if (iprint.ge.1) then write(*,999) atemp,ftemp,fcnt write(10,999) atemp,ftemp,fcnt end if c If the functional value decreases then set the current c point and continue the extrapolation if (ftemp.lt.fplus) then alpha= atemp fplus= ftemp do i= 1, nind xplus(i)= xtemp(i) end do go to 120 c If the functional value does not decrease then discard the c last trial and finish the extrapolation with the previous point else f= fplus do i= 1, nind x(i)= xplus(i) end do call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.1) then write(*,1001) flag write(10,1001) flag end if go to 500 end if if (f.lt.fbext) then exgcnt= exgcnt + 1 else exbcnt= exbcnt + 1 end if flag= 0 if (iprint.ge.1) then write(*,990) flag write(10,990) flag end if go to 500 end if c ======================================================= c End of extrapolation c ======================================================= c ======================================================= c Interpolation c ======================================================= 200 continue intcnt= intcnt + 1 c Initialization interp= 0 c Test f going to -inf 210 if (fplus.le.fmin) then c Finish the interpolation with the current point f= fplus do i= 1, nind x(i)= xplus(i) end do call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.1) then write(*,1001) flag write(10,1001) flag end if go to 500 end if flag= 4 if (iprint.ge.1) then write(*,994) flag write(10,994) flag end if go to 500 end if c Test maximum number of functional evaluations if (fcnt.ge.maxfc) then c As this is an abrupt termination then the current point of c the interpolation may be worst than the intial one c If the current point is better than the initial one then c finish the interpolation with the current point else c discard all we did inside this line search and finish with c the initial point if (fplus.lt.f) then f= fplus do i= 1, nind x(i)= xplus(i) end do call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.1) then write(*,1001) flag write(10,1001) flag end if go to 500 end if end if flag= 8 if (iprint.ge.1) then write(*,998) flag write(10,998) flag end if go to 500 end if c Compute new step interp= interp + 1 if (alpha.lt.sigma1) then alpha= alpha/nint else atemp= (-gtd*alpha**2)/(2.0d0*(fplus-f-alpha*gtd)) if (atemp.lt.sigma1 .or. atemp.gt.sigma2*alpha) then alpha= alpha/nint else alpha= atemp end if end if c Compute new trial point do i= 1, nind xplus(i)= x(i) + alpha*d(i) end do call calcf(nind,ind,xplus,n,x,fplus,inform) fcnt= fcnt + 1 if (inform.ne.0) then flag = 9 if (iprint.ge.1) then write(*,1000) flag write(10,1000) flag end if go to 500 end if c Print information of this iteration if (iprint.ge.1) then write(*,999) alpha,fplus,fcnt write(10,999) alpha,fplus,fcnt end if c Test Armijo condition if (fplus.le.f+gamma*alpha*gtd) then c Finish the line search f= fplus do i= 1, nind x(i)= xplus(i) end do call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt= gcnt + 1 if (inform.ne.0) then flag = 10 if (iprint.ge.1) then write(*,1001) flag write(10,1001) flag end if go to 500 end if flag= 0 if (iprint.ge.1) then write(*,990) flag write(10,990) flag end if go to 500 end if c Test whether at least mininterp interpolations were made and c two consecutive iterates are much close samep= .true. do i= 1, nind if (dabs(alpha*d(i)).gt.dmax1(epsrel*dabs(x(i)),epsabs)) then samep= .false. end if end do if (interp.gt.mininterp.and.samep) then c As this is an abrupt termination then the current point of c the interpolation may be worst than the intial one c If the current point is better than the initial one then c finish the interpolation with the current point else c discard all we did inside this line search and finish with c the initial point c if (fplus.lt.f) then c f= fplus c do i= 1, nind c x(i)= xplus(i) c end do c call calcg(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) c gcnt= gcnt + 1 c if (inform.ne.0) then c flag = 10 c if (iprint.ge.1) then c write(*,1001) flag c write(10,1001) flag c end if c go to 500 c end if c end if c The previous lines were commented because, as it is been c used, this subroutine must return with the intial point c in case of finding a very small interpolation step flag= 6 if (iprint.ge.1) then write(*,996) flag write(10,996) flag end if go to 500 end if c Else, iterate go to 210 c ======================================================= c End of interpolation c ======================================================= 500 continue c ======================================================= c Return c ======================================================= return c Non-executable statements 980 format(/,6x,'TN Line search (alphamax= ',1PD11.4,')') 999 format(6x,'Alpha= ',1PD11.4,' F= ',1PD11.4,' FE= ',I5) 990 format(6x,'Flag of TN Line search= ',I2, *' (Convergence with an Armijo-like criterion)') 994 format(6x,'Flag of TN Line search= ',I2, *' (Small functional value, smaller than ',/,6X,'parameter fmin)') 996 format(6x,'Flag of TN Line search= ',I2, *' (Too small step in the interpolation process)') 998 format(6x,'Flag of TN Line search= ',I2, *' (Too many functional evaluations)') 1000 format(6x,'Flag of TN Line search= ',I2, *' (Error in evalf subroutine)') 1001 format(6x,'Flag of TN Line search= ',I2, *' (Error in evalg subroutine)') end c ****************************************************** c ****************************************************** subroutine Gshrink(nind,ind,n,x,g,l,u) integer n, nind, ind(nind) double precision x(n), g(n), l(n), u(n) c Shrink all vectors x, g, l and u from the full dimension c space (dimension n) to the reduced space (dimension nind). integer i, indi double precision temp c Shrink x, g, l and u for the reduced space do i= 1, nind indi= ind(i) if (i.ne.indi) then temp= x(indi) x(indi)= x(i) x(i)= temp temp= g(indi) g(indi)= g(i) g(i)= temp temp= l(indi) l(indi)= l(i) l(i)= temp temp= u(indi) u(indi)= u(i) u(i)= temp end if end do return end c ****************************************************** c ****************************************************** subroutine Gexpand(nind,ind,n,x,g,l,u) integer n, nind, ind(nind) double precision x(n), g(n), l(n), u(n) c Expands vectors x, g, l and u from the reduced space c (dimension nind) to the full space (dimension n). integer i, indi double precision temp c Expand x, g, l and u to the full space do i= nind, 1, -1 indi= ind(i) if (i.ne.indi) then temp= x(indi) x(indi)= x(i) x(i)= temp temp= g(indi) g(indi)= g(i) g(i)= temp temp= l(indi) l(indi)= l(i) l(i)= temp temp= u(indi) u(indi)= u(i) u(i)= temp end if end do return end c ****************************************************** c ****************************************************** subroutine calcf(nind,ind,x,n,xc,al,inform) integer nind,ind(nind),n,inform double precision x(n),xc(n),al c This subroutines computes the augmented lagrangean. It is c called from the reduced space (dimension nind), expands c the point x where the function will be evaluated and call c the subroutine evalf to compute the augmented lagrangean c function. Finally, returns vector x to the reduced space c (shrinking). c c About subroutines named calc[something]. The subroutines c whos names start with ``calc'' work in (are called from) c the reduced space. Their task is (i) expand the arguments c to the full space, (ii) call the corresponding ``eval'' c subroutine (which works in the full space), and (iii) c shrink the parameters again and also shrink an eventualy c output of the ``eval'' subroutines. Subroutines of this c type are: calcf, calcg, calchd and calchddiff (the last c one called from calchd) and the corresponding subroutines c in the full space are evalf, evalg and evalhd. integer i,indi double precision temp c Complete x do i= nind+1, n x(i)= xc(i) end do c Expand x to the full space do i= nind, 1, -1 indi= ind(i) if (i.ne.indi) then temp= x(indi) x(indi)= x(i) x(i)= temp end if end do c Compute f call evalq(n,x,al,inform) c Shrink x to the reduced space do i= 1, nind indi= ind(i) if (i.ne.indi) then temp= x(indi) x(indi)= x(i) x(i)= temp end if end do return end c ****************************************************** c ****************************************************** subroutine calcg(nind,ind,x,n,xc,gtype,nal,sterel,steabs,inform) integer nind,ind(nind),n,gtype,inform double precision x(n),xc(n),nal(n),sterel,steabs c This subroutine computes the augmented lagrangean gradient c vector nal(x). The way it is computed depends on parameter gtype: c c gtype = 0: a user defined function, called evalg, is used, c gtype = 1: central finite differences are used. c c It is called from the reduced space (dimension nind), expands c the point x where the gradient will be evaluated and call c the subroutine evalg to compute the augmented lagrangean c gradient vector. Finally, shrinks vectors x and nal to the c reduced space. c c About subroutines named calc[something]. The subroutines c whos names start with ``calc'' work in (are called from) c the reduced space. Their task is (i) expand the arguments c to the full space, (ii) call the corresponding ``eval'' c subroutine (which works in the full space), and (iii) c shrink the parameters again and also shrink an eventualy c output of the ``eval'' subroutines. Subroutines of this c type are: calcf, calcg, calchd and calchddiff (the last c one called from calchd) and the corresponding subroutines c in the full space are evalf, evalg and evalhd. integer i,indi double precision temp c Complete x do i= nind+1, n x(i)= xc(i) end do c Expand x to the full space do i= nind, 1, -1 indi= ind(i) if (i.ne.indi) then temp= x(indi) x(indi)= x(i) x(i)= temp end if end do c Compute gradient vector c True gradient if (gtype.eq.0) then call evalnq(n,x,nal,inform) c Central finite diferences approximation else if (gtype.eq.1) then call evalgdiff(n,x,nal,sterel,steabs,inform) end if c shrink x and g to the reduced space do i= 1, nind indi= ind(i) if (i.ne.indi) then temp= x(indi) x(indi)= x(i) x(i)= temp temp= nal(indi) nal(indi)= nal(i) nal(i)= temp end if end do return end c ****************************************************** c ****************************************************** subroutine calchd(nind,ind,x,d,g,n,xc,gtype,htvtype,hd,xtemp, *sterel,steabs,goth,inform) logical goth integer nind,ind(nind),n,gtype,htvtype,inform double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n),sterel,steabs c This subroutine computes the product Hessian times vector d. c The way it is computed depends on parameter htvtype: c c mvptype = 0: a user subroutine, called evalhd, for computing c the hessian times vector d product is called, c c mvptype = 1: a incremental quotients aproximation is used. c c If a user subroutine will be used then this subroutine (which c is called from the reduced space) expands vectors x and d, c calls the user supplied subroutine to compute the hessian times c vector d product, and shrinks vectors x, d and hd. c c About subroutines named calc[something]. The subroutines c whos names start with ``calc'' work in (are called from) c the reduced space. Their task is (i) expand the arguments c to the full space, (ii) call the corresponding ``eval'' c subroutine (which works in the full space), and (iii) c shrink the parameters again and also shrink an eventualy c output of the ``eval'' subroutines. Subroutines of this c type are: calcf, calcg, calchd and calchddiff (the last c one called from calchd) and the corresponding subroutines c in the full space are evalf, evalg and evalhd. integer i,indi double precision temp c ======================================================= c User defined subroutine to compute the Hessian times c vector d product (it works in the full dimension space) c ======================================================= if (htvtype.eq.0) then c Complete d with zeroes do i= nind+1, n d(i)= 0.0d0 end do c Complete x do i= nind+1, n x(i)= xc(i) end do c Expand x, d and g to the full space do i= nind, 1, -1 indi= ind(i) if (i.ne.indi) then temp= x(indi) x(indi)= x(i) x(i)= temp temp= d(indi) d(indi)= d(i) d(i)= temp temp= g(indi) g(indi)= g(i) g(i)= temp end if end do c Call the user defined subroutine call evalhd(nind,ind,x,d,g,n,xc,gtype,htvtype,hd,xtemp, * sterel,steabs,goth,inform) c Shrink x, d, g and hd to the reduced space do i= 1, nind indi= ind(i) if (i.ne.indi) then temp= x(indi) x(indi)= x(i) x(i)= temp temp= d(indi) d(indi)= d(i) d(i)= temp temp= g(indi) g(indi)= g(i) g(i)= temp temp= hd(indi) hd(indi)= hd(i) hd(i)= temp end if end do c ======================================================= c Incremental quotients approximation c (it works in the reduced space) c ======================================================= else if (htvtype.eq.1) then call calchddiff(nind,ind,x,d,g,n,xc,gtype,hd,xtemp,sterel, * steabs,inform) end if return end c **************************************************** c **************************************************** subroutine calchddiff(nind,ind,x,d,g,n,xc,gtype,hd,xtemp,sterel, *steabs,inform) integer n,nind,ind(n),gtype,inform double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n),sterel,steabs c This subroutine computes the Hessian times vector d product c by means of a ``directional finite difference''. The idea is c that, at the current point x, the product H d is the limit of c c [ Gradient(x + t d) - Gradient(x) ] / t c c In this implementation we use c c t = max(steabs, sterel ||x||_\infty) / ||d||_\infty c c provided that d not equal 0, of course. c c So, we evaluate the Gradient at the auxiliary point x + t d c and use the quotient above to approximate H d. c c About subroutines named calc[something]. The subroutines c whos names start with ``calc'' work in (are called from) c the reduced space. Their task is (i) expand the arguments c to the full space, (ii) call the corresponding ``eval'' c subroutine (which works in the full space), and (iii) c shrink the parameters again and also shrink an eventualy c output of the ``eval'' subroutines. Subroutines of this c type are: calcf, calcg, calchd and calchddiff (the last c one called from calchd) and the corresponding subroutines c in the full space are evalf, evalg and evalhd. c c On Entry c c n integer c order of the x c c x double precision x(n) c point for which Hessian(x) times d will be approximated c c d double precision d(n) c vector for which the Hessian times vetor product will c be approximated c c g double precision g(n) c gradient at x c c xtemp double precision xtemp(n) c working vector c c sterel, steabs double precision c this constants mean a ``relative small number'' and c ``an absolute small number'' c c On Return c c hd double precision hd(n) c approximation of H d integer i double precision xinfn,dinfn,step c Compute incremental quotients step xinfn= 0.0d0 dinfn= 0.0d0 do i= 1, nind xinfn= dmax1(xinfn, dabs(x(i))) dinfn= dmax1(dinfn, dabs(d(i))) end do step= dmax1(sterel*xinfn,steabs) / dinfn c Set the point in which the gradient will be c evaluated: xtemp = x + step * d do i= 1, nind xtemp(i)= x(i) + step*d(i) end do c Evaluate the gradient at (x + step * d) call calcg(nind,ind,xtemp,n,xc,gtype,hd,sterel,steabs,inform) c Compute incremental quotients do i= 1, nind hd(i)= (hd(i)-g(i)) / step end do return end c **************************************************** c **************************************************** subroutine evalgdiff(n,x,nal,sterel,steabs,inform) integer n,inform double precision x(n),nal(n),sterel,steabs c Computes the augmented lagrangean gradient vector by central c finite differences. This subroutine, which works in the full c space, is prepared to replace the subroutine evalg (to c evaluate the gradient vector) in the case of the lastest have c not being provided by the user. integer j double precision tmp,step,fplus,fminus do j= 1, n tmp= x(j) step= dmax1(steabs, sterel*dabs(tmp)) x(j)= tmp + step call evalq(n,x,fplus,inform) x(j)= tmp - step call evalq(n,x,fminus,inform) nal(j)= (fplus-fminus) / (2.0d0*step) x(j)= tmp end do return end c **************************************************** c **************************************************** double precision function norm2s(n,x) integer n double precision x(n) double precision hsldnrm2 norm2s= hsldnrm2(n,x,1) ** 2 return end c **************************************************** c **************************************************** DOUBLE PRECISION FUNCTION HSLDNRM2(N,DX,INCX) DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) DOUBLE PRECISION CUTLO,CUTHI PARAMETER (CUTLO=8.232D-11,CUTHI=1.304D19) INTEGER INCX,N DOUBLE PRECISION DX(*) DOUBLE PRECISION HITEST,SUM,XMAX INTEGER I,J,NEXT,NN INTRINSIC DABS,DSQRT,FLOAT IF (N.GT.0) GO TO 10 HSLDNRM2 = ZERO GO TO 300 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N*INCX I = 1 20 GO TO NEXT 30 IF (DABS(DX(I)).GT.CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO 50 IF (DX(I).EQ.ZERO) GO TO 200 IF (DABS(DX(I)).GT.CUTLO) GO TO 85 ASSIGN 70 TO NEXT GO TO 105 100 I = J ASSIGN 110 TO NEXT SUM = (SUM/DX(I))/DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 70 IF (DABS(DX(I)).GT.CUTLO) GO TO 75 110 IF (DABS(DX(I)).LE.XMAX) GO TO 115 SUM = ONE + SUM* (XMAX/DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 75 SUM = (SUM*XMAX)*XMAX 85 HITEST = CUTHI/DFLOAT(N) DO 95 J = I,NN,INCX IF (DABS(DX(J)).GE.HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 HSLDNRM2 = DSQRT(SUM) GO TO 300 200 CONTINUE I = I + INCX IF (I.LE.NN) GO TO 20 HSLDNRM2 = XMAX*DSQRT(SUM) 300 CONTINUE RETURN END c **************************************************** c **************************************************** c Modifications introduced from March 1st to March 21th of 2002 c in ocassion of the ISPG development: c c 1) Comments of some new parameters introduced in the previous c modification c c 2) As it was, in the first iteration of GENCAN (when kappa c takes value equal 1) and for one-dimensional faces, cgmaxit c (the maximum number of Conjugate Gradient iterations to compute c the internal to the face truncated-Newton direction) was being 0. c As it is obviously wrong, we add a max between what was being c computed and one to allow at least one CG iteration. c c 3) Parameter inform in subroutines evalf, evalg and evalhd c supplied by the user was added c c Modifications introduced from May 31th to November 2nd of 2001 c in ocassion of the ALGENCAN development: c c Fixed bugs: c c 1) The first spectral steplength was not been projected in the c [lammin,lammax] interval. c c 2) The conjugate gradients accuracy (cgeps) which is linearly c dependent of the euclidian-norm of the projected gradient, was c also not been projected in the interval [cgepsi,cgepsf]. c c 3) Conjugate gradients said that it was being used an euclidian c norm trust region when it has really being used an infinite norm c trust region and viceversa. c c 4) Sometimes, the analytic gradient has been used although the c user choose the finite differences option. c c Modifications: c c 1) To avoid roundoff errors, an explicit detection of at least one c variable reaching its bound when a maximum step is being made was c added. c c 2) The way in which two points were considered very similar in, for c example, the interpolations and the extrapolations (which was c dependent of the infinity norm of the points) showed to be very c scale dependent. A new version which test the difference coordinate c to coordinate was done. In this was the calculus of the current c point x and the descent direction sup-norm is not done any c more. c c 3) The same constants epsrel and epsabs were used as small c relative and absolute values for, for example, detecting similar c points and for finite diferences. Now, c epsrel and epsabs are used for detecting similar points c (and the recommended values are 10^{-10} and 10^{-20}, respectively) c and new constants sterel ans steabs were introduced for finite c differences (and the recommended values are 10^{-7} and 10^{-10}, c respectively). c c 4) Two new stopping criteria for CG were added: (i) we stop if c two consecutive iterates are too close; and (ii) we also c stop if there is no enough quadratic model progress during c maxitnqmp iterations. c c 5) The linear relation between the conjugate gradient accuracy c and the norm of the projected gradient can be computed using c the euclidian- and the sup-norm of the projected gradient (only c euclidian-norm version was present in the previous version. The c linear relation is such that the CG accuracy is cgepsi when the c projected gradient norm value is equal to the value corresponding c to the initial guess and the CG accuracy is cgepsf when the c projected gradient norm value is cgrelf). c c 6) Inside Conjugate Gradients, the euclidian-norm is been c computed using an algorithm developed by C.L.LAWSON, 1978 JAN 08. c Numerical experiments showed that the performance of GENCAN depends c basically on the conjugate gradients performance and stopping c criteria and that the conjugate gradients depends on the way the c euclidian-norm is been computed. These things deserve further c research. c c 7) In the Augmented Lagrangean algorithm ALGENCAN, which uses c GENCAN to solve the bounded constrained subproblems, the maximum c number of Conjugate Gradients iterations (cgmaxit), which in c this version is linearly dependent of the projected gradient norm, c was set to 2 * (# of free variables). As CG is not using restarts c we do not know very well what this means. On the other hand, c the accuracy (given by cgeps) continues being more strict when c we are near to the solution and less strict when we ar far from c the solution. c c 8) Many things in the output were changed. c **************************************************** c **************************************************** subroutine evalhd(nind,ind,y,dy,nl,m,yc,gtype,htvtype,hd,temp, *sterel,steabs,goth,inform) C Computing H d, where H = lambda A A^t C PARAMETERS integer mmax parameter (mmax = 100 000) integer nmax parameter (nmax = 4 000) integer anzemax parameter (anzemax = 3 * mmax) C SCALAR ARGUMENTS logical goth integer nind,m,gtype,htvtype,inform double precision sterel,steabs C ARRAY ARGUMENTS integer ind(nind) double precision y(m),dy(m),nl(m),yc(m),hd(m),temp(m) C LOCAL SCALARS integer i C LOCAL ARRAYS double precision Atd(nmax) C COMMON SCALARS integer anze,n double precision lambda C COMMON ARRAYS integer arow(anzemax),acol(anzemax) double precision g(nmax),bbar(mmax),aval(anzemax) C COMMON BLOCKS common /ldata1/ arow,acol,aval,anze,n common /ldata2/ g,bbar,lambda inform = 0 C lambda A A^t d C A^t d do i = 1,n Atd(i) = 0.0d0 end do do i = 1,anze Atd(acol(i)) = Atd(acol(i)) + aval(i) * dy(arow(i)) end do C A A^t d do i = 1,m hd(i) = 0.0d0 end do do i = 1,anze hd(arow(i)) = hd(arow(i)) + aval(i) * Atd(acol(i)) end do C lambda A A^t d do i = 1,m hd(i) = lambda * hd(i) end do end subroutine evalq(m,y,l,inform) C Computing -L(y) = 0.5 lambda || (gk + A^t y) ||_2^2 + bbar^t y C PARAMETERS integer mmax parameter (mmax = 100 000) integer nmax parameter (nmax = 4 000) integer anzemax parameter (anzemax = 3 * mmax) C SCALAR ARGUMENTS integer m,inform double precision l C ARRAY ARGUMENTS double precision y(m) C LOCAL SCALARS integer i C LOCAL ARRAYS double precision d(nmax) C COMMON SCALARS integer anze,n double precision lambda C COMMON ARRAYS integer arow(anzemax),acol(anzemax) double precision g(nmax),bbar(mmax),aval(anzemax) C COMMON BLOCKS common /ldata1/ arow,acol,aval,anze,n common /ldata2/ g,bbar,lambda inform = 0 C d = g + A^t y do i = 1,n d(i) = g(i) end do do i = 1,anze d(acol(i)) = d(acol(i)) + aval(i) * y(arow(i)) end do C l = 0.5 lambda d^t d + bbar^t y l = 0.0d0 do i = 1,n l = l + d(i) ** 2 end do l = 0.5d0 * lambda * l do i = 1,m l = l + bbar(i) * y(i) end do return end subroutine evalnq(m,y,nl,inform) C Computing - \nabla L(y) = lambda A (gk + A^t y) + bbar C PARAMETERS integer mmax parameter (mmax = 100 000) integer nmax parameter (nmax = 4 000) integer anzemax parameter (anzemax = 3 * mmax) C SCALAR ARGUMENTS integer m,inform C ARRAY ARGUMENTS double precision y(m),nl(m) C LOCAL SCALARS integer i C LOCAL ARRAYS double precision d(nmax) C COMMON SCALARS integer anze,n double precision lambda C COMMON ARRAYS integer arow(anzemax),acol(anzemax) double precision g(nmax),bbar(mmax),aval(anzemax) C COMMON BLOCKS common /ldata1/ arow,acol,aval,anze,n common /ldata2/ g,bbar,lambda inform = 0 C d = g + A^t y do i = 1,n d(i) = g(i) end do do i = 1,anze d(acol(i)) = d(acol(i)) + aval(i) * y(arow(i)) end do C nl = bbar + lambda A d do i = 1,m nl(i) = bbar(i) end do do i = 1,anze nl(arow(i)) = nl(arow(i)) + lambda * aval(i) * d(acol(i)) end do return end