C Last update: May 14th, 2003. C ALGENCAN main subroutine. This subroutine implements the 67 Augmented C Lagrangian methods for inequaltity constraints and bound-constrained C problems described in: C C E. G. Birgin, R. Castillo and J. M. Martínez, "Numerical comparison C of Augmented Lagrangian algorithms for nonconvex problems", C Computational Optimization and Applications 31, pp. 31-56, 2005. subroutine algencan(pname,n,x,l,u,m,c,mu,mubar,rho,zeta,epsfc, +maxit,maxfc,epsnqmp,maxitnqmp,al,nal,gpeucn2,gpsupn,constrv, +outiter,totiter,totfcnt,totgcnt,totcgit,totspgit,tottnit, +inform,nearlyq,epsgpen,epsgpsn,epsnfp,fmin,delta0,cgmaxit,cgscre, +cggpnf,cgepsi,cgepsf,maxitnfp,maxitngp,gtype,htvtype,trtype, +iprint,w1,w2,cprev,eta,delmin,lammax,lammin,theta,gamma,beta, +sigma1,sigma2,nint,next,mininterp,maxextrap,ncomp,sterel,steabs, +epsrel,epsabs,infrel,infabs) C SCALAR ARGUMENTS logical nearlyq character * 10 pname integer n,m,maxit,maxfc,maxitnqmp,outiter,totiter,totfcnt,totgcnt, + totcgit,totspgit,tottnit,inform,cgmaxit,maxitnfp,maxitngp, + gtype,htvtype,trtype,iprint,mininterp,maxextrap,ncomp, + cgscre double precision epsfc,epsnqmp,al,gpeucn2,gpsupn, + constrv,epsgpen,epsgpsn,epsnfp,fmin,delta0,cggpnf, + cgepsi,cgepsf,eta,delmin,lammax,lammin,theta,gamma,beta, + sigma1,sigma2,nint,next,sterel,steabs,epsrel,epsabs, + infrel,infabs C ARRAY ARGUMENTS integer w2(n) double precision x(n),l(n),u(n),c(m),mu(m),mubar(m),zeta(m), + rho(m),nal(n),w1(10*n+maxitngp),cprev(m) C Solves the general-constrained minimization problem C C Minimize f(x) C C subject to C C c(x) <= 0 C l <= x <= u C C using different Augmented Lagrangean methods (as part of a working C paper by E. G. Birgin, R. Castillo and J. M. Martinez). C inform integer C This output parameter tells what happened in this C subroutine, according to the following conventions: C C 0 = convergence with feasibility and complementarity; C C 1 = it was achieved the maximum allowed number of C outer iterations (maxoutit); C C 2 = it was achieved the maximum allowed total number of C inner iterations (maxtotit); C C 3 = it was achieved the maximum allowed total number of C functional evaluations (maxtotfc); C C 4 = the algorithm stopped by ``lack of feasibility progress'', C that means that c(x_k) > epsfc and c(x_k) >= c(x_{k-1}) C during maxitncp consecutive iterations; C C -90 = error in evalfc subroutine; C C -91 = error in evalgjc subroutine; C C -92 = error in evalp subroutine; C C -93 = error in evaldpdy subroutine; C C -94 = error in evalhd subroutine. C COMMON SCALAR integer ptype,ttype double precision tconst,umuini,urhoini,urhomul,urhofrac,umumin, + umumax C LOCAL SCALARS integer i,iter,fcnt,gcnt,cgcnt,spgiter,spgfcnt,maxtotit,maxtotfc, + tniter,tnfcnt,tnstpcnt,tnintcnt,tnexgcnt,tnexbcnt,tnintfe, + tnexgfe,tnexbfe,maxoutit,nprint,mprint,icrit,maxitncp,flag double precision constrvprev,complv,f,v,vprev,uw,uwprev,muprev C logical interior C COMMON BLOCKS common /p/ ptype,ttype,tconst common /user/umuini,urhoini,urhomul,urhofrac,umumin,umumax C CONSTANTS C just for printing nprint = min0(n,ncomp) mprint = min0(m,ncomp) maxoutit = 12 maxitncp = maxoutit maxtotit = 100000 maxtotfc = 100000 C INTIALIZATION icrit = 0 outiter = 0 totiter = 0 totspgit = 0 tottnit = 0 totfcnt = 0 totgcnt = 0 totcgit = 0 if (iprint.ge.1) then write(*, 900) ptype,ttype,pname,n,m,epsfc,epsfc,epsgpsn write(10,900) ptype,ttype,pname,n,m,epsfc,epsfc,epsgpsn end if C COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS call evalfc(n,x,f,m,c,flag) if ( flag .ne. 0 ) then ! ERROR IN THE OBJECTIVE FUNCTION OR THE CONSTRAINTS ! EVALUATION, STOP inform = -90 if (iprint.ge.1) then write(*, 1000) inform write(10,1000) inform end if go to 500 end if C SET INITIAL LAGRANGEAN MULTIPLIERS APPROXIMATION, INITIAL PENALTY C PARAMETER AND ZETA do i = 1,m mu(i) = umuini mubar(i) = max( umumin, min( umumax, mu(i) ) ) rho(i) = urhoini if ( ptype .eq. 9 .and. ttype .ge. 16 .and. ttype .le. 18 ) + then rho(i) = max( rho(i), 2.0d0 * mubar(i) ) end if zeta(i) = rho(i) end do C COMPUTE MAXIMUM FEASIBILITY VIOLATION AND COMPLEMENTARITY VIOLATION constrv = 0.0d0 complv = 0.0d0 do i = 1,m constrv = max( constrv, c(i) ) if ( c(i) .lt. -epsfc .and. mu(i) .gt. epsfc ) then complv = max( complv, - mu(i) * c(i) ) end if end do C PRINT INITIAL INFORMATION if (iprint.ge.1) then write(*,910) outiter,nprint,(x(i),i=1,nprint) write(*,920) f,constrv,complv,mprint,(mu(i),i=1,mprint) write(10,910) outiter,nprint,(x(i),i=1,nprint) write(10,920) f,constrv,complv,mprint,(mu(i),i=1,mprint) end if C MAIN LOOP 100 continue C TEST STOPPING CRITERIA if (outiter .ge. maxoutit) then ! MAXIMUM NUMBER OF OUTER ITERATIONS EXCEEDED, STOP inform = 1 if (iprint.ge.1) then write(*, 961) inform,maxoutit write(10,961) inform,maxoutit end if go to 500 end if if (totiter .gt. maxtotit) then ! MAXIMUM NUMBER OF INNER ITERATIONS EXCEEDED, STOP inform = 2 if (iprint.ge.1) then write(*, 962) inform,maxtotit write(10,962) inform,maxtotit end if go to 500 end if if (totfcnt .gt. maxtotfc) then ! MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED, STOP inform = 3 if (iprint.ge.1) then write(*, 963) inform,maxtotfc write(10,963) inform,maxtotfc end if go to 500 end if C DO AN OUTER ITERATION outiter= outiter + 1 if (iprint.ge.1) then write(*, 940) outiter write(10,940) outiter end if C CALL THE INNER-SOLVER if (iprint.ge.1) then write(*, 950) epsgpsn write(10,950) epsgpsn end if call gencan(n,x,l,u,epsgpen,epsgpsn,epsnfp,maxitnfp,maxitngp,fmin, +maxit,maxfc,delta0,cgmaxit,cgscre,cggpnf,cgepsi,cgepsf, +epsnqmp,maxitnqmp,nearlyq,gtype,htvtype,trtype,iprint,al,nal, +gpeucn2,gpsupn,iter,fcnt,gcnt,cgcnt,spgiter,spgfcnt,tniter,tnfcnt, +tnstpcnt,tnintcnt,tnexgcnt,tnexbcnt,tnintfe,tnexgfe,tnexbfe, +inform,w1(1),w1(n+1),w1(2*n+1),w2(1),w1(3*n+1),w1(3*n+maxitngp+1), +eta,delmin,lammax,lammin,theta,gamma,beta,sigma1,sigma2,nint,next, +mininterp,maxextrap,ncomp,sterel,steabs,epsrel,epsabs,infrel, +infabs) totiter = totiter + iter totspgit = totspgit + spgiter tottnit = tottnit + tniter totfcnt = totfcnt + fcnt totgcnt = totgcnt + gcnt totcgit = totcgit + cgcnt if ( inform .lt. 0 ) then if (iprint.ge.1) then write(*, 1000) inform write(10,1000) inform end if go to 500 end if C COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS do i = 1,m cprev(i) = c(i) end do call evalfc(n,x,f,m,c,flag) if ( flag .ne. 0 ) then ! ERROR IN THE OBJECTIVE FUNCTION OR THE CONSTRAINTS ! EVALUATION, STOP inform = -90 if (iprint.ge.1) then write(*, 1000) inform write(10,1000) inform end if go to 500 end if C COMPUTE NEW LAGRANGEAN MULTIPLIERS APPROXIMATION, UPDATE PENALTY C PARAMETERS AND ZETA do i = 1,m muprev = mu(i) call evaldpdy(c(i),rho(i),mubar(i),zeta(i),mu(i),inform) mubar(i) = max( umumin, min( umumax, mu(i) ) ) if ( inform .lt. 0 ) then if (iprint.ge.1) then write(*, 1000) inform write(10,1000) inform end if go to 500 end if v = max( 0.0d0, c(i) ) vprev = max( 0.0d0, cprev(i) ) if ( c(i) .lt. -epsfc .and. mu(i) .gt. epsfc ) then uw = - mu(i) * c(i) else uw = 0.0d0 end if if ( cprev(i) .lt. -epsfc .and. muprev .gt. epsfc ) then uwprev = - muprev * cprev(i) else uwprev = 0.0d0 end if if ( v .gt. urhofrac * vprev .or. uw .gt. urhofrac * uwprev ) + then rho(i) = urhomul * rho(i) end if if ( ptype .eq. 9 .and. ttype .ge. 16 .and. ttype .le. 18 ) + then rho(i) = max( rho(i), 2.0d0 * mubar(i) ) end if zeta(i) = rho(i) end do C COMPUTE MAXIMUM FEASIBILITY VIOLATION AND COMPLEMENTARITY VIOLATION constrvprev= constrv constrv = 0.0d0 complv = 0.0d0 do i = 1,m constrv = max( constrv, c(i) ) if ( c(i) .lt. -epsfc .and. mu(i) .gt. epsfc ) then complv = max( complv, - mu(i) * c(i) ) end if end do C PRINT ITERATION INFORMATION if (iprint.ge.1) then write(*, 910) outiter,nprint,(x(i),i=1,nprint) write(*, 920) f,constrv,complv,mprint,(mu(i),i=1,mprint) write(*, 925) mprint,(rho(i),i=1,mprint) write(*, 930) al,gpsupn,totiter,totfcnt,totgcnt,totcgit write(10,910) outiter,nprint,(x(i),i=1,nprint) write(10,920) f,constrv,complv,mprint,(mu(i),i=1,mprint) write(10,925) mprint,(rho(i),i=1,mprint) write(10,930) al,gpsupn,totiter,totfcnt,totgcnt,totcgit end if C TEST FEASIBILITY AND COMPLEMENTARITY if ( constrv .le. epsfc .and. complv .le. epsfc ) then ! THE POINT IS FEASIBLE AND COMPLEMENTARITY IS SATISFIED flag= 0 if (iprint.ge.1) then write(*, 960) flag,epsfc,epsfc,gpsupn,epsgpsn write(10,960) flag,epsfc,epsfc,gpsupn,epsgpsn end if go to 500 end if C TEST IF THE INFEASIBILITY DECREASED SUFFICIENTLY if ( constrv .gt. epsfc .and. constrv .ge. constrvprev ) then icrit = icrit + 1 if ( icrit .ge. maxitncp) then ! LACK OF FEASIBILITY PROGRESS, STOP flag = 4 if ( iprint .ge. 1 ) then write(*, 964) flag,maxitncp write(10,964) flag,maxitncp end if go to 500 end if else icrit = 0 end if C ITERATE go to 100 C STOP 500 continue return 900 format( /,' ****** Output of the Augmented Lagrangian', + ' defined by P = ',I1,' Theta = ',I2,' ******'/, + /,' Problem: ',A10,/, + /,' Main characteristics of the problem: ',/, + /,' Number of variables = ',4X,I7, + /,' Number of inequality constraints = ',4X,I7, + /,' Required constraints infeasibility = ',1PD11.4, + /,' Required complementarity violation = ',1PD11.4, + /,' Required projected gradient sup-norm = ',1PD11.4) 910 format( /,' OUTER ITERATION ',I11,/, + /,' Current point (first ',I6, ' components): ', + /,6(1X,1PD11.4)) 920 format( /,' Objective function value = ',1PD11.4, + /,' Constraints infeasibility = ',1PD11.4, + /,' Complementarity violation = ',1PD11.4,/, + /,' Current Lagrange multipliers estimates', + ' (first ',I6, ' components): ',/,6(1X,1PD11.4)) 925 format( /,' Current penalty parameters (first ',I6, + ' components): ',/,6(1X,1PD11.4)) 930 format( /,' Augmented lagrangian value = ',1PD11.4, + /,' Projected gradient sup-norm = ',1PD11.4,/, + /,' No. of inner iterations = ',I11, + /,' No. of functional evaluations = ',I11, + /,' No. of gradient evaluations = ',I11, + /,' No. of CG iterations = ',I11) 940 format( /,' Computing Outer Iteration ',I2,' ...') 950 format( /,' MINIMIZATION OF THE AUGMENTED LAGRANGIAN ',/, + /,' Required gradient sup-norm',1PD11.4) 960 format( /,' Flag of ALGENCAN = ',I3,/, + /,' Solution was found !!! ',/, + ' The infeasibility of the constraints is smaller', + ' than ',1PD11.4,' and the',/,' complementarity', + ' violation is smaller than ',1PD11.4,'. The', + ' sup-norm of',/,' the projected continuous gradient', + ' is ',1PD11.4,' (and something smaller',/, + ' than ',1PD11.4,' was requested).') 961 format( /,' Flag of ALGENCAN = ',I3,/, + /,' It was exceeded the maximum allowed number of outer', + ' iterations',/,' (maxoutit = ',I7,')') 962 format( /,' Flag of ALGENCAN = ',I3,/, + /,' It was exceeded the maximum allowed total number of', + ' inner iterations',/,' (maxtotit = ',I7,')') 963 format( /,' Flag of ALGENCAN = ',I3,/, + /,' It was exceeded the maximum allowed total number of', + ' functional evaluations',/,' (maxtotfc = ',I7,')') 964 format( /,' Flag of ALGENCAN = ',I3,/, + /,' Constraint violations have not decreased', + ' substantially over',I7,' major',/,' iterations.') 1000 format( /,' Flag of ALGENCAN = ',I3,' Fatal Error',/, + /,' The following codes means: ',/, + /,' -90 : error in evalfc subroutine', + /,' -91 : error in evalgjc subroutine', + /,' -92 : error in evalp subroutine', + /,' -93 : error in evaldpdy subroutine', + /,' -94 : error in evalhd subroutine',/) end c **************************************************** c **************************************************** c Last update: March 22th, 2002 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,inform,s,y,d,ind,lastgpns,w,eta,delmin,lammax,lammin, *theta,gamma,beta,sigma1,sigma2,nint,next,mininterp,maxextrap, *ncomp,sterel,steabs,epsrel,epsabs,infrel,infabs) C SCALAR ARGUMENTS 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,inform,mininterp, + maxextrap,ncomp double precision epsgpen,epsgpsn,epsnfp,fmin,udelta0,cggpnf, + cgepsi,cgepsf,epsnqmp,f,gpeucn2,gpsupn,eta,delmin,lammax, + lammin,theta,gamma,beta,sigma1,sigma2,nint,next,sterel, + steabs,epsrel,epsabs,infrel,infabs C ARRAY ARGUMENTS integer ind(n) double precision x(n),l(n),u(n),g(n),s(n),y(n),d(n), + lastgpns(0:maxitngp-1),w(5*n) 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", to appear in Computational Optimization and c Applications, 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 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 maxextrap integer c constant to limit the number of extrapolations in the truncated c Newton direction c c RECOMMENDED: mininterp = maxfc 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, infrel, infabs 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}, c infrel = 10^{+20} and infabs = 10^{+99} 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 inform 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 < 0 = error in evalal or evalnal subroutines. C LOCAL SCALARS character * 3 ittype integer i,itnfp,nind,nprint,cgmaxit,cgiter,fcntprev, + tnintprev,tnexgprev,tnexbprev,rbdtype,rbdind,infotmp double precision lambda,gieucn2,gpi,fprev,sts,sty,epsgpen2,delta, + currprog,bestprog,cgeps,xnorm,ometa2,amax,amaxx,acgeps, + bcgeps,kappa,gpeucn20,gpsupn0,gpnmax 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, and fprev= infabs bestprog= 0.0d0 itnfp= 0 c for testing progress in the projected gradient norm. do i = 0,maxitngp-1 lastgpns(i) = infabs end do c Print problem information if(iprint.ge.2) 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 evalal(n,x,f,inform) fcnt= fcnt + 1 if ( inform .lt. 0 ) then if (iprint.ge.2) then write(*, 1000) inform write(10,1000) inform end if return end if if (gtype.eq.0) then call evalnal(n,x,g,inform) else if (gtype.eq.1) then call evalnaldiff(n,x,g,sterel,steabs,inform) end if gcnt= gcnt + 1 if ( inform .lt. 0 ) then if (iprint.ge.2) then write(*,1000) inform write(10,1000) inform end if return 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.2) 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 inform= 4 if (iprint.ge.2) then write(*,994) inform,fmin write(10,994) inform,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 inform= 8 if (iprint.ge.2) then write(*,998) inform,maxfc write(10,998) inform,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 inform= 0 if (iprint.ge.2) then write(*,990) inform,epsgpen write(10,990) inform,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 inform= 1 if (iprint.ge.2) then write(*,991) inform,epsgpsn write(10,991) inform,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 inform= 2 if (iprint.ge.2) then write(*,992) inform,epsnfp,maxitnfp write(10,992) inform,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 gpnmax = 0.0d0 do i = 0,maxitngp-1 gpnmax = max( gpnmax, lastgpns(i) ) end do lastgpns(mod( iter,maxitngp )) = gpeucn2 if ( gpeucn2 .ge. gpnmax ) then inform = 3 if ( iprint .ge. 2 ) then write(*, 993) inform,maxitngp write(10,993) inform,maxitngp end if go to 500 endif c Test whether the number of iterations is exhausted if (iter.ge.maxit) then inform= 7 if (iprint.ge.2) then write(*,997) inform,maxit write(10,997) inform,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,inform, * w(1),w(n+1),gamma,sigma1,sigma2,nint,mininterp,sterel,steabs, * epsrel,epsabs,infrel,infabs) spgfcnt= spgfcnt + (fcnt-fcntprev) if ( inform .lt. 0 ) then if (iprint.ge.2) then write(*, 1000) inform write(10,1000) inform end if return end if c Compute the gradient of the new iterate if (gtype.eq.0) then call evalnal(n,x,g,inform) else if (gtype.eq.1) then call evalnaldiff(n,x,g,sterel,steabs,inform) end if gcnt= gcnt + 1 if ( inform .lt. 0 ) then if (iprint.ge.2) then write(*, 1000) inform write(10,1000) inform end if return 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 c delta= dmax1(delmin,100.0d0*dmax1(1.0d0,xnorm)) delta= dmax1(delmin,10.0d0*dmax1(1.0d0,xnorm)) else delta= udelta0 end if else c delta= dmax1(delmin,1.0d0*dsqrt(sts)) delta= dmax1(delmin,10.0d0*dsqrt(sts)) end if c Shrink the point, its gradient and the bounds call shrink(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,inform,w(1),w(n+1),w(2*n+1),w(3*n+1),w(4*n+1), * theta,sterel,steabs,epsrel,epsabs,infrel,infabs) cgcnt= cgcnt + cgiter if ( inform .lt. 0 ) then if (iprint.ge.2) then write(*, 1000) inform write(10,1000) inform end if return end if c Compute maximum step if ( inform .eq. 2 ) then amax = 1.0d0 else amax = infabs do i = 1,nind c if (d(i).gt.0.0d0.and.u(i).lt.infrel) then if ( d(i) .gt. 0.0d0 ) then amaxx = ( u(i) - x(i) ) / d(i) if ( amaxx .lt. amax ) then amax = amaxx rbdind = i rbdtype = 2 end if c else if (d(i).lt.0.0d0.and.l(i).gt.-infrel) then else if ( d(i) .lt. 0.0d0 ) 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, + inform,w(1),w(n+1),w(2*n+1),gamma,beta,sigma1,sigma2,nint, + next,mininterp,maxextrap,sterel,steabs,epsrel,epsabs,infrel, + infabs) if ( inform .lt. 0 ) then if (iprint.ge.2) then write(*, 1000) inform write(10,1000) inform end if return 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 expand(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 ( inform .eq. 6 ) then if (iprint.ge.2) 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, * inform,w(1),w(n+1),gamma,sigma1,sigma2,nint,mininterp, * sterel,steabs,epsrel,epsabs,infrel,infabs) spgfcnt= spgfcnt + (fcnt-fcntprev) if ( inform .lt. 0 ) then if (iprint.ge.2) then write(*, 1000) inform write(10,1000) inform end if return end if c Compute the gradient in the new iterate infotmp = inform if (gtype.eq.0) then call evalnal(n,x,g,inform) else if (gtype.eq.1) then call evalnaldiff(n,x,g,sterel,steabs,inform) end if gcnt= gcnt + 1 if ( inform .lt. 0 ) then if (iprint.ge.2) then write(*,1000) inform write(10,1000) inform end if return end if inform = infotmp 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.2) 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 ( inform .eq. 4 ) then if (iprint.ge.2) then write(*, 994) inform,fmin write(10,994) inform,fmin end if go to 500 else if ( inform .eq. 6 ) then if (iprint.ge.2) then write(*, 996) inform,mininterp,epsrel,epsabs write(10,996) inform,mininterp,epsrel,epsabs end if go to 500 else if ( inform .eq. 8 ) then if (iprint.ge.2) then write(*, 998) inform,maxfc write(10,998) inform,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.2) 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 = ',I3, *' (convergence with Euclidian-norm of the projected gradient', */,1X,'smaller than ',1PD11.4,')') 991 format(/1X,'Flag of GENCAN = ',I3, *' (convergence with sup-norm of the projected gradient', */,1X,'smaller than ',1PD11.4,')') 992 format(/1X,'Flag of GENCAN= ',I3, *' (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 = ',I3, *' (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 = ',I3, *' (The algorithm stopped because the functional value is', */,1X,'smaller than ',1PD11.4) 996 format(/1X,'Flag of GENCAN = ',I3, *' (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 = ',I3, *' (It was exceeded the maximum allowed number of iterations', */,1X,'(maxit=',I7,')') 998 format(/1X,'Flag of GENCAN = ',I3, *' (It was exceeded the maximum allowed number of functional', */,1X,'evaluations (maxfc=',I7,')') 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) 1000 format(/1X,'Flag of GENCAN = ',I3,' Fatal Error') end c ****************************************************** c ****************************************************** subroutine spgls(n,x,f,g,l,u,lambda,fmin,maxfc,iprint,fcnt,inform, *xtrial,d,gamma,sigma1,sigma2,nint,mininterp,sterel,steabs,epsrel, *epsabs,infrel,infabs) C SCALAR ARGUMENTS integer n,maxfc,iprint,fcnt,inform,mininterp double precision f,lambda,fmin,gamma,sigma1,sigma2,nint,sterel, + steabs,epsrel,epsabs,infrel,infabs C ARRAY ARGUMENTS double precision x(n),g(n),l(n),u(n),xtrial(n),d(n) c Safeguarded quadratic interpolation, used in Spectral Projected c Gradient directions. 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, infrel, infabs 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}, c infrel = 10^{+20} and infabs = 10^{+99} 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 inform 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 < 0 = error in evalal subroutine. C LOCAL SCALARS logical samep integer i,interp double precision ftrial,gtd,alpha,atemp c Print presentation information if (iprint.ge.3) then write(*, 980) lambda write(10,980) lambda end if 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 evalal(n,xtrial,ftrial,inform) fcnt= fcnt + 1 if ( inform .lt. 0 ) then if (iprint.ge.3) then write(*, 1000) inform write(10,1000) inform end if return end if c Print information of the first trial if (iprint.ge.3) then write(*, 999) alpha,ftrial,fcnt 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 inform = 0 if (iprint.ge.3) then write(*, 990) inform write(10,990) inform 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 inform = 4 if (iprint.ge.3) then write(*, 994) inform write(10,994) inform 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 inform = 8 if (iprint.ge.3) then write(*, 998) inform write(10,998) inform 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 evalal(n,xtrial,ftrial,inform) fcnt= fcnt + 1 if ( inform .lt. 0 ) then if (iprint.ge.3) then write(*, 1000) inform write(10,1000) inform end if return end if c Print information of the current trial if (iprint.ge.3) 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 inform = 6 if (iprint.ge.3) then write(*, 996) inform write(10,996) inform 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 = ',I3, * ' (Convergence with an Armijo-like criterion)') 994 format(6x,'Flag of SPG Line search = ',I3, * ' (Small functional value, smaller than ',/, * 6X,'parameter fmin)') 996 format(6x,'Flag of SPG Line search = ',I3, * ' (Too small step in the interpolation process)') 998 format(6x,'Flag of SPG Line search = ',I3, * ' (Too many functional evaluations)') 1000 format(6x,'Flag of SPG Line search = ',I3,' Fatal Error') 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,inform,w,y,r,d,sprev,theta,sterel,steabs,epsrel,epsabs, *infrel,infabs) C SCALAR ARGUMENTS logical nearlyq integer nind,n,maxitnqmp,maxit,gtype,htvtype,trtype,iprint,ncomp, + iter,rbdtype,rbdind,inform double precision delta,eps,epsnqmp,theta,sterel,steabs,epsrel, + epsabs,infrel,infabs C ARRAY ARGUMENTS integer ind(nind) double precision x(n),g(n),l(n),u(n),s(n),w(n),y(n),r(n),d(n), + sprev(n) 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, infrel, infabs 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}, c infrel = 10^{+20} and infabs = 10^{+99} 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 inform 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 < 0 = error in evalhd subroutine. C LOCAL SCALARS character * 5 rbdtypea logical goth,samep integer i,itnqmp,rbdposatype,rbdnegatype,rbdposaind,rbdnegaind 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 = infabs 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. 3 ) 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 inform = 0 if ( iprint .ge. 3 ) then write(*, 990) inform write(10,990) inform end if go to 500 end if c if the maximum number of iterations was achieved then stop if ( iter .ge. maxit ) then inform = 8 if ( iprint .ge. 3 ) then write(*, 998) inform write(10,998) inform 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 .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return 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 = sqrt( 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 = infabs amax1n = -infabs do i = 1,nind if ( d(i) .gt. 0.0d0 ) then amax1 = min( amax1, ( delta - s(i) ) / d(i) ) amax1n = max( amax1n, ( - delta - s(i) ) / d(i) ) else if ( d(i) .lt. 0.0d0 ) then amax1 = min( amax1, ( - delta - s(i) ) / d(i) ) amax1n = max( 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 = infabs amax2n = - infabs do i = 1,nind if ( d(i) .gt. 0.0d0 ) then c if (u(i).lt.infrel) then amax2x = ( u(i) - x(i) - s(i) ) / d(i) if ( amax2x .lt. amax2 ) then amax2 = amax2x rbdposaind = i rbdposatype = 2 end if c end if c if (l(i).gt.-infrel) then amax2nx = ( l(i) - x(i) - s(i) ) / d(i) if ( amax2nx .gt. amax2n ) then amax2n = amax2nx rbdnegaind = i rbdnegatype = 1 end if c end if else if ( d(i) .lt. 0.0d0 ) then c if (l(i).gt.-infrel) then amax2x = ( l(i) - x(i) - s(i) ) / d(i) if ( amax2x .lt. amax2 ) then amax2 = amax2x rbdposaind = i rbdposatype = 1 end if c end if c if (u(i).lt.infrel) then amax2nx = ( u(i) - x(i) - s(i) ) / d(i) if ( amax2nx .gt. amax2n ) then amax2n = amax2nx rbdnegaind = i rbdnegatype = 2 end if c 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 = min( amax1 , amax2 ) amaxn = max( 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 = min( 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 qamaxn = q + 0.5d0 * amaxn ** 2 * dtw + amaxn * dtr if ( nearlyq .and. + ( qamax .lt. q .or. qamaxn .lt. q ) ) then if ( qamax .lt. qamaxn ) then alpha = amax q = qamax else alpha = amaxn q = qamaxn end if c Else, stop at the current point else inform = 7 if ( iprint .ge. 3 ) then write(*, 997) inform write(10,997) inform end if go to 500 end if 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. 3 ) 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 inform = 3 if ( iprint .ge. 3 ) then write(*, 993) inform write(10,993) inform 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 C HABIA! c ACA HAY UN ERROR SI ALPHA = 1.0D+20 Y AMAX2 TAMBIEN ENTONCES c NO ES VERDAD QUE ESTAMOS EN EL BORDE DE LA CAJA Y RBDIND, QUE c SERA USADO COMO INDICE DE VETOR EN BREVE, NO ESTA DEFINIDO !!! c PARA RESOLVERLO TUVE QUE USAR 1.0D+99 Y -1.0D+99 PARA INICIALIZAR c EL CALCULO DE LOS AMAX, EN LUGAR DE 1.0D+20 Y -1.0D+20. c OTRA SOLUCION ES QUE EL RADIO DE LA REGION DE CONFIANZA NO SEA c NUNCA MAYOR QUE 1.0D+20. 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 inform = 2 if ( iprint .ge. 3 ) then write(*, 992) inform,ind(rbdind),rbdtypea write(10,992) inform,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 inform = 1 if ( iprint .ge. 3 ) then write(*, 991) inform write(10,991) inform end if go to 500 end if c If two consecutive iterates are much close then stop samep = .true. do i = 1,nind if ( abs( alpha * d(i) ) .gt. + max( epsrel * abs( s(i) ), epsabs ) ) then samep = .false. end if end do if ( samep ) then inform = 6 if ( iprint .ge. 3 ) then write(*, 996) inform write(10,996) inform 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 inform = 4 c if (iprint.ge.3) then c write(*,994) inform,itnqmp c write(10,994) inform,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 = max( currprog, bestprog ) if ( currprog .le. epsnqmp * bestprog ) then itnqmp = itnqmp + 1 if ( itnqmp .ge. maxitnqmp ) then inform = 4 if ( iprint .ge. 3 ) then write(*, 994) inform,itnqmp,epsnqmp,bestprog write(10,994) inform,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. 3 ) 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 = ',I3,' (Convergence with small residual)') 991 format(6x,'Flag of CG = ',I3, *' (Convergence to the trust region boundary)') 992 format(6x,'Flag of CG = ',I3, *' (Convergence to the boundary of the box constraints,',/,6x, *'taking step >= 1, variable ',I6,' will reaches its ',A5, *' bound)') 993 format(6x,'Flag of CG = ',I3, *' (The next CG iterate will not satisfy the angle condition)') 994 format(6x,'Flag of CG = ',I3, *' (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 = ',I3, *' (Very near consecutive iterates)') 997 format(6x,'Flag of CG= ',I3, *' (d such that d^T H d = 0 and g^T d = 0 was found)') 998 format(6x,'Flag of CG = ',I3,' (Too many GC iterations)') 1000 format(6x,'Flag of CG = ',I3,' Fatal Error') 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,inform,xplus, +xtemp,xbext,gamma,beta,sigma1,sigma2,nint,next,mininterp, +maxextrap,sterel,steabs,epsrel,epsabs,infrel,infabs) C SCALAR ARGUMENTS integer nind,n,rbdtype,rbdind,maxfc,gtype,iprint,fcnt,gcnt,intcnt, + exgcnt,exbcnt,inform,mininterp,maxextrap double precision f,amax,fmin,gamma,beta,sigma1,sigma2,nint,next, + sterel,steabs,epsrel,epsabs,infrel,infabs C ARRAY ARGUMENTS integer ind(nind) double precision x(n),g(n),d(n),l(n),u(n),xplus(n),xtemp(n), + xbext(n) 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, xbext double precision xplus(nind),xtemp(nind),xbext(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 maxextrap integer c constant to limit the number of extrapolations c c RECOMMENDED: maxextrap = 1000 (a big number) 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, infrel, infabs 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}, c infrel = 10^{+20} and infabs = 10^{+99} 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 inform 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 < 0 = error in evalal or evalnal subroutines. c Armijo parameter and beta parameter. You can modify them if c necessary C LOCAL SCALARS logical samep integer i,interp,extrap 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 = min( 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 calcal(nind,ind,xplus,n,x,fplus,inform) fcnt = fcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if c Print initial information if ( iprint .ge. 3 ) 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 if ( iprint .ge. 3 ) 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. 3 ) then write(*, *) ' Armijo condition holds' write(10,*) ' Armijo condition holds' end if call calcnal(nind,ind,xplus,n,x,gtype,g,sterel,steabs, * inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return 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. 3 ) 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 inform = 0 if ( iprint .ge. 3 ) then write(*, 990) inform write(10,990) inform end if go to 500 else c Extrapolate if ( iprint .ge. 3 ) 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 and x before extrapolation fbext = fplus do i = 1,nind xbext(i) = xplus(i) end do go to 100 end if else c Interpolate if ( iprint .ge. 3 ) 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. 3 ) 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. 3 ) 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 and x before extrapolation fbext = fplus do i = 1,nind xbext(i) = xplus(i) end do go to 100 else c Interpolate if ( iprint .ge. 3 ) 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 extrap = 0 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 if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then call calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if if ( f .lt. fbext ) then exgcnt = exgcnt + 1 else exbcnt = exbcnt + 1 end if end if inform = 4 if ( iprint .ge.3 ) then write(*, 994) inform write(10,994) inform 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 c If extrap=0 and amax>1 the gradient was computed for testing c the beta condition and it is not necessary to compute it again if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then call calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if if ( f .lt. fbext ) then exgcnt = exgcnt + 1 else exbcnt = exbcnt + 1 end if end if inform = 8 if ( iprint .ge. 3 ) then write(*, 998) inform write(10,998) inform end if go to 500 end if c Test if the maximum number of extrapolations was exceeded if ( extrap .ge. maxextrap ) then c Finish the extrapolation with the current point f = fplus do i = 1,nind x(i) = xplus(i) end do c If extrap=0 and amax>1 the gradient was computed for testing c the beta condition and it is not necessary to compute it again if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then call calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if if ( f .lt. fbext ) then exgcnt = exgcnt + 1 else exbcnt = exbcnt + 1 end if end if inform = 7 if ( iprint .ge. 3 ) then write(*, 997) inform write(10,997) inform 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) = max( l(i), min( 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 .gt. amax ) then samep = .true. do i = 1,nind if ( abs( xtemp(i) - xplus(i) ) .gt. + max( epsrel * abs( 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 c If extrap=0 and amax>1 the gradient was computed for testing c the beta condition and it is not necessary to compute it again if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then call calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs, + inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if if ( f .lt. fbext ) then exgcnt = exgcnt + 1 else exbcnt = exbcnt + 1 end if end if inform = 0 if ( iprint .ge. 3 ) then write(*, 990) inform write(10,990) inform end if go to 500 end if end if c Evaluate function call calcal(nind,ind,xtemp,n,x,ftemp,inform) fcnt = fcnt + 1 if ( inform .lt. 0 ) then c if ( iprint .ge. 3 ) then c write(*, 1000) inform c write(10,1000) inform c end if c return c If the objective function is not well defined in an c extrapolated point, we discard all the extrapolated points c and return to a safe region (where the point before c starting the extrapolations is) f = fbext do i = 1,nind x(i) = xbext(i) end do c If extrap=0 and amax>1 the gradient was computed for testing c the beta condition and it is not necessary to compute it again if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then call calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if exbcnt = exbcnt + 1 end if inform = 0 if ( iprint .ge. 3 ) then write(*, 1010) inform write(10,1010) inform end if go to 500 end if c Print information of this iteration if ( iprint .ge. 3 ) 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 extrap = extrap + 1 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 c If extrap=0 and amax>1 the gradient was computed for testing c the beta condition and it is not necessary to compute it again if ( extrap .ne. 0 .or. amax .le. 1.0d0 ) then call calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if if ( f .lt. fbext ) then exgcnt = exgcnt + 1 else exbcnt = exbcnt + 1 end if end if inform = 0 if ( iprint .ge.3 ) then write(*, 990) inform write(10,990) inform 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 calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if inform = 4 if ( iprint .ge. 3 ) then write(*, 994) inform write(10,994) inform 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 calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if end if inform = 8 if ( iprint .ge. 3 ) then write(*, 998) inform write(10,998) inform 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 calcal(nind,ind,xplus,n,x,fplus,inform) fcnt = fcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if c Print information of this iteration if ( iprint .ge. 3 ) 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 calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if inform = 0 if ( iprint .ge. 3 ) then write(*, 990) inform write(10,990) inform 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 ( abs( alpha * d(i) ) .gt. + max( epsrel * abs( 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 calcnal(nind,ind,x,n,x,gtype,g,sterel,steabs,inform) c gcnt = gcnt + 1 c if ( inform .lt. 0 ) then c if ( iprint .ge. 3 ) then c write(*, 1000) inform c write(10,1000) inform c end if c return 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. From c that intial point, something different will be tried. inform = 6 if ( iprint .ge. 3 ) then write(*, 996) inform write(10,996) inform 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= ',I3, + ' (Convergence with an Armijo-like criterion)') 994 format(6X,'Flag of TN Line search= ',I3, + ' (Small functional value, smaller than ',/, + 6X,'parameter fmin)') 996 format(6X,'Flag of TN Line search= ',I3, + ' (Too small step in the interpolation process)') 997 format(6X,'Flag of TN Line search= ',I3, + ' (Too many extrapolations)') 998 format(6X,'Flag of TN Line search= ',I3, + ' (Too many functional evaluations)') 1000 format(6X,'Flag of TN Line search = ',I3,' Fatal Error') 1010 format(6X,'Flag of TN Line search= ',I3, + ' (Fatal Error in an extrapolated point)') end c ****************************************************** c ****************************************************** subroutine shrink(nind,ind,n,x,g,l,u) C SCALAR ARGUMENTS integer n,nind C ARRAY ARGUMENTS integer 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). C LOCAL SCALARS 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 expand(nind,ind,n,x,g,l,u) C SCALAR ARGUMENTS integer n, nind C ARRAY ARGUMENTS integer 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). C LOCAL SCALARS 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 calcal(nind,ind,x,n,xc,al,inform) C SCALAR ARGUMENTS integer nind,n,inform double precision al C ARRAY ARGUMENTS integer ind(nind) double precision x(n),xc(n) 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 evalal 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: calcal, calcnal, 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 LOCAL SCALARS 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 evalal(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 calcnal(nind,ind,x,n,xc,gtype,nal,sterel,steabs,inform) C SCALAR ARGUMENTS integer nind,n,gtype,inform double precision sterel,steabs C ARRAY ARGUMENTS integer ind(nind) double precision x(n),xc(n),nal(n) 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 evalnal, 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 evalnal 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: calcal, calcnal, 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 LOCAL SCALARS 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 evalnal(n,x,nal,inform) c Central finite diferences approximation else if (gtype.eq.1) then call evalnaldiff(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) C SCALAR ARGUMENTS logical goth integer nind,n,gtype,htvtype,inform double precision sterel,steabs C ARRAY ARGUMENTS integer ind(nind) double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n) 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: calcal, calcnal, 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 LOCAL SCALARS 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) C SCALAR ARGUMENTS integer n,nind,gtype,inform double precision sterel,steabs C ARRAY ARGUMENTS integer ind(n) double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n) 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: calcal, calcnal, 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 C LOCAL SCALARS 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 calcnal(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 evalal(n,x,al,inform) C SCALAR ARGUMENTS double precision al integer n,inform C ARRAY ARGUMENTS double precision x(n) C This subroutine computes the augmented lagrangian function C C On Entry: C C n integer, C size of the problem, C C x double precision x(n), C point at which the function will be evaluated. C C On Return C C al double precision, C augmented lagrangean value at x, C C inform integer, C termination parameter: C 0= the function was succesfully evaluated, C 1= some error occurs in the function evaluation. C PARAMETERS integer mmax parameter (mmax=100000) C COMMON SCALARS integer m C COMMON ARRAYS double precision mubar(mmax),zeta(mmax),rho(mmax) C LOCAL SCALARS integer i,flag double precision f,p C LOCAL ARRAYS double precision c(mmax) C COMMON BLOCKS common /auglag/mubar,zeta,rho,m C INITIALIZE inform = 0 C COMPUTE OBJECTIVE FUNTION AND CONSTRAINTS call evalfc(n,x,f,m,c,flag) if ( flag .ne. 0 ) then inform = -90 return end if C SAVE THE CONSTRAINTS FOR THE GRADIENT EVALUATION call savec(n,x,m,c) C COMPUTES AL(x) = f(x) + \sum_i=1^m P(c(i),rho(i),mu(i)) al = f do i = 1,m call evalp(c(i),rho(i),mubar(i),zeta(i),p,inform) if ( inform .lt. 0 ) then return end if al = al + p end do end c **************************************************** c **************************************************** subroutine evalnal(n,x,nal,inform) C SCALAR ARGUMENTS integer n,inform C ARRAY ARGUMENTS double precision nal(n),x(n) C This subroutine computes the gradient of the C augmented lagrangian function. C C On Entry: C C n integer, C size of the problem, C C x double precision x(n), C point at which the gradient will be evaluated. C C On Return C C nal double precision nal(n), C gradient of the augmented lagrangean at x, C C inform integer, C termination parameter: C 0= the gradient was succesfully evaluated, C 1= some error occurs in the gradient evaluation. C PARAMETERS integer mmax,jcnzemax parameter (mmax=100000,jcnzemax=100000) C COMMON SCALARS integer m C COMMON ARRAYS double precision mubar(mmax),zeta(mmax),rho(mmax) C LOCAL SCALARS logical got integer i,jcnze,flag double precision f C LOCAL ARRAYS integer jcvar(jcnzemax),jcfun(jcnzemax) double precision c(mmax),dpdc(mmax),jcval(jcnzemax) C COMMON BLOCKS common /auglag/mubar,zeta,rho,m C INITIALIZE inform = 0 C COMPUTE CONSTRAINTS call gotc(n,x,m,c,got) if ( .not. got ) then call evalfc(n,x,f,m,c,flag) if ( flag .ne. 0 ) then inform = -90 return end if end if C COMPUTE dP/dc do i = 1,m call evaldpdy(c(i),rho(i),mubar(i),zeta(i),dpdc(i),inform) if ( inform .lt. 0 ) then return end if end do C COMPUTE \nabla L(x) = \nabla f(x) + \sum_{i=1}^m dPdc * dcdx do i = 1,n nal(i) = 0.0d0 end do call evalgjc(n,x,m,jcvar,jcfun,jcval,jcnze,flag) if ( flag .ne. 0 ) then inform = -91 return end if do i = 1,jcnze if ( jcfun(i).eq.0 ) then nal(jcvar(i)) = nal(jcvar(i)) + jcval(i) else nal(jcvar(i)) = + nal(jcvar(i)) + dpdc(jcfun(i)) * jcval(i) end if end do end c **************************************************** c **************************************************** subroutine evalnaldiff(n,x,nal,sterel,steabs,inform) C SCALAR ARGUMENTS integer n,inform double precision sterel,steabs C ARRAY ARGUMENTS double precision x(n),nal(n) 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 evalnal (to c evaluate the gradient vector) in the case of the lastest have c not being provided by the user. C LOCAL SCALARS integer j double precision tmp,step,fplus,fminus inform = 0 do j= 1, n tmp= x(j) step= dmax1(steabs, sterel*dabs(tmp)) x(j)= tmp + step call evalal(n,x,fplus,inform) if ( inform .lt. 0 ) then return end if x(j)= tmp - step call evalal(n,x,fminus,inform) if ( inform .lt. 0 ) then return end if nal(j)= (fplus-fminus) / (2.0d0*step) x(j)= tmp end do return end c **************************************************** c **************************************************** subroutine evalp(y,rho,mu,zeta,p,inform) C SCALAR ARGUMENTS integer inform double precision y,rho,mu,zeta,p C COMMON SCALARS integer ptype,ttype double precision tconst C LOCAL SCALARS double precision t,v,tmp,ytildv,thetav,thetav2,c C COMMON BLOCKS common /p/ ptype,ttype,tconst inform = 0 if ( ptype .eq. 0 ) then c This PHR but replacing max(0,a) = 0.5 ( a + |a| ) c and approximating |a| by sqrt( a^2 + 1/rho ). if ( ttype .eq. 1 ) then c = 1.0d0 else if ( ttype .eq. 2 ) then c = 1.0d-1 else if ( ttype .eq. 3 ) then c = 1.0d-2 else if ( ttype .eq. 4 ) then c = 1.0d-3 else if ( ttype .eq. 5 ) then c = 1.0d-4 else if ( ttype .eq. 6 ) then c = 1.0d-5 else if ( ttype .eq. 7 ) then c = 1.0d-6 else if ( ttype .eq. 8 ) then c = 0.0d0 end if tmp = mu + rho * y p = ( 0.5d0 / rho ) * + ( 0.25d0 * ( tmp + sqrt( tmp ** 2 + c / rho ) ) ** 2 - + mu ** 2 ) else if ( ptype .eq. 1 ) then if ( y .ge. 0.0d0 ) then p = mu * y + 0.5d0 * rho * y ** 2 + rho ** 2 * y ** 3 else if ( y .ge. -mu / rho ) then p = mu * y + 0.5d0 * rho * y ** 2 else p = -0.5d0 * mu ** 2 / rho end if else if ( ptype .eq. 2 ) then if ( y .ge. 0.0d0 ) then p = mu * y + mu * rho * y ** 2 + rho ** 2 * y ** 3 / 6.0d0 else p = mu * y / ( 1.0d0 - rho * y ) end if else if ( ptype .eq. 3 ) then if ( y .ge. 0.0d0 ) then p = mu * y + mu * rho * y ** 2 else p = ( mu * y ) / ( 1.0d0 - rho * y ) end if else if ( ptype .eq. 4 .and. ttype .eq. 1 ) then if ( mu + rho * y .ge. 0.0d0 ) then p = mu * y + 0.5d0 * rho * y ** 2 else p = - 0.5d0 * mu ** 2 / rho end if else if ( ptype .eq. 4 .and. ttype .eq. 2 ) then tmp = sign( 1.0d0, y ) * abs( y ) ** ( 1.0d0 / 3.0d0 ) * + rho ** ( 1.0d0 / 3.0d0 ) + mu if ( tmp .ge. 0.0d0 ) then p = (3.0d0 * tmp + mu) * y / 4.0d0 else p = - 0.25d0 * mu ** 4 / rho end if else if ( ptype .eq. 4 .and. ttype .eq. 3 ) then if ( mu + sinh( rho * y ) .ge. 0.0d0 ) then p = mu * y + ( cosh( rho * y ) - 1.0d0 ) / rho else tmp = log ( -mu + sqrt( mu ** 2 + 1.0d0 ) ) p = mu * tmp / rho + ( cosh( tmp ) - 1.0d0 ) / rho end if else if ((ptype .eq. 5 .or. ptype .eq. 8) .and. ttype .eq. 4) then t = rho * y v = exp( t ) - 1.0d0 p = mu * tconst * v / rho else if ( ptype .eq. 5 .and. ttype .eq. 5 ) then t = rho * y if ( t .le. 0.5d0 ) then v = - log( 1.0d0 - t ) p = mu * tconst * v / rho else v = exp( 2.0d0 * t - 1.0d0 ) + log(2.0d0) - 1.0d0 p = mu * tconst * v / rho end if else if ( ptype .eq. 5 .and. ttype .eq. 6 ) then t = rho * y if ( t .le. 0.5d0 ) then v = - log( 1.0d0 - t ) p = mu * tconst * v / rho else v = 2.0d0 * rho * y ** 2 + ( log( 2.0d0 ) - 0.5d0 ) / rho p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 7 ) then t = rho * y if ( t .le. 0.5d0 ) then v = y / ( 1.0d0 - t ) p = mu * tconst * v else v = exp ( 4.0d0 * t - 2.0d0 ) p = mu * tconst * v / rho end if else if ( ptype .eq. 5 .and. ttype .eq. 8 ) then t = rho * y if ( t .le. 0.5d0 ) then v = y / ( 1.0d0 - t ) p = mu * tconst * v else v = 8.0d0 * rho * y ** 2 - 4.0d0 * y + 1.0d0 / rho p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 9 ) then t = rho * y if ( t .le. -0.5d0 ) then v = - 0.25d0 * log ( -2.0d0 * t ) - 0.375d0 p = mu * tconst * v / rho else v = 0.5d0 * rho * y ** 2 + y p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 10 ) then t = rho * y tmp = ( 1.0d0 + t ) + sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 ) v = 0.0625d0 * tmp ** 2 + log ( 0.25d0 * tmp ) - 1.0d0 p = mu * tconst * v / rho else if ( ptype .eq. 5 .and. ttype .eq. 11 ) then t = rho * y if ( t .le. 0.5d0 ) then v = max( 0.0d0, t + 0.5d0 ) ** 3 / 6.0d0 - + 1.0d0 / 24.0d0 p = mu * tconst * v / rho else v = 0.5d0 * rho * y ** 2 p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 12 ) then t = rho * y if ( t .le. 0.5d0 ) then v = exp ( t ) p = mu * tconst * v / rho else v = 0.5d0 * exp( 0.5d0 ) * + ( rho * y ** 2 + y + 1.25d0 / rho ) p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 13 ) then t = rho * y if ( t .le. -0.5d0 ) then v = - log ( - t ) - 1.0d0 p = mu * tconst * v / rho else v = 2.0d0 * rho * y ** 2 + 4.0d0 * y + + ( 0.5d0 + log ( 2.0d0 ) ) / rho p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 14 ) then t = rho * y if ( t .le. -0.5d0 ) then v = - 1.0d0 / t p = mu * tconst * v / rho else v = 8.0d0 * rho * y ** 2 + 12.0d0 * y + 6.0d0 / rho p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 15 ) then t = rho * y if ( t .le. -1.0d0 ) then v = 4.0d0 / ( 1.0d0 - t ) - 2.0d0 p = mu * tconst * v / rho else if ( t .le. -0.25d0 ) then v = - log ( - t ) p = mu * tconst * v / rho else v = 8.0d0 * rho * y ** 2 + 8.0d0 * y + * ( 1.5d0 + log ( 4.0d0 ) ) / rho p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 16 ) then t = rho * y v = 0.5d0 * ( y + sqrt( y ** 2 + 4.0d0 / rho ** 2 ) ) p = mu * tconst * v else if ( ptype .eq. 5 .and. ttype .eq. 17 ) then t = rho * y if ( t .le. 0.0d0 ) then v = log ( 1.0d0 + exp ( t ) ) p = mu * tconst * v / rho else v = y + log ( 1.0d0 + exp ( - t ) ) / rho p = mu * tconst * v end if else if ( ptype .eq. 5 .and. ttype .eq. 18 ) then t = rho * y if ( t .le. 0.0d0 ) then v = 0.5d0 * exp ( t ) p = mu * tconst * v / rho else v = y + 0.5d0 * exp ( - t ) / rho p = mu * tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 4 ) then t = rho * mu * y v = exp( t ) - 1.0d0 p = tconst * v / rho else if ( ptype .eq. 6 .and. ttype .eq. 5 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then v = -log( 1.0d0 - t ) p = tconst * v / rho else v = exp( 2.0d0 * t - 1.0d0 ) + log( 2.0d0 ) - 1.0d0 p = tconst * v / rho end if else if ( ptype .eq. 6 .and. ttype .eq. 6 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then v = -log( 1.0d0 - t ) p = tconst * v / rho else v = 2.0d0 * rho * ( mu * y ) ** 2 + + ( log( 2.0d0 ) - 0.5d0 ) / rho p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 7 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then v = mu * y / ( 1.0d0 - t ) p = tconst * v else v = exp ( 4.0d0 * t - 2.0d0 ) p = tconst * v / rho end if else if ( ptype .eq. 6 .and. ttype .eq. 8 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then v = mu * y / ( 1.0d0 - t ) p = tconst * v else v = 8.0d0 * rho * ( mu * y ) ** 2 - 4.0d0 * mu * y + + 1.0d0 / rho p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 9 ) then t = rho * mu * y if ( t .le. -0.5d0 ) then v = -0.25d0 * log ( -2.0d0 * t ) - 0.375d0 p = tconst * v / rho else v = 0.5d0 * rho * ( mu * y ) ** 2 + mu * y p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 10 ) then t = rho * mu * y tmp = ( 1.0d0 + t ) + sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 ) v = 0.0625d0 * tmp ** 2 + log ( 0.25d0 * tmp ) - 1.0d0 p = tconst * v / rho else if ( ptype .eq. 6 .and. ttype .eq. 11 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then v = max( 0.0d0, t + 0.5d0 ) ** 3 / 6.0d0 - 1.0d0 / 24.0d0 p = tconst * v / rho else v = 0.5d0 * rho * ( mu * y ) ** 2 p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 12 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then v = exp( t ) p = tconst * v / rho else v = exp( 0.5d0 ) * ( 0.5d0 * rho * ( mu * y ) ** 2 + + 0.5d0 * mu * y + 0.625d0 / rho ) p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 13 ) then t = rho * mu * y if ( t .le. -0.5d0 ) then v = - log ( - t ) - 1.0d0 p = tconst * v / rho else v = 2.0d0 * rho * ( mu * y ) ** 2 + 4.0d0 * mu * y + + ( 0.5d0 + log ( 2.0d0 ) ) / rho p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 14 ) then t = rho * mu * y if ( t .le. -0.5d0 ) then v = - 1.0d0 / t p = tconst * v / rho else v = 8.0d0 * rho * ( mu * y ) ** 2 + 12.0d0 * mu * y + + 6.0d0 / rho p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 15 ) then t = rho * mu * y if ( t .le. -1.0d0 ) then v = 4.0d0 / ( 1.0d0 - t ) - 2.0d0 p = tconst * v / rho else if ( t .le. -0.25d0 ) then v = - log ( -t ) p = tconst * v / rho else v = 8.0d0 * mu * y * ( rho * mu * y + 1.0d0 ) + + ( 1.5d0 + log ( 4.0d0 ) ) / rho p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 16 ) then t = rho * mu * y v = 0.5d0 * ( mu * y + sqrt( ( mu * y ) ** 2 + + 4.0d0 / rho ** 2 ) ) p = tconst * v else if ( ptype .eq. 6 .and. ttype .eq. 17 ) then t = rho * mu * y if ( t .le. 0.0d0 ) then v = log ( 1.0d0 + exp ( t ) ) p = tconst * v / rho else v = mu * y + log ( 1.0d0 + exp ( -t ) ) / rho p = tconst * v end if else if ( ptype .eq. 6 .and. ttype .eq. 18 ) then t = rho * mu * y if ( t .le. 0.0d0 ) then v = 0.5d0 * exp ( t ) p = tconst * v / rho else v = mu * y + 0.5d0 * exp ( -t ) / rho p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 4 ) then t = rho * y / mu v = exp( t ) - 1.0d0 p = mu ** 2 * tconst * v / rho else if ( ptype .eq. 7 .and. ttype .eq. 5 ) then if ( rho * y .le. 0.5d0 * mu ) then t = rho * y / mu v = - log( 1.0d0 - t ) p = mu ** 2 * tconst * v / rho else t = rho * y / mu v = exp( 2.0d0 * t - 1.0d0 ) + log(2.0d0) - 1.0d0 p = mu ** 2 * tconst * v / rho end if else if ( ptype .eq. 7 .and. ttype .eq. 6 ) then if ( rho * y .lt. 0.5d0 * mu ) then t = rho * y / mu v = - log( 1.0d0 - t ) p = mu ** 2 * tconst * v / rho else v = 2.0d0 * rho * y ** 2 + + ( log( 2.0d0 ) - 0.5d0 ) * mu ** 2 / rho p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 7 ) then if ( rho * y .le. 0.5d0 * mu ) then v = mu ** 2 * y / ( mu - rho * y ) p = tconst * v else t = rho * y / mu v = exp ( 4.0d0 * t - 2.0d0 ) p = mu ** 2 * tconst * v / rho end if else if ( ptype .eq. 7 .and. ttype .eq. 8 ) then if ( rho * y .lt. 0.5d0 * mu ) then v = mu ** 2 * y / ( mu - rho * y ) p = tconst * v else v = 8.0d0 * rho * y ** 2 - 4.0d0 * y * mu + mu ** 2 / rho p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 9 ) then if ( rho * y .lt. - 0.5d0 * mu ) then t = rho * y / mu v = - 0.25d0 * log ( -2.0d0 * t ) - 0.375d0 p = mu ** 2 * tconst * v / rho else v = 0.5d0 * rho * y ** 2 + y * mu p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 10 ) then t = rho * y / mu tmp = ( 1.0d0 + t ) + sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 ) v = 0.0625d0 * tmp ** 2 + log ( 0.25d0 * tmp ) - 1.0d0 p = mu ** 2 * tconst * v / rho else if ( ptype .eq. 7 .and. ttype .eq. 11 ) then if ( rho * y .lt. 0.5d0 * mu ) then t = rho * y / mu v = max( 0.0d0, t + 0.5d0 ) ** 3 / 6.0d0 - 1.0d0 / 24.0d0 p = mu ** 2 * tconst * v / rho else v = 0.5d0 * rho * y ** 2 p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 12 ) then if ( rho * y .lt. 0.5d0 * mu ) then t = rho * y / mu v = exp( t ) p = mu ** 2 * tconst * v / rho else v = exp( 0.5d0 ) * ( 0.5d0 * rho * y ** 2 + + 0.5d0 * y * mu + 0.625d0 * mu ** 2 / rho ) p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 13 ) then if ( rho * y .lt. - 0.5d0 * mu ) then t = rho * y / mu v = - log ( - t ) - 1.0d0 p = mu ** 2 * tconst * v / rho else v = 2.0d0 * rho * y ** 2 + 4.0d0 * y * mu + + ( 0.5d0 + log ( 2.0d0 ) ) * mu ** 2 / rho p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 14 ) then if ( rho * y .lt. - 0.5d0 * mu ) then v = - mu ** 3 / ( rho ** 2 * y ) p = tconst * v else v = 8.0d0 * rho * y ** 2 + 12.0d0 * y * mu + + 6.0d0 * mu ** 2 / rho p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 15 ) then if ( rho * y .lt. - mu ) then t = rho * y / mu v = 4.0d0 / ( 1.0d0 - t ) - 2.0d0 p = mu ** 2 * tconst * v / rho else if ( rho * y .lt. - 0.25d0 * mu ) then t = rho * y / mu v = - log ( - t ) p = mu ** 2 * tconst * v / rho else v = 8.0d0 * y * ( rho * y + mu ) + + ( 1.5d0 + log ( 4.0d0 ) ) * mu ** 2 / rho p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 16 ) then if ( y .eq. 0.0d0 ) then v = 1.0d0 p = mu ** 2 * tconst * v / rho else v = 0.5d0 * mu * ( y + abs( y ) * + sqrt( 1.0d0 + 4.0d0 * mu ** 2 / ( rho * y ) ** 2 ) ) p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 17 ) then if ( y .le. 0.0d0 ) then t = rho * y / mu v = log ( 1.0d0 + exp ( t ) ) p = mu ** 2 * tconst * v / rho else t = rho * y / mu v = mu * y + mu ** 2 * log ( 1.0d0 + exp ( -t ) ) / rho p = tconst * v end if else if ( ptype .eq. 7 .and. ttype .eq. 18 ) then if ( y .le. 0.0d0 ) then t = rho * y / mu v = 0.5d0 * exp ( t ) p = mu ** 2 * tconst * v / rho else t = rho * y / mu v = mu * y + 0.5d0 * mu ** 2 * exp ( -t ) / rho p = tconst * v end if else if ( ptype .eq. 8 .and. ttype .eq. 5 ) then if ( mu .le. 2.0d0 ) then t = rho * y + 1.0d0 - 1.0d0 / mu if ( t .le. 0.5d0 ) then p = -( log( 1.0d0 / mu - rho * y ) + log( mu ) ) / rho else p = ( exp( 2.0d0 * ( rho * y - 1.0d0 / mu + 0.5d0 ) ) + + log( 2.0d0 ) - 1.0d0 - log( mu ) ) / rho end if else t = rho * y + 0.5d0 * ( log ( 0.5d0 * mu ) + 1.0d0 ) if ( t .le. 0.5d0 ) then p = ( 1.0d0 - log( 1.0d0 - t ) - 0.5d0 * mu - + log( 2.0d0 ) ) / rho else p = ( 0.5d0 * mu * ( exp( 2.0d0 * rho * y ) - 1.0d0 ) + + log( 2.0d0 ) + 1.0d0 ) / rho end if end if else if ( ptype .eq. 8 .and. ttype .eq. 7 ) then if ( mu .le. 4.0d0 ) then t = rho * y + 1.0d0 - 1.0d0 / sqrt ( mu ) if ( t .le. 0.5d0 ) then p = mu * y / ( 1.0d0 - sqrt( mu ) * rho * y ) else p = ( exp ( 4.0d0 * t - 2.0d0 ) - + sqrt ( mu ) + 1.0d0 ) / rho end if else t = rho * y + 0.25d0 * log ( 0.25d0 * mu ) + 0.5d0 if ( t .le. 0.5d0 ) then p = ( 1.0d0 / ( 1.0d0 - t ) - 1.0d0 - + 0.25d0 * mu ) / rho else p = ( exp ( 4.0d0 * rho * y ) - 1.0d0 ) * + 0.25d0 * mu / rho end if end if else if ( ptype .eq. 8 .and. ttype .eq. 9 ) then if ( mu .le. 0.5d0 ) then t = rho * y - 0.25d0 / mu if ( t .le. - 0.5d0 ) then p = log( - 2.0d0 * t ) + log( 2.0d0 * mu ) p = - 0.25d0 * p / rho else p = ( 0.5d0 * t ** 2 + t - 0.25d0 * + ( log( mu ) - 1.5d0 - log ( 0.5d0 ) ) ) / rho end if else t = rho * y + mu - 1.0d0 if ( t .le. - 0.5d0 ) then p = log ( - 2.0d0 * t ) + 2.0d0 * mu ** 2 - 0.5d0 p = - 0.25d0 * p / rho else p = 0.5d0 * rho * y ** 2 + y * mu end if end if else if ( ( ptype .eq. 8 .and. ttype .eq. 6 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 8 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 10 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 11 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 12 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 13 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 14 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 15 ) ) then call ytild(mu,zeta,ytildv,inform) call theta(rho * y + ytildv,thetav,inform) call theta(ytildv,thetav2,inform) p = ( thetav - thetav2 ) / rho else if ( ( ptype .eq. 9 .and. ttype .eq. 16 ) .or. + ( ptype .eq. 9 .and. ttype .eq. 17 ) .or. + ( ptype .eq. 9 .and. ttype .eq. 18 ) ) then call ytild(mu,zeta,ytildv,inform) call theta(rho * y + ytildv,thetav,inform) call theta(ytildv,thetav2,inform) p = zeta * ( thetav - thetav2 ) / rho else inform = -92 return end if c This just true if p = NaN or |p| > 1.0d+20. In this case, c project onto the box [-1.0d+20, 1.0d+20]. if ( .not. p .gt. -1.0d+20 .or. .not. p .lt. 1.0d+20 ) then p = max( -1.0d+20, min( p, 1.0d+20 ) ) end if end c **************************************************** c **************************************************** subroutine evaldpdy(y,rho,mu,zeta,dpdy,inform) C SCALAR ARGUMENTS integer inform double precision y,rho,mu,zeta,dpdy C COMMON SCALARS integer ptype,ttype double precision tconst C LOCAL SCALARS double precision t,dtdy,tmp,dtmpdy,aux,ytildv,c C COMMON BLOCKS common /p/ ptype,ttype,tconst inform = 0 if ( ptype .eq. 0 ) then if ( ttype .eq. 1 ) then c = 1.0d0 else if ( ttype .eq. 2 ) then c = 1.0d-1 else if ( ttype .eq. 3 ) then c = 1.0d-2 else if ( ttype .eq. 4 ) then c = 1.0d-3 else if ( ttype .eq. 5 ) then c = 1.0d-4 else if ( ttype .eq. 6 ) then c = 1.0d-5 else if ( ttype .eq. 7 ) then c = 1.0d-6 else if ( ttype .eq. 8 ) then c = 0.0d0 end if tmp = mu + rho * y aux = sqrt( tmp ** 2 + c / rho ) dpdy = 0.25d0 * aux * ( 1.0d0 + tmp / aux ) ** 2 else if ( ptype .eq. 1 ) then if ( y .ge. 0.0d0 ) then dpdy = mu + rho * y + 3.0d0 * rho ** 2 * y ** 2 else if ( y .ge. - mu / rho ) then dpdy = mu + rho * y else dpdy = 0.0d0 end if else if ( ptype .eq. 2 ) then if ( y .ge. 0.0d0 ) then dpdy = mu + 2.0d0 * mu * rho * y + + ( rho * y ) ** 2 / 3.0d0 else tmp = 1.0d0 - rho * y dpdy = ( mu / tmp ) * ( 1.0d0 + ( rho * y ) / tmp ) end if else if ( ptype .eq. 3 ) then if ( y .ge. 0.0d0 ) then dpdy = mu + 2.0d0 * mu * rho * y else tmp = 1.0d0 - rho * y dpdy = ( mu / tmp ) * ( 1.0d0 + ( rho * y ) / tmp ) end if else if ( ptype .eq. 4 .and. ttype .eq. 1 ) then dpdy = max( 0.0d0, mu + rho * y ) else if ( ptype .eq. 4 .and. ttype .eq. 2 ) then dpdy = max( 0.0d0, sign( 1.0d0, y ) * + abs( y ) ** ( 1.0d0 / 3.0d0) * + rho ** ( 1.0d0 / 3.0d0 ) + mu) else if ( ptype .eq. 4 .and. ttype .eq. 3 ) then dpdy = max ( 0.0d0, mu + sinh(rho * y) ) else if ((ptype .eq. 5 .or. ptype .eq. 8) .and. ttype .eq. 4) then t = rho * y dtdy = exp ( t ) dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 5 ) then t = rho * y if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) else dtdy = 2.0d0 * exp( 2.0d0 * t - 1.0d0 ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 6 ) then t = rho * y if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) else dtdy = 4.0d0 * t end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 7 ) then t = rho * y if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2 else dtdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 8 ) then t = rho * y if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2 else dtdy = 16.0d0 * t - 4.0d0 end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 9 ) then t = rho * y if ( t .le. -0.5d0 ) then dtdy = - 0.25d0 / t else dtdy = 1.0d0 + t end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 10 ) then t = rho * y aux = sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 ) tmp = ( 1.0d0 + t ) + aux dtmpdy = 1.0d0 + ( 1.0d0 + t ) / aux dtdy = 0.125d0 * tmp * dtmpdy + dtmpdy / tmp dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 11 ) then t = rho * y if ( t .le. 0.5d0 ) then dtdy = 0.5d0 * max(0.0d0, t + 0.5d0) ** 2 else dtdy = t end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 12 ) then t = rho * y if ( t .le. 0.5d0 ) then dtdy = exp ( t ) else dtdy = exp( 0.5d0 ) * ( t + 0.5d0 ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 13 ) then t = rho * y if ( t .le. -0.5d0 ) then dtdy = - 1.0d0 / t else dtdy = 4.0d0 * t + 4.0d0 end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 14 ) then t = rho * y if ( t .le. -0.5d0 ) then dtdy = 1.0d0 / t ** 2 else dtdy = 16.0d0 * t + 12.0d0 end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 15 ) then t = rho * y if ( t .le. -1.0d0 ) then dtdy = 4.0d0 / ( 1.0d0 - t ) ** 2 else if ( t .le. -0.25d0 ) then dtdy = - 1.0d0 / t else dtdy = 16.0d0 * t + 8.0d0 end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 16 ) then t = rho * y if ( abs( t ) .le. 1.0d0 ) then tmp = sqrt( t ** 2 + 4.0d0 ) else tmp = abs( t ) * sqrt( 1.0d0 + 4.0d0 / t ** 2 ) end if dtdy = 0.5d0 * ( 1.0d0 + t / tmp ) dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 17 ) then t = rho * y if ( t .le. 0.0d0 ) then dtdy = exp ( t ) / ( 1.0d0 + exp ( t ) ) else dtdy = 1.0d0 / ( 1.0d0 + exp ( -t ) ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 5 .and. ttype .eq. 18 ) then t = rho * y if ( t .le. 0.0d0 ) then dtdy = 0.5d0 * exp ( t ) else dtdy = 1.0d0 - 0.5d0 * exp ( -t ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 4 ) then t = rho * mu * y dtdy = exp ( t ) dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 5 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) else dtdy = 2.0d0 * exp( 2.0d0 * t - 1.0d0 ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 6 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) else dtdy = 4.0d0 * t end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 7 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2 else dtdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 8 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2 else dtdy = 16.0d0 * t - 4.0d0 end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 9 ) then t = rho * mu * y if ( t .le. -0.5d0 ) then dtdy = - 0.25d0 / ( rho * y ) dpdy = tconst * dtdy else dtdy = 1.0d0 + t dpdy = mu * tconst * dtdy end if else if ( ptype .eq. 6 .and. ttype .eq. 10 ) then t = rho * mu * y aux = sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 ) tmp = ( 1.0d0 + t ) + aux dtmpdy = 1.0d0 + ( 1.0d0 + t ) / aux dtdy = 0.125d0 * tmp * dtmpdy + dtmpdy / tmp dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 11 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then dtdy = 0.5d0 * max(0.0d0, t + 0.5d0) ** 2 else dtdy = t end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 12 ) then t = rho * mu * y if ( t .le. 0.5d0 ) then dtdy = exp ( t ) else dtdy = exp( 0.5d0 ) * ( t + 0.5d0 ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 13 ) then t = rho * mu * y if ( t .le. -0.5d0 ) then dtdy = - 1.0d0 / ( rho * y ) dpdy = tconst * dtdy else dtdy = 4.0d0 * t + 4.0d0 dpdy = mu * tconst * dtdy end if else if ( ptype .eq. 6 .and. ttype .eq. 14 ) then t = rho * mu * y if ( t .le. -0.5d0 ) then dtdy = 1.0d0 / ( mu * ( rho * y ) ** 2 ) dpdy = tconst * dtdy else dtdy = 16.0d0 * t + 12.0d0 dpdy = mu * tconst * dtdy end if else if ( ptype .eq. 6 .and. ttype .eq. 15 ) then t = rho * mu * y if ( t .le. -1.0d0 ) then dtdy = 4.0d0 / ( 1.0d0 - t ) ** 2 dpdy = mu * tconst * dtdy else if ( t .le. -0.25d0 ) then dtdy = - 1.0d0 / ( rho * y ) dpdy = tconst * dtdy else dtdy = 16.0d0 * t + 8.0d0 dpdy = mu * tconst * dtdy end if else if ( ptype .eq. 6 .and. ttype .eq. 16 ) then t = rho * mu * y if ( abs( t ) .le. 1.0d0 ) then tmp = sqrt( t ** 2 + 4.0d0 ) else tmp = abs( t ) * sqrt( 1.0d0 + 4.0d0 / t ** 2 ) end if dtdy = 0.5d0 * ( 1.0d0 + t / tmp ) dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 17 ) then t = rho * mu * y if ( t .le. 0.0d0 ) then dtdy = exp ( t ) / ( 1.0d0 + exp ( t ) ) else dtdy = 1.0d0 / ( 1.0d0 + exp ( -t ) ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 6 .and. ttype .eq. 18 ) then t = rho * mu * y if ( t .le. 0.0d0 ) then dtdy = 0.5d0 * exp ( t ) else dtdy = 1.0d0 - 0.5d0 * exp ( -t ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 7 .and. ttype .eq. 4 ) then t = rho * y / mu dtdy = exp ( t ) dpdy = mu * tconst * dtdy else if ( ptype .eq. 7 .and. ttype .eq. 5 ) then t = rho * y / mu if ( rho * y .le. 0.5d0 * mu ) then dtdy = 1.0d0 / ( 1.0d0 - t ) else dtdy = 2.0d0 * exp( 2.0d0 * t - 1.0d0 ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 7 .and. ttype .eq. 6 ) then if ( rho * y .lt. 0.5d0 * mu ) then t = rho * y / mu dtdy = 1.0d0 / ( 1.0d0 - t ) dpdy = mu * tconst * dtdy else dtdy = 4.0d0 * rho * y dpdy = tconst * dtdy end if else if ( ptype .eq. 7 .and. ttype .eq. 7 ) then t = rho * y / mu if ( rho * y .le. 0.5d0 * mu ) then dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2 else dtdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 7 .and. ttype .eq. 8 ) then if ( rho * y .lt. 0.5d0 * mu ) then t = rho * y / mu dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2 dpdy = mu * tconst * dtdy else dtdy = 16.0d0 * rho * y - 4.0d0 * mu dpdy = tconst * dtdy end if else if ( ptype .eq. 7 .and. ttype .eq. 9 ) then if ( rho * y .lt. - 0.5d0 * mu ) then dtdy = - 0.25d0 * mu ** 2 / ( rho * y ) dpdy = tconst * dtdy else dtdy = mu + rho * y dpdy = tconst * dtdy end if else if ( ptype .eq. 7 .and. ttype .eq. 10 ) then t = rho * y / mu aux = sqrt( ( 1.0d0 + t ) ** 2 + 8.0d0 ) tmp = ( 1.0d0 + t ) + aux dtmpdy = 1.0d0 + ( 1.0d0 + t ) / aux dtdy = 0.125d0 * tmp * dtmpdy + dtmpdy / tmp dpdy = mu * tconst * dtdy else if ( ptype .eq. 7 .and. ttype .eq. 11 ) then if ( rho * y .lt. 0.5d0 * mu ) then t = rho * y / mu dtdy = 0.5d0 * max(0.0d0, t + 0.5d0) ** 2 dpdy = mu * tconst * dtdy else dtdy = rho * y dpdy = tconst * dtdy end if else if ( ptype .eq. 7 .and. ttype .eq. 12 ) then if ( rho * y .lt. 0.5d0 * mu ) then t = rho * y / mu dtdy = exp ( t ) dpdy = mu * tconst * dtdy else dtdy = exp( 0.5d0 ) * ( rho * y + 0.5d0 * mu ) dpdy = tconst * dtdy end if else if ( ptype .eq. 7 .and. ttype .eq. 13 ) then if ( rho * y .lt. - 0.5d0 * mu ) then dtdy = - mu ** 2 / ( rho * y ) dpdy = tconst * dtdy else dtdy = 4.0d0 * rho * y + 4.0d0 * mu dpdy = tconst * dtdy end if else if ( ptype .eq. 7 .and. ttype .eq. 14 ) then if ( rho * y .lt. - 0.5d0 * mu ) then dtdy = mu ** 3 / ( rho * y ) ** 2 dpdy = tconst * dtdy else dtdy = 16.0d0 * rho * y + 12.0d0 * mu dpdy = tconst * dtdy end if else if ( ptype .eq. 7 .and. ttype .eq. 15 ) then if ( rho * y .lt. - mu ) then t = rho * y / mu dtdy = 4.0d0 / ( 1.0d0 - t ) ** 2 dpdy = mu * tconst * dtdy else if ( rho * y .lt. - 0.25d0 * mu ) then dtdy = - mu ** 2 / ( rho * y ) dpdy = tconst * dtdy else dtdy = 16.0d0 * rho * y + 8.0d0 * mu dpdy = tconst * dtdy end if else if ( ptype .eq. 7 .and. ttype .eq. 16 ) then dtdy = 0.5d0 * ( 1.0d0 + sign( 1.0d0, y ) / + sqrt( 1.0d0 + 4.0d0 * mu ** 2 / ( rho * y ) ** 2 ) ) dpdy = mu * tconst * dtdy else if ( ptype .eq. 7 .and. ttype .eq. 17 ) then t = rho * y / mu if ( y .le. 0.0d0 ) then dtdy = exp ( t ) / ( 1.0d0 + exp ( t ) ) else dtdy = 1.0d0 / ( 1.0d0 + exp ( -t ) ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 7 .and. ttype .eq. 18 ) then t = rho * y / mu if ( y .le. 0.0d0 ) then dtdy = 0.5d0 * exp ( t ) else dtdy = 1.0d0 - 0.5d0 * exp ( -t ) end if dpdy = mu * tconst * dtdy else if ( ptype .eq. 8 .and. ttype .eq. 5 ) then if ( mu .le. 2.0d0 ) then t = 1.0d0 + rho * y - 1.0d0 / mu if ( t .le. 0.5d0 ) then dpdy = 1.0d0 / ( 1.0d0 / mu - rho * y ) else dpdy = 2.0d0 * + exp( 2.0d0 * ( rho * y - 1.0d0 / mu + 0.5d0 ) ) end if else t = rho * y + 0.5d0 * ( log ( 0.5d0 * mu ) + 1.0d0 ) if ( t .le. 0.5d0 ) then dpdy = 1.0d0 / ( 1.0d0 - t ) else dpdy = mu * exp( 2.0d0 * rho * y ) end if end if else if ( ptype .eq. 8 .and. ttype .eq. 7 ) then if ( mu .le. 4.0d0 ) then t = rho * y + 1.0d0 - 1.0d0 / sqrt ( mu ) if ( t .le. 0.5d0 ) then dpdy = 1.0d0 / ( 1.0d0 / sqrt( mu ) - rho * y ) ** 2 else dpdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 ) end if else t = rho * y + 0.25d0 * log ( 0.25d0 * mu ) + 0.5d0 if ( t .le. 0.5d0 ) then dpdy = 1.0d0 / ( 1.0d0 - t ) ** 2 else dpdy = exp ( 4.0d0 * rho * y ) * mu end if end if else if ( ptype .eq. 8 .and. ttype .eq. 9 ) then if ( mu .le. 0.5d0 ) then t = rho * y - 0.25d0 / mu if ( t .le. - 0.5d0 ) then dpdy = - 0.25d0 / t else dpdy = t + 1.0d0 end if else tmp = rho * y + mu if ( tmp .le. 0.5d0 ) then dpdy = - 0.125d0 / ( tmp - 1.0d0 ) else dpdy = tmp end if end if else if ( ( ptype .eq. 8 .and. ttype .eq. 6 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 8 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 10 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 11 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 12 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 13 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 14 ) .or. + ( ptype .eq. 8 .and. ttype .eq. 15 ) ) then call ytild(mu,zeta,ytildv,inform) call evaldtdy(rho * y + ytildv,dtdy,inform) dpdy = dtdy else if ( ( ptype .eq. 9 .and. ttype .eq. 16 ) .or. + ( ptype .eq. 9 .and. ttype .eq. 17 ) .or. + ( ptype .eq. 9 .and. ttype .eq. 18 ) ) then call ytild(mu,zeta,ytildv,inform) call evaldtdy(rho * y + ytildv,dtdy,inform) dpdy = zeta * dtdy else inform = -93 return end if c This just true if p = NaN or |p| > 1.0d+20. In this case, c project onto the box [-1.0d+20, 1.0d+20]. if ( .not. dpdy .gt. -1.0d+20 .or. .not. dpdy .lt. 1.0d+20 ) +then dpdy = max( -1.0d+20, min( dpdy, 1.0d+20 ) ) end if end c **************************************************** c **************************************************** subroutine theta(t,thetav,inform) C SCALAR ARGUMENTS integer inform double precision t,thetav C COMMON SCALARS integer ptype,ttype double precision tconst C LOCAL SCALARS double precision tpl,tmp C COMMON BLOCKS common /p/ ptype,ttype,tconst inform = 0 if ( ttype .eq. 1 ) then thetav = 0.5d0 * t**2 else if ( ttype .eq. 2 ) then thetav = 0.75d0 * abs(t) ** (4.0d0/3.0d0) else if ( ttype .eq. 3 ) then thetav = cosh(t) - 1.0d0 else if ( ttype .eq. 4 ) then thetav = exp(t) - 1.0d0 else if ( ttype .eq. 5 ) then if ( t .le. 0.5d0 ) then thetav = -log( 1.0d0 - t ) else thetav = exp( 2.0d0 * t - 1.0d0 ) + log(2.0d0) - 1.0d0 end if else if ( ttype .eq. 6 ) then if ( t .le. 0.5d0 ) then thetav = -log( 1.0d0 - t ) else thetav = 2.0d0 * t ** 2 + log(2.0d0) - 0.5d0 end if else if ( ttype .eq. 7 ) then if ( t .le. 0.5d0 ) then thetav = 1.0d0 / ( 1.0d0 - t ) - 1.0d0 else thetav = exp ( 4.0d0 * t - 2.0d0 ) end if else if ( ttype .eq. 8 ) then if ( t .le. 0.5d0 ) then thetav = 1.0d0 / ( 1.0d0 - t ) - 1.0d0 else thetav = 8.0d0 * t ** 2 - 4.0d0 * t + 1.0d0 end if else if ( ttype .eq. 9 ) then if ( t .le. -0.5d0 ) then thetav = -0.25d0 * log ( -2.0d0 * t ) - 0.375d0 else thetav = 0.5d0 * t ** 2 + t end if else if ( ttype .eq. 10 ) then tpl = 1.0d0 + t if ( abs( tpl ) .le. 1.0d0 ) then tmp = tpl + sqrt( tpl ** 2 + 8.0d0 ) else tmp = tpl + abs( tpl ) * sqrt( 1.0d0 + 8.0d0 / tpl ** 2 ) end if thetav = 0.0625d0 * tmp ** 2 + log ( 0.25d0 * tmp ) - 1.0d0 else if ( ttype .eq. 11 ) then if ( t .le. 0.5d0 ) then thetav = max(0.0d0, t + 0.5d0)**3 / 6.0d0 - 1.0d0 / 24.0d0 else thetav = 0.5d0 * t ** 2 end if else if ( ttype .eq. 12 ) then if ( t .le. 0.5d0 ) then thetav = exp (t) else thetav = exp(0.5d0) * (0.5d0*t**2 + 0.5d0*t + 0.625d0) end if else if ( ttype .eq. 13 ) then if ( t .le. -0.5d0 ) then thetav = - log ( -t ) - 1.0d0 else thetav = 2.0d0 * t**2 + 4.0d0 * t + 0.5d0 + log ( 2.0d0 ) end if else if ( ttype .eq. 14 ) then if ( t .le. -0.5d0 ) then thetav = - 1.0d0 / t else thetav = 8.0d0 * t**2 + 12.0d0 * t + 6.0d0 end if else if ( ttype .eq. 15 ) then if ( t .le. -1.0d0 ) then thetav = 4.0d0 / ( 1.0d0 - t ) - 2.0d0 else if ( t .le. -0.25d0 ) then thetav = - log ( -t ) else thetav = 8.0d0 * t**2 + 8.0d0 * t + 1.5d0 + log ( 4.0d0 ) end if else if ( ttype .eq. 16 ) then if ( abs( t ) .le. 1.0d0 ) then tmp = sqrt( t ** 2 + 4.0d0 ) else tmp = abs( t ) * sqrt( 1.0d0 + 4.0d0 / t ** 2 ) end if thetav = 0.5d0 * ( t + tmp ) else if ( ttype .eq. 17 ) then if ( t .le. 0.0d0 ) then thetav = log ( 1.0d0 + exp ( t ) ) else thetav = t + log ( 1.0d0 + exp ( -t ) ) end if else if ( ttype .eq. 18 ) then if ( t .le. 0.0d0 ) then thetav = 0.5d0 * exp ( t ) else thetav = t + 0.5d0 * exp ( -t ) end if end if thetav = thetav * tconst end c **************************************************** c **************************************************** subroutine evaldtdy(t,dtdy,inform) C SCALAR ARGUMENTS integer inform double precision t,dtdy C COMMON SCALARS integer ptype,ttype double precision tconst C LOCAL SCALARS double precision tpl,tmp,dtmpdc,aux C COMMON BLOCKS common /p/ ptype,ttype,tconst inform = 0 if ( ttype .eq. 1 ) then dtdy = t else if ( ttype .eq. 2 ) then dtdy = sign(1.0d0,t) * abs(t) ** (1.0d0/3.0d0) else if ( ttype .eq. 3 ) then dtdy = sinh ( t ) else if ( ttype .eq. 4 ) then dtdy = exp ( t ) else if ( ttype .eq. 5 ) then if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) else dtdy = 2.0d0 * exp( 2.0d0 * t - 1.0d0 ) end if else if ( ttype .eq. 6 ) then if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) else dtdy = 4.0d0 * t end if else if ( ttype .eq. 7 ) then if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2 else dtdy = 4.0d0 * exp ( 4.0d0 * t - 2.0d0 ) end if else if ( ttype .eq. 8 ) then if ( t .le. 0.5d0 ) then dtdy = 1.0d0 / ( 1.0d0 - t ) ** 2 else dtdy = 16.0d0 * t - 4.0d0 end if else if ( ttype .eq. 9 ) then if ( t .le. -0.5d0 ) then dtdy = - 0.25d0 / t else dtdy = 1.0d0 + t end if else if ( ttype .eq. 10 ) then tpl = 1.0d0 + t aux = sqrt( 1.0d0 + 8.0d0 / tpl ** 2 ) tmp = tpl + abs( tpl ) * aux dtmpdc = 1.0d0 + sign(1.0d0,tpl) / aux dtdy = 0.125d0 * tmp * dtmpdc + dtmpdc / tmp else if ( ttype .eq. 11 ) then if ( t .le. 0.5d0 ) then dtdy = 0.5d0 * max(0.0d0, t + 0.5d0) ** 2 else dtdy = t end if else if ( ttype .eq. 12 ) then if ( t .le. 0.5d0 ) then dtdy = exp ( t ) else dtdy = exp( 0.5d0 ) * ( t + 0.5d0 ) end if else if ( ttype .eq. 13 ) then if ( t .le. -0.5d0 ) then dtdy = - 1.0d0 / t else dtdy = 4.0d0 * t + 4.0d0 end if else if ( ttype .eq. 14 ) then if ( t .le. -0.5d0 ) then dtdy = 1.0d0 / t ** 2 else dtdy = 16.0d0 * t + 12.0d0 end if else if ( ttype .eq. 15 ) then if ( t .le. -1.0d0 ) then dtdy = 4.0d0 / ( 1.0d0 - t ) ** 2 else if ( t .le. -0.25d0 ) then dtdy = - 1.0d0 / t else dtdy = 16.0d0 * t + 8.0d0 end if else if ( ttype .eq. 16 ) then if ( abs( t ) .le. 1.0d0 ) then tmp = sqrt( t ** 2 + 4.0d0 ) else tmp = abs( t ) * sqrt( 1.0d0 + 4.0d0 / t ** 2 ) end if dtdy = 0.5d0 * ( 1.0d0 + t / tmp ) else if ( ttype .eq. 17 ) then if ( t .le. 0.0d0 ) then dtdy = exp ( t ) / ( 1.0d0 + exp ( t ) ) else dtdy = 1.0d0 / ( 1.0d0 + exp ( -t ) ) end if else if ( ttype .eq. 18 ) then if ( t .le. 0.0d0 ) then dtdy = 0.5d0 * exp ( t ) else dtdy = 1.0d0 - 0.5d0 * exp ( -t ) end if end if dtdy = dtdy * tconst end c **************************************************** c **************************************************** subroutine ytild(mu,zeta,ytildv,inform) C SCALAR ARGUMENTS integer inform double precision mu,zeta,ytildv C COMMON SCALARS integer ptype,ttype double precision tconst C COMMON BLOCKS common /p/ ptype,ttype,tconst inform = 0 if ( ttype .eq. 4 ) then ytildv = log ( mu ) else if ( ttype .eq. 5 ) then if ( mu .le. 2.0d0 ) then ytildv = 1.0d0 - 1.0d0 / mu else ytildv = 0.5d0 * ( log ( 0.5d0 * mu ) + 1.0d0 ) end if else if ( ttype .eq. 6 ) then if ( mu .le. 2.0d0 ) then ytildv = 1.0d0 - 1.0d0 / mu else ytildv = 0.25d0 * mu end if else if ( ttype .eq. 7 ) then if ( mu .le. 4.0d0 ) then ytildv = 1.0d0 - 1.0d0 / sqrt ( mu ) else ytildv = 0.25d0 * log ( 0.25d0 * mu ) + 0.5d0 end if else if ( ttype .eq. 8 ) then if ( mu .le. 4.0d0 ) then ytildv = 1.0d0 - 1.0d0 / sqrt ( mu ) else ytildv = 0.0625d0 * mu + 0.25d0 end if else if ( ttype .eq. 9 ) then if ( mu .le. 0.5d0 ) then ytildv = - 0.25d0 / mu else ytildv = mu - 1.0d0 end if else if ( ttype .eq. 10 ) then ytildv = 2.0d0 * mu - 1.0d0 / mu - 1.0d0 else if ( ttype .eq. 11 ) then if ( mu .le. 0.5d0 ) then ytildv = sqrt ( 2.0d0 * mu ) - 0.5d0 else ytildv = mu end if else if ( ttype .eq. 12 ) then if ( mu .le. exp (0.5d0) ) then ytildv = log ( mu ) else ytildv = mu * exp( -0.5d0 ) - 0.5d0 end if else if ( ttype .eq. 13 ) then if ( mu .le. 2.0d0 ) then ytildv = - 1.0d0 / mu else ytildv = 0.25d0 * mu - 1.0d0 end if else if ( ttype .eq. 14 ) then if ( mu .le. 4.0d0 ) then ytildv = - 1.0d0 / sqrt ( mu ) else ytildv = 0.0625d0 * mu - 0.75d0 end if else if ( ttype .eq. 15 ) then if ( mu .le. 1.0d0 ) then ytildv = 1.0d0 - 2.0d0 / sqrt ( mu ) else if ( mu .le. 4.0d0 ) then ytildv = - 1.0d0 / mu else ytildv = 0.0625d0 * mu - 0.5d0 end if else if ( ttype .eq. 16 ) then ytildv = ( 2.0d0 * mu - zeta ) / sqrt( mu * zeta - mu ** 2 ) else if ( ttype .eq. 17 ) then ytildv = log ( mu / ( zeta - mu ) ) else if ( ttype .eq. 18 ) then if ( mu .le. 0.5d0 * zeta ) then ytildv = log ( 2.0d0 * mu / zeta ) else ytildv = - log ( ( 2.0d0 * ( zeta - mu ) ) / zeta ) end if end if end c **************************************************** c **************************************************** subroutine gotc(n,x,m,c,got) C SCALAR ARGUMENTS logical got integer n,m C ARRAY ARGUMENTS double precision c(n),x(n) C PARAMETERS integer nmax,mmax parameter (nmax=100000,mmax=100000) C LOCAL SCALARS integer k,i logical samep C COMMON SCALAR integer nc C COMMON ARRAYS double precision xs(nmax,2),cs(mmax,2) C COMMON BLOCKS common /constr/xs,cs,nc got= .false. k= 1 100 if (k.le.min0(2,nc)) then i= 1 samep= .true. 110 if (i.le.n.and.samep) then if (x(i).ne.xs(i,k)) then samep= .false. end if i= i + 1 go to 110 end if if (samep) then got= .true. do i= 1, m c(i)= cs(i,k) end do return end if k= k + 1 go to 100 end if return end c **************************************************** c **************************************************** subroutine savec(n,x,m,c) C SCALAR ARGUMENTS integer n,m C ARRAY ARGUMENTS double precision c(n),x(n) C PARAMETERS integer nmax,mmax parameter (nmax=100000,mmax=100000) C LOCAL SCALARS integer i C COMMON SCALAR integer nc C COMMON ARRAYS double precision xs(nmax,2),cs(mmax,2) C COMMON BLOCKS common /constr/xs,cs,nc if (nc.lt.2) then nc= nc + 1 end if do i= 1, n xs(i,2)= xs(i,1) end do do i= 1, m cs(i,2)= cs(i,1) end do do i= 1, n xs(i,1)= x(i) end do do i= 1, m cs(i,1)= c(i) end do return end c **************************************************** c **************************************************** double precision function norm2s(n,x) C SCALAR ARGUMENTS integer n C ARRAY ARGUMENTS 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 May 14th, 2003. c c 1) The way amax2 and amax2n were being computing may casued a c segmentation fault. Its initialization was changed from infty and c -infty to 1.0d+99 and -1.0d+99, respectively. Using infty, when c combined with a big trust region radius, the final value of amax2 c or amax2n may cause the impression that a bound is being attained, c when it is not. "Redundant" ifs inside the amax2 and anax2n calculation c were deleted. It should considered the possibility of using two c conatants, namely, bignum = 1.0d+20 and infty = 1.0d+99, instead c of just infty. c Modifications introduced near April 1st 2003 in ocassion of the c PHR and inequality-constraints Augemented Lagrangean methods c development: c c 1) The use of iprint was redefined and iprint2 was deleted. c c 2) The way to detect no progress in the log of the projected c gradient norm was changed. As it was, "no progess means no c reduction in the projected gradient norm over M iterations". c But this criterion implicitely assumed that the projected c gradient norm msut decrease monotonicaly. Is it is clearly not c true, the criterion was changed by a non-monotone decrease c criterion. Now, progress means that the projected gradient c norm is, at each iteration, smaller than the maximum over the c last M iterations. And "no porgess" means the it does not occurs c during not smaller than the c c 3 ) The computation of qamaxn inside cg subroutine was in the wrong c place (it was being used before computed) and it was maybe the c reason for which the option nerlyq = .true. never worked properly. c With this correction this option should be tested again. 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,x,d,g,n,xc,gtype,htvtype,hd,xtemp, *sterel,steabs,goth,inform) C SCALAR ARGUMENTS logical goth integer nind,n,gtype,htvtype,inform double precision sterel,steabs C ARRAY ARGUMENTS integer ind(nind) double precision x(n),d(n),g(n),xc(n),hd(n),xtemp(n) C PARAMETERS integer mmax,jcnzemax parameter (mmax=100000,jcnzemax=100000) C COMMON SCALARS integer m C COMMON ARRAYS double precision mubar(mmax),zeta(mmax),rho(mmax) C LOCAL SCALARS logical got integer i,jcnze,flag double precision xinfn,dinfn,step,f C LOCAL ARRAYS integer jcvar(jcnzemax),jcfun(jcnzemax) double precision c(mmax),ctemp(mmax),dpdc(mmax),dpdctemp(mmax), + jcval(jcnzemax) C COMMON BLOCKS common /auglag/mubar,zeta,rho,m inform = 0 c Compute incremental quotients step xinfn = 0.0d0 dinfn = 0.0d0 do i = 1,n xinfn = max( xinfn, abs( x(i) ) ) dinfn = max( dinfn, abs( d(i) ) ) end do step = max( sterel * xinfn, steabs ) / dinfn c Set the point in which the gradient will be c evaluated: xtemp = x + step * d do i = 1,n xtemp(i) = x(i) + step * d(i) end do c Evaluate the gradient at (x + step * d) C COMPUTE CONSTRAINTS AT x AND xtemp call gotc(n,x,m,c,got) if ( .not. got ) then call evalfc(n,x,f,m,c,flag) if ( flag .ne. 0 ) then inform = -90 return end if end if call evalfc(n,xtemp,f,m,ctemp,flag) if ( flag .ne. 0 ) then inform = -90 return end if C COMPUTE dP/dc do i = 1,m dpdc(i) = mubar(i) + rho(i) * c(i) dpdctemp(i) = mubar(i) + rho(i) * ctemp(i) end do C COMPUTE \nabla L(x) = \nabla f(x) + \sum_{i=1}^m dPdc * dcdx do i = 1,n hd(i) = 0.0d0 end do call evalgjc(n,xtemp,m,jcvar,jcfun,jcval,jcnze,flag) if ( flag .ne. 0 ) then inform = -91 return end if do i = 1,jcnze if ( jcfun(i) .eq. 0 ) then hd(jcvar(i)) = hd(jcvar(i)) + jcval(i) else if ( dpdc(jcfun(i)) .gt. 0.0d0 ) then hd(jcvar(i)) = + hd(jcvar(i)) + dpdctemp(jcfun(i)) * jcval(i) end if end if end do c Compute incremental quotients do i = 1,n hd(i) = ( hd(i) - g(i) ) / step end do return end