C IFFCO -- A bound constrained optimization software package. C Author: Paul Gilmore, C North Carolina Sate University C Department of Mathematics C Box 8205 C Raleigh, N.C. 27695-8205 C E-mail: paul@latour.math.ncsu.edu C Date Written 02/11/93 MM/DD/YY C C 1997 update by Tony Choi C C the user-supplied function format matches the manual. C C func is initialized to zero. C C equality is tested for between 2 double precision #'s C x,y by testing |x - y| < eps where eps = 10^{-16}. C C A function evaluation is no longer done every time the scale C is reduced. The current point has not moved so there is no need. C C The line search no longer evaluates the function if it has already C evaluated it at the same point. C C subroutine iffco(fcn,x,u,l,fscale,rminh,rmaxh,n, + maxit,restart,writ,tol,f,ncuts) integer maxor C C Set the maximum order for the arrays used. C parameter(maxor=100) integer n,maxit,restart,writ,ncuts,oops integer func,linf,cutnum,i,k,minscal,jpvt(maxor) real*8 fcn,x(n),u(n),l(n),fscale,rminh,rmaxh,tol,f external fcn real*8 y(maxor),fp,gf2,xs1(maxor),xs2(maxor),p(maxor) real*8 ifc2nor,s(maxor),engf,sr1(maxor,maxor),lp(maxor) real*8 ox(maxor),ogf(maxor),enp,v(maxor),h real*8 b(maxor,maxor),t(maxor),work(maxor),gf(maxor) real*8 meps CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFFCO C C Subroutine iffco is a projected quasi-Newton algorithm for constrained C nonlinear minimization. Iffco is designed to solve problems of the C form: C C min f: Q --> R, C C where f is the function to be minimized C and Q is an n-dimensional hyper-box given by the the following C equation: C C Q={ x : l(i) <= x(i) <= u(i), i = 1,...,n }. C C A brief outline of the algorithm for iffco follows. C C The main input variables used by iffco are: C 1) An initial point in the hyper-box Q. C 2) An initial finite difference step (scale) used to calculate C the gradient. C 3) A lower bound for the scales used to calculate the gradient. C C Iffco calculates the function and the gradient of the C function at the given point. It then tests for convergence C at the current point with the current scale. If the convergence criteria C is met, the scale is reduced by a half and the convergence criteria is again C checked. This process continues until either the convergence criteria C is not satisfied or until the scale is reduced to below the lower bound for C scales. If the scale has been reduced below the lower bound, the algorithm C terminates. Otherwise, iffco calculates a descent direction (see the C subroutine ifcstep for comments on how this descent direction is chosen) C and then attempts to find an acceptable new iterate in that direction C using a cutback line search algorithm (see subroutine ifclines for details). C If a new point cannot be found, the scale is reduced by half and the C process continues. If a new point is found, iffco continues the process C of checking for termination and looking for a new acceptable point at the C current point, scale pair. C C Iffco is designed to solve problems where the function f has the form C C f(x) = g(x) + o(x). C C The function g(x) is a smooth function with a simple form, such as C a convex quadratic. The function o(x) is a low amplitude high frequency C function. Iffco is especially effective on problems where the amplitude C of o(x) decays near minima of g(x). However iffco has proven to be C effective on more complex problems then those described above. C C On entry C C fcn -- The argument containing the name of the user-supplied C subroutine that returns values for the function to be C minimized. C C x -- A double-precision vector of length n. The initial point C in the iterative process. X must be within the hyper-box C defined by the constraints. C C u -- A double-precision vector of length n. The vector C containing the upper bounds for the n independent C variables. C C l -- A double-precision vector of length n. The vector C containing the lower bounds for the n independent C variables. C C fscale -- A user-supplied double-precision constant used to scale C the function. Fscale should be an approximation to the C maximum size of the absolute value of the function f in C the hyper-box. The default value for fscale is 1.d0. C C rminh -- A user-supplied double-precision constant. Rminh is a C lower bound for the last and smallest scale the algorithm C uses to calculate the finite difference gradient. C As soon as the scales are reduced to be less than rminh, C the algorithm either terminates or restarts. However, C an optimization run is always done with the scale C rmaxh regardless of the value of rminh. Rminh must C satisfy: C C 0 < rminh. C C If the above inequality is not satisfied rminh is set C to 1.d-4. See the users guid for help in determing rminh. C C rmaxh -- A user-supplied double-precision constant. Rmaxh is the C first and largest scale the algorithm uses to calculate C the finite difference gradient. Rmaxh must satisfy: C C 0 < rmaxh < 0.5. C C If the above inequality is not satisfied rminh is set C to 5.d-1. See the users guid for help in determing rmaxh. C C n -- An integer. The dimension of the problem. If n > 100 the C user must set the parameter maxor to be greater than n. C C maxit -- A user-supplied integer. Maxit specifies the maximum number C of iterations the algorithm is allowed to take at a given C scale. The default value for maxit is 200. C C restart -- A user-supplied integer. If restart is set to 1, C the algorithm restarts after the smallest scale C has been used. The restart uses the last point C generated by the algorithm as the new initial point C and rmaxh as the first scale in the sequence of finite C difference steps. C C writ -- A user-supplied integer. Writ specifies the verbosity C of the algorithm: C C writ = 1 the algorithm writes error C messages and the values of C user-supplied parameters C to standard output. C C writ = 2 the algorithm produces the C same information as for write=1 C plus a history of the iterative C process. C C If writ is equal to any other number no output is produced. C C tol -- A user-supplied double-precision constant used for determining C convergence of the algorithm at a given point for a given scale. C The determination of appropriate values of tol is C problem-dependent. See the user's guide for help determining C tol. The default value for tol is 1. C C ncuts -- A user-supplied integer. Ncuts specifies the maximum number C cutbacks in the line search the algorithm may take. C The default value for ncuts is 10. C C On return C C x -- The final point obtained in the optimization process. C X should be a good approximation to the global minimum C for the function within the hyper-box. C C f -- The value of the function at x. C C Subroutines and Functions C C ifcpreprc ifcinfcn, ifcgrad, ifcstep, ifclines, ifcprint, ifc2nor, ifcpstprc C C The user must also have access to LINPACK, because the two LINPACK C subroutines dchdc and dposl are called in the subroutine ifcstep. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Initialize func & machine epsilon func = 0 meps = 1.d-16 C C Scale the initial iterate x for use by iffco. C Check for errors in the user supplied input variables. C call ifcpreprc(x,u,l,n,fscale,tol,maxit,rmaxh,rminh +,ncuts,writ,restart,xs1,xs2,oops) if(oops.eq.1) return C C If the verbose printing option (writ = 2) has been selected, C write the following header to standard output. C if(writ.eq.2)then write(6,150) 150 format(3x,'K ',3x,'| X |',8x,'F', + 9x,'| GF |',8x,'H',9x,'CUTS',/) end if C C Set the scale to the largest allowable scale, rmaxh. C h=rmaxh C C Set the check for a minimum at all scales to on. C 160 minscal=1 C C Calculate the function at the current point with the new scale. C call ifcinfcn(fcn,x,xs1,xs2,fscale,n,func,f) C C Set the iteration counter to zero. C 170 k=0 C C Calculate the gradient of the function at the current point with the new C scale. If the ith constraint is active and the ith component of the C gradient gives a direction out of the hyper-box, set the ith component of C the gradient to zero. C call ifcgrad(fcn,x,xs1,xs2,fscale,f,n,h,func,gf) do 190 i=1,n if(dabs(x(i)).le.meps.and.gf(i).gt.0.d0) then gf(i)=0.d0 end if if(dabs(x(i)-1.d0).le.meps.and.gf(i).lt.0.d0) then gf(i)=0.d0 end if 190 continue C C Begin the iterative process, using the current scale. C do 500 k=1,maxit C C Calculate the projection of the steepest descent step, x-gf, C onto the hyper-box. C do 200 i=1,n if(x(i)-gf(i).ge.1.d0) then y(i)=1.d0 else if(x(i)-gf(i).le.0.d0) then y(i)=0.d0 else y(i)=x(i)-gf(i) end if end if 200 continue C C Calculate the projected gradient. This is done by subtracting the projected C steepest descent step from the current point in the iterative process. C In iffco, the projected gradient takes the place of the gradient for C the purpose of checking for convergence. C do 210 i=1,n 210 s(i)=x(i)-y(i) C C Calculate the two norm of the projected gradient. C engf=ifc2nor(s,n) C C Test for convergence. This is done by checking if the two norm of the projected C gradient is less than the product of the tolerance tol and the current scale. C if(engf.lt.tol*h) then if(writ.eq.2) then call ifcprint(x,f,gf,n,k,h,0,1) end if goto 1000 end if C C Calculate a negative descent direction p. The length of this vector is C called the step length. C call ifcstep(x,ox,gf,ogf,sr1,s,v,t,b,work,jpvt,n,k-1,p) C C If the step length is less than the length of the current scale, C force the negative descent direction to have length equal to the current C scale. C enp=ifc2nor(p,n) if(enp.lt.h) then do 249 i=1,n 249 p(i)=(h/enp)*p(i) end if C C Calculate the inner product of the gradient and the negative descent direction. C gf2=0.d0 do 250 i=1,n 250 gf2=gf2+gf(i)*p(i) C C Perform a backtracking line search using the descent direction -p. C call ifclines(fcn,f,x,gf2,p,xs1,xs2,fscale,h,ncuts,n,func,lp +,linf,y,fp,cutnum) C C If the line search was successful, save the old iterate x, the old function C value f, and the old gradient gf. Update x and the function value at x. C Update the counter minscal to indicate that a step has been taken. C if(linf.eq.1)then if(writ.eq.2)then call ifcprint(x,f,gf,n,k,h,cutnum,3) end if do 330 i=1,n ogf(i)=gf(i) ox(i)=x(i) 330 x(i)=y(i) f=fp minscal=0 goto 370 end if C C If the line search failed, reduce the scale and check for termination C of the algorithm. C if(writ.eq.2)then call ifcprint(x,f,gf,n,k,h,0,0) end if goto 1000 C C Calculate the gradient of f at the new point and zero out the appropriate C elements as described above. C 370 call ifcgrad(fcn,x,xs1,xs2,fscale,f,n,h,func,gf) do 390 i=1,n if(dabs(x(i)).le.meps.and.gf(i).gt.0.d0)then gf(i)=0.d0 end if if(dabs(x(i)-1.d0).le.meps.and.gf(i).lt.0.d0)then gf(i)=0.d0 end if 390 continue 500 continue C C If the maximum number of iterations has been exceeded, reduce the scale C and check for termination of the algorithm. C if(writ.eq.2)then write(6,*)'Maxit exceeded.' end if 1000 continue C C Reduce the scale. C h=(5.d-1)*h C C Check if the scale has been reduced to the lower bound rminh. C if(h.lt.rminh) then C C If no restarts are to be done, then terminate. C if(restart.le.0) goto 2000 C If no steps have been taken at any of the scales, then terminate. C if(minscal.eq.1)then goto 2000 C C Set the scale to the maximum value rmaxh, and restart the algorithm at the C current point. C else h=rmaxh if(writ.eq.2)then write(6,*)'restart' end if end if end if C C If the scale is not less then rminh, continue the iterative process with C the new scale. C goto 170 2000 continue C C Unscale the final soution x and the value of the function at x. C call ifcpstprc(x,xs1,xs2,f,fscale,n,func,writ) return end subroutine ifcpreprc(x,u,l,n,fscale,tol,maxit,rmaxh,rminh +,ncuts,writ,restart,xs1,xs2,oops) integer n,i,maxit,ncuts,writ,restart,oops real*8 x(n),u(n),l(n),xs1(n),xs2(n),rminh,rmaxh real*8 fscale,tol,meps CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCPREPRC C C Subroutine ifcpreprc is the preprocessing subroutine for iffco. C Subroutine ifcpreprc uses an afine mapping to map the hyper-box given by C the constraints on the variable x onto the n-dimensional unit cube. C This mapping is done using the following equation: C C x(i)=x(i)/(u(i)-l(i))-l(i)/(u(i)-l(i)). C C Ifcpreprc checks the following: C C 1. if the bounds l and u are well-defined. That is, if C C l(i) < u(i) forevery i. C C 2. if x is in the hyper-box defined by l and u. That is, if C C l(i) <= x(i) <= u(i). C C 3. if the termination constant tol is set. If not, tol is set C to 1. C C 4. if the maximum scale rmaxh has an acceptable value. That is, C if 0 < rmaxh <= 0.5. If not, rmaxh is set to 0.5. C C 5. if the minimum scale rminh has an acceptable value. That is, C if 0 < rminh. If not, rminh is set to 1.d-4. C C 6. if the maximum number of iterates for each scale is set. C If not, maxit is set to 200. C C 7. if the constant for scaling the function is set to an C acceptable value. That is, if fscale > 0. If not, fscale is C set to 1.0. C C 8. if the parameter for the maximum number of cutbacks in the C line search, ncuts, is set. If not, ncuts is set to 10. C C After the optimization, ifcpreprc unscales the final soloution x and the C value of the function, f, at x. C C On entry C C fcn -- The argument containing the name of the user-supplied C subroutine that returns values for the function to be C minimized. C C x -- A double-precision vector of length n. The point at which C the derivative is to be evaluated. C C u -- A double-precision vector of length n. The vector C containing the upper bounds for the n independent C variables. C C l -- A double-precision vector of length n. The vector C containing the lower bounds for the n independent C variables. C C n -- An integer. The dimension of the problem. C C fscale -- A user-supplied double-precision constant used to scale C the function. Fscale should be an approximation to the C maximum size of the absolute value of the function f in C the hyper-box. C C tol -- A user-supplied double-precision constant used for determining C convergence of the algorithm at a given point for a given scale. C The determination of appropriate values of tol is C problem dependent. See the user's guide for help determining C tol. The default value for tol is 1. C C maxit -- A user-supplied integer. Maxit specifies the maximum number C of iterations the algorithm is allowed to take at a given C scale. C C rmaxh -- A user-supplied double-precision constant. Rmaxh is the C first and largest scale the algorithm uses to calculate C the finite difference gradient. Rmaxh must satisfy: C C 0 < rmaxh < 0.5. C C rminh -- A user-supplied double-precision constant. Rminh is a C lower bound for the last and smallest scale the algorithm C uses to calculate the finite difference gradient. C As soon as the scales are reduced to be less than rminh C the algorithm either terminates or restarts. However, C an optimization run is always done with the scale C rmaxh regardless of the value of rminh. Rminh must C satisfy: C C 0 < rminh. C C ncuts -- A user-supplied integer. Ncuts specifies the maximum number C cutbacks in the line search the algorithm may take. C C writ -- A user-supplied integer. Writ specifies the verbosity C of the algorithm: C C writ = 1 the algorithm writes error C messages and the values of C user-supplied parameters C to standard output. C C writ = 2 the algorithm produces the C same information as for write=1 C plus a history of the iterative C process. C C If writ is equal to any other number no output is C produced. C C restart -- A user-supplied integer. If restart is set to 1, C the algorithm restarts after the smallest scale C has been used. The restart uses the last point C generated by the algorithm as the new initial point C and rmaxh as the first scale in the sequence of finite C difference steps. C C C On return C C x -- The initial point scaled for use by iffco. C C xs1 -- A double-precision vector of length n, used for scaling C and unscaling the vector x. C C xs2 -- A double-precision vector of length n, used for scaling C and unscaling the vector x. C C C oops -- An integer. If an upper bound is less than a lower C bound or if the initial point is not in the C hyper-box oops is set to 1 and iffco terminates. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC oops=0 meps = 1.d-16 do 20 i=1,n C C Check if the hyper-box is well-defined. C if(u(i).le.l(i))then if(writ.eq.1.or.writ.eq.2)then write(6,10) i 10 format('Upper bound less than or equal to lower bound ', +'for constraint number',i5.1) end if oops=1 return end if C C Check if the initial iterate is in the hyper-box. C Allow some leeway for nearly equal components. C if(x(i).lt.l(i)-meps.or.x(i).gt.u(i)+meps) then if(writ.eq.1.or.writ.eq.2)then write(6,15)i 15 format('Initial iterate not in bounds for ', +'constraint number',i5.1) end if oops=1 return end if 20 continue C C Check if termination constant tol is set. If not, set tol = 1. C if(tol.le.0.d0)then tol=1.d0 if(writ.eq.1.or.writ.eq.2)then write(6,*)'Default termination constant used. tol=1' end if end if C C Check if maximum scale rmaxh is set. If not, set rmaxh = 0.5. C if(rmaxh.le.0.d0.or.rmaxh.gt.5.d-1)then rmaxh=5.d-1 if(writ.eq.1.or.writ.eq.2)then write(6,*)'Default maximum scale used. rmaxh=0.5' end if end if C C Check if minimum scale rminh is set. If not, set rminh = 1.d-4. C if(rminh.le.0.d0)then rminh=1.d-4 if(writ.eq.1.or.writ.eq.2)then write(6,*)'Default minimum scale used. rminh=1.d-4' end if end if C C Check if maximum number of iterates for each scale maxit is set. C If not, set maxit = 200. C if(maxit.le.0)then maxit=200 if(writ.eq.1.or.writ.eq.2)then write(6,*)'Default maximum number of iterates ', +'used. maxit=200' end if end if C C Check if the constant for scaling the function is set. C If not, set fscale = 1.0. C if(fscale.le.0.d0)then fscale=1.d0 if(writ.eq.1.or.writ.eq.2)then write(6,*)'Default constant for scaling the function' +,' used. fscale=1.d0' end if end if C C Check if the number of cutbacks in the line search has been set. C If not, set ncuts = 10. C if(ncuts.le.0)then ncuts=10 if(writ.eq.1.or.writ.eq.2)then write(6,*)'Default number of cutbacks in ', +'the line search used. ncuts=10' end if end if C C Check if restarts are to be done. C if(restart.ge.1)then if(writ.eq.1.or.writ.eq.2)then write(6,*)'Restarts will be done.' end if else if(writ.eq.1.or.writ.eq.2)then write(6,*)'No restarts will be done.' end if end if C C Write out the parameters being used. C if(writ.eq.1.or.writ.eq.2)then write(6,35) write(6,*) 'Function scale ===> ',fscale write(6,*) 'Tolerance ===> ',tol write(6,*) 'Maximum iterations ===> ',maxit write(6,*) 'Ncuts ===> ',ncuts write(6,*) 'Maximum scale ===> ',rmaxh write(6,*) 'Minimum scale ===> ',rminh write(6,35) 35 format(//) end if C C Description of the output for the optimization. C Only printed if writ = 2. C if(writ.eq.2)then write(6,*)'All data is for the scaled problem.' write(6,*)' ' write(6,*)'K = iteration counter for given scale.' write(6,*)' ' write(6,*)'| X | = two norm of current point.' write(6,*)' ' write(6,*)'F = function value at current point.' write(6,*)' ' write(6,*)'| GF | = two norm of the gradient at the' write(6,*)'current point.' write(6,*)' ' write(6,*)'H = current scale being used.' write(6,*)' ' write(6,*)'CUTS = number of cutbacks in the line search.' write(6,*)' ' write(6,*)'Convergence - The current point satisfies the' write(6,*)'termination criteria at the current scale.' write(6,*)' ' write(6,*)'Line search failed - The line search algorithm' write(6,*)'was unable to find an acceptable new point' write(6,*)'using the current scale.' write(6,35) end if C C Scale the initial iterate so that it is in the unit cube. C do 50 i=1,n xs1(i)=(u(i)-l(i)) xs2(i)=l(i)/xs1(i) x(i)=x(i)/xs1(i)-xs2(i) 50 continue return end subroutine ifcinfcn(fcn,x,c1,c2,fscale,n,func,f) integer n,i,func real*8 x(n),c1(n),c2(n),f,fscale external fcn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCINFCN C C Subroutine ifcinfcn unscales the variable x for use in the user-supplied C function-evaluation subroutine fcn. After fcn returns to ifcinfcn, ifcinfcn C then rescales x for use by iffco. Ifcinfcn also scales the function value C f returned by fcn by dividing f with the user-supplied constant fscale. C Ifcinfcn also updates the function-evaluation counter func. C C On entry C C fcn -- The argument containing the name of the user-supplied C subroutine that returns values for the function to be C minimized. C C x -- A double-precision vector of length n. The point at which C the derivative is to be evaluated. C C xs1 -- A double-precision vector of length n. Used for scaling C and unscaling the vector x by ifcinfcn. C C xs2 -- A double-precision vector of length n. Used for scaling C and unscaling the vector x by ifcinfcn. C C fscale -- A user-supplied double-precision constant used to scale C the function. Fscale should be an approximation to the C maximum size of the absolute value of the function f in C the hyper-box. C C n -- An integer. The dimension of the problem. C C On return C C func -- An integer. The function-evaluation counter. Func is updated C by 1 every time ifcinfcn is called. C C f -- A double-precision scalar. The function value at x scaled for use C by iffco, ifcgrad, or ifclines. C C C Subroutines and Functions C C The subroutine whose name is passed through the argument fcn. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Add one to the function-evaluation counter func. C func=func+1 C C Unscale the variable x. C do 20 i=1,n x(i)=x(i)*c1(i)+c1(i)*c2(i) 20 continue C C Call the function-evaluation subroutine fcn. C call fcn(x,n,f) C C Scale the function value just returned by fcn. C f=f/fscale C C Rescale the variable x. C do 30 i=1,n x(i)=x(i)/c1(i)-c2(i) 30 continue return end subroutine ifcgrad(fcn,x,c1,c2,fscale,fox,n,h,func,gf) integer n,j,func real*8 x(n),fox,gf(n),foxph,z,temp,foxmh +,fscale,c1(n),c2(n),h,fcn external fcn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCGRAD C C Subroutine ifcgrad is a finite difference gradient evaluation subroutine. C Ifcgrad will not take a finite difference step outside of the hyper-box C defined by the constraints on x. Thus, if a forward difference step would C require a function evaluation outside the hyper-box, ifcgrad will not take C it but instead will take the backward difference step. Ifcgrad will always C try to use the centered difference formula if possible. If this is not C possible, it will use either the forward difference formula or the C backward difference formula depending on which one does not violate C the constraints. C C On entry C C x -- A double-precision vector of length n. The point at which C the derivative is to be evaluated. C C xs1 -- A double-precision vector of length n. Used for scaling C and unscaling the vector x by ifcinfcn. C C xs2 -- A double-precision vector of length n. Used for scaling C and unscaling the vector x by ifcinfcn. C C fscale -- A user-supplied double-precision constant used to scale C the function. Fscale should be an approximation to the C maximum size of the absolute value of the function f in C the hyper-box. C C f -- A double-precision scalar. The function value at x. C C n -- An integer. The dimension of the problem. C C h -- A double-precision scalar. The finite difference step being C used to calculate the gradient. C C func -- An integer. The function-evaluation counter. C C On return C C func -- The function-evaluation counter, updated by C the number of calls to the function-evaluation C subroutine needed to calculate the gradient. C C gf -- A double-precision vector of length n. The finite C difference gradient at the point x. C C Subroutines and Functions C C ifcinfcn C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC do 30 j=1,n C C Calculate the forward difference step. C temp=x(j)+h h=temp-x(j) z=x(j) x(j)=x(j)+h C C Check if the forward difference step is in bounds. C if(x(j).gt.1.d0)then C C If the forward difference step is not in bounds, calculate the C backward difference partial for the jth variable. C x(j)=x(j)-h-h call ifcinfcn(fcn,x,c1,c2,fscale,n,func,foxmh) gf(j)=(fox-foxmh)/h goto 20 end if C C Evaluate the function at the forward difference step. C call ifcinfcn(fcn,x,c1,c2,fscale,n,func,foxph) C C Calculate the backward difference step. C x(j)=x(j)-h-h C C Check if the backward difference step is in bounds. C if(x(j).lt.0.d0)then C C If the backward difference step is not in bounds, calculate the C backward difference partial for the jth variable. C gf(j)=(foxph-fox)/h goto 20 end if C C Calculate the centered difference partial for the jth variable. C call ifcinfcn(fcn,x,c1,c2,fscale,n,func,foxmh) gf(j)=(foxph-foxmh)/(2.d0*h) 20 x(j)=z 30 continue return end double precision function ifc2nor(x,n) real*8 x(n),sum integer i,n CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FUNCTION IFC2NOR C C Function ifc2nor calculates the euclidean norm of the vector x. C Eucnor is a double-precision function. C C On entry C C x -- A double-precision vector of length n. The vector C whose norm is being calculated. C C n -- An integer. The dimension of the problem. C C On return C C ifc2nor -- The euclidean norm of x. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Calculate the sum of squared elements of x. C sum=0.d0 do 20 i=1,n 20 sum=x(i)**2+sum C C Calculate the square root of the sum of squares. C C ifc2nor=dsqrt(sum)/dsqrt(dble(n)) ifc2nor=dsqrt(sum/(dble(n))) return end subroutine ifcstep(x,ox,gf,ogf,sr1,s,v,t,b,work,jpvt, +n,fststp,p) integer n,fststp,i,j,info,jpvt(n),denom real*8 x(n),ox(n),gf(n),ogf(n),sr1(n,n),s(n) real*8 b(n,n),work(n),v(n),p(n),t(n) real*8 meps CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCSTEP C C Subroutine ifcstep calculates the negative of a descent direction for the C object function at the current point x in the iterative process. If the C current scale is being used for the first time to calculate the gradient, C the gradient is used and the SR1 approximate Hessian, sr1, is set to the C identity. Otherwise, the SR1 approximate Hessian, sr1, is calculated. sr1 is C then tested for positive definitness by the LINPACK subroutine dchdc. C If the matrix is positive definite, dchdc factors the matrix using a C Cholesky decomposition. The LINPACK subroutine dposl is then used to C determine the negative quasi-Newton step. If the matrix sr1 is not positive C definite, the gradient is used as the negative descent direction. C C On entry C C x -- A double-precision vector of length n. The current point C in the iterative process. C C ox -- A double-precision vector of length n. The previous point C in the iterative process. C C gf -- A double-precision vector of length n. The finite C difference gradient at the point x. C C ogf -- A double-precision vector of length n. The finite C difference gradient at the point ox. C C sr1 -- An upper triangular double-precision matrix of dimension nxn. C The previous SR1 approximate Hessian. C C s -- A double-precision vector of length n. S is a work array C used in ifcstep. C C v -- A double-precision vector of length n. V is a work array C used in ifcstep. C C t -- A double-precision vector of length n. T is a work array C used in ifcstep. C C b -- A double-precision matrix of dimension nxn. B is a work array C used in ifcstep. C C work -- A double-precision vector of length n. Work is a work array C used in ifcstep. C C jpvt -- An integer array of length n. Jpvt is a work array C used in ifcstep. C C n -- An integer. The dimension of the problem. C C fststp -- An integer. Fststp indicates if the scale being used C to calculate the gradient has been changed. If so, C fststp is set to zero. Otherwise, fststp is set to C one. C C On return C C sr1 -- The updated SR1 approximate Hessian. C C p -- A double-precision vector of length n. P is the negative descent C direction for f at x. C C Subroutines C C ifcimat, ifcsr1, ifcswapm, dchdc, dposl C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC meps = 1.d-16 if(fststp.eq.0)then C C If this is the first step to be taken at the current scale, set the C negative descent direction to the gradient, set the SR1 approximate c Hessian to the identity, and return. C call ifcimat(n,sr1) do 20 i=1,n 20 p(i)=gf(i) return else C C Calculate the SR1 update. C do 30 i=1,n s(i)=x(i)-ox(i) 30 v(i)=gf(i)-ogf(i) call ifcsr1(v,s,t,sr1,n,denom) end if C C If the ith constraint is active, zero out the nondiagonal elements of C the SR1 matrix on the ith row and column. C do 40 i=1,n-1 do 40 j=i+1,n if(((dabs(x(i)).le.meps.and.gf(i).ge.0.d0).or.(dabs(x(i)- + 1.d0).le.meps.and.gf(i).le.0.d0)).or.((dabs(x(j)).le. + meps.and.gf(j).ge.0.d0).or.(dabs(x(j)-1.d0).le.meps + .and.gf(j).le.0.d0))) then sr1(i,j)=0.d0 end if 40 continue C C Copy the SR1 matrix to a temporary matrix b for use by the LINPACK C subroutines dchdc and dposl. C call ifcswapm(sr1,n,b) C C Factor the matrix b. C call dchdc(b,n,n,work,jpvt,0,info) C C If the SR1 matrix is not positive definite, set it to the identity C and use the gradient as the negative descent direction. C if(info.lt.n)then call ifcimat(n,sr1) do 50 i=1,n 50 p(i)=gf(i) return end if do 60 i=1,n C C Calculate the new negative quasi-Newton direction. C 60 p(i)=gf(i) call dposl(b,n,n,p) return end subroutine ifcswapm(m1,n,om) integer n,i,k real*8 m1(n,n),om(n,n) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCSWAPM C C Subroutine ifcswapm copies the matrix m1 onto the matrix om. C C C On entry C C m1 -- An upper triangular nxn double-precision matrix. C C n -- An integer. The order of the matrix m1. C C On return C C om -- An upper triangular nxn double-precision matrix. om = m1. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC do 30 k=1,n do 20 i=k,n om(k,i)=m1(k,i) 20 continue 30 continue return end subroutine ifcsr1(y,s,t,sr1,n,denom) integer n,i,j,denom real*8 y(n),s(n),sr1(n,n),t(n),st,ast CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCSR1 C C Subroutine ifcsr1 forms the SR1 approximate Hessian. The SR1 update C is calculated by the following formula: C C sr1=sr1+[(y-sr1*s)(y-sr1*s)^T]/(s^T(y-sr1*s)), C C where sr1 is the previous sr1 matrix, y is the difference between the C gradient at the current point and the gradient at the last point, and C s is the previous step. C C On entry C C y -- A double-precision vector of length n. The difference C between the gradient at the current point and at the C previous point. C C s -- A double-precision vector of length n. The difference C between the current point and the previous point. C C sr1 -- An upper triangular double-precision matrix of dimension nxn. C The previous SR1 approximate Hessian. C C t -- A double-precision vector of length n. T is a work array C used in ifcsr1. C C n -- An integer. The dimension of the problem. C C On return C C sr1 - An upper triangular double-precision matrix of dimension nxn. C The new SR1 approximate Hessian. C C denom -- An integer. If the inner product |s^T(y-sr1*s)| < 10^(-8), C denom is set to 0 and the SR1 update is skipped. C Otherwise, denom = 1 and the SR1 update is computed. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC denom=1 st=0.d0 C C Calculate the matrix vector product sr1*s. C do 45 i=1,n t(i)=0.d0 do 30 j=1,i 30 t(i)=sr1(j,i)*s(j)+t(i) do 40 j=i+1,n 40 t(i)=sr1(i,j)*s(j)+t(i) 45 continue C C Calculate the difference y-sr1*s. C do 50 i=1,n 50 t(i)=y(i)-t(i) C C Calculate the inner product s^T(y-sr1*s). C do 60 i=1,n 60 st=s(i)*t(i)+st C C Check if the inner product s^T(y-sr1*s) is too small in C absolute value. If so, skip the sr1 update. C ast=dabs(st) if(ast.le.1.d-15)then denom=0 goto 999 end if C C Calculate the new sr1 approximate Hessian. C do 70 i=1,n do 70 j=i,n 70 sr1(i,j)=sr1(i,j)+(1.d0/st)*t(i)*t(j) 999 return end subroutine ifcimat(n,m) integer n,i,j real*8 m(n,n) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCIMAT C C Subroutine ifcimat sets the matrix m to the identity matrix. C C On entry C C n -- An integer. The order of the matrix. C C On return C C m -- An nxn double-precision matrix. M is the identity matrix. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Overwrite m with the identity matrix. C do 40 i=1,n do 30 j=1,n m(j,i)=0.d0 30 continue m(i,i)=1.d0 40 continue return end subroutine ifclines(fcn,f,x,gf2,p,xs1,xs2,fscal,h,ncuts,n,func, +lp,linf,y,fp,cutnum) integer maxor C C Set the maximum order for the arrays used. C parameter(maxor=100) integer n,ncuts,linf,i,j,func,onwall,offwall,cutnum real*8 f,x(n),gf2,p(n),xs1(n),xs2(n),alp,y(n),fp,fscal real*8 alpp,alp2p,fpp,talp,enlp,rh,tmpx,h real*8 ifc2nor,fcn,lp(n) external fcn real*8 yold(maxor),meps CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCLINES C C Subroutine ifclines is a cutback line search algorithm that uses a cubic C model of the objective function to determine the step length reduction C parameter. Below is a general outline of a cutback line search algorithm. C C Given a descent direction P for the function to be minimized, f C at the point x: C C 1) Set the step length parameter alp = 1. C 2) if x+alp*P is acceptable, set x = x+alp*P and return to the main C program. C 3) Set alp = r*alp, 0 < r < 1, and return to 2. C C Ifclines determines if x+alp*P is acceptable by using the Armijo rule: C C f(x+alp*P) <= f(x)+10^(-4)*alp*gf(x)^T P, C C where gf(x)^T is the transpose of the gradient at x. C C The reduction factor r is determined in the following manner: C C 1) If the previous step length parameter activated a constraint C that had not been previously activated, l = l/2. C 2) If the previous step length parameter was the first step length C parameter not to activate a new constraint, use a quadratic C model of the function to determine the new reduction factor. C 3) Otherwise, use a cubic model of the function to determine C the new reduction factor. See the user's manual for details C on forming the quadratic and cubic reduction factors. C C The line search is halted in three ways: C C 1) The Armijo rule is satisfied. In this case, the line search C terminates successfully. C 2) The algorithm does not find an acceptable point after the C maximum allowed number of reductions, ncuts, in the line search C have been taken. In this case, the line search terminates C unsuccessfully. C 3) The step length is reduced below 0.1*h, where h is the current C scale being used to determine the gradient. In this case, C the line search terminates unsuccessfully. C C On entry C C fcn -- The argument containing the name of the user-supplied C subroutine that returns values for the function to be C minimized. C C x -- A double-precision vector of length n. The current point C in the iterative process. C C f -- A double-precision scalar. The function value at x. C C gft -- A double-precision scalar. Gft equals the inner product C of the descent direction for f at x and the gradient. C C p -- A double-precision vector of length n. P is the negative descent C direction for f at x. C C xs1 -- A double-precision vector of length n, used for scaling C and unscaling the vector x. C C xs2 -- A double-precision vector of length n, used for scaling C and unscaling the vector x. C C fscale -- A user-supplied double-precision constant used to scale C the function. Fscale should be an approximation to the C maximum size of the absolute value of the function f in C the hyper-box. C C h -- A double-precision scalar. The current scale being used to C calculate the finite difference gradient. C C ncuts -- The maximum number of cutbacks allowed in the line search. C C n -- An integer. The dimension of the problem. C C func -- An integer. The function-evaluation counter. C C lp -- A double-precision vector of length n. Lp is a work C array used by ifclines. C C On return C C func -- The function evaluation-counter, updated by C the number of calls to the function-evaluation C subroutine needed to do the line search. C C linf -- An integer. Linf reports whether the line search succeeded C or failed. If linf = 0 the line search failed. C If linf = 1 the line search succeeded. C C y -- A double-precision vector of length n. If the line search C succeeds, y contains the new point in the iterative C process. If the line search fails, y contains the last C point tried in the line search. In this case y is C overwritten in iffco. C C fp -- A double-precision scalar. Fp contains the value C of the objective function evaluated at y. C C cutnum -- An integer. Cutnum contains the number of cutbacks C taken in the backtracking line search. If the step C x+p is accepted, cutnum = 0; if one reduction is C taken, cutnum = 1, and so on. C C Subroutines C C ifccubcalc, ifcinfcn, ifc2nor C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC meps = 1.d-16 C Initialize yold to something different than y. do 10 i=1,n yold(i) = y(i) + 1 10 continue C C Set the counter for steps that do not activate a new constraint. C offwall=0 C C Set the check for line search success to yes. C linf = 1 C C Begin the line search. C do 100 j=0,ncuts cutnum=j C C Set the step length parameter alp to check the full step. C if(j.eq.0)then alp=1.d0 else C C If a new constraint has been activated, set alp = 0.5*alp. C if(onwall.eq.1) then alp2p=alpp alpp=alp alp=5.d-1*alp else C C If the previous step length parameter was the first step length parameter C not to activate a new constraint, calculate the new step length parameter C using a quadratic model of the function. C if(offwall.eq.1)then talp=alp*gf2/(2.d0*(fp-f+alp*gf2)) if(talp.lt.1.d-1) talp=1.d-1 if(talp.gt.5.d-1) talp=5.d-1 alp2p=alpp alpp=alp alp=talp*alp else C C Calculate the new step length parameter using a cubic model of the function. C alp2p=alpp alpp=alp call ifccubcalc(f,gf2,fp,fpp,alp,alpp,alp2p) end if end if end if C C If the step length is reduced to less than 0.1*h, then return and report C failure in the line search. C if (j.gt.0) then do 20 i=1,n 20 lp(i)=alp*p(i) enlp=ifc2nor(lp,n) rh=1.d-2*h if (enlp.lt.rh) then goto 101 end if end if C C Set the check for new active constraints to off. C onwall = 0 C C Calculate the projection of the new iterate onto the hyper-box. C do 30 i=1,n if(x(i)-alp*p(i).ge.1.d0)then y(i)=1.d0 else if(x(i)-alp*p(i).le.0.d0)then y(i)=0.d0 else y(i)=x(i)-alp*p(i) end if end if C C Check if new constraints have been activated. If so, set the check for new C active constraints to on. C tmpx=x(i)-alp*p(i) if(tmpx.lt.0.d0.and.x(i).ne.0.d0) onwall = 1 if(tmpx.gt.1.d0.and.x(i).ne.1.d0) onwall = 1 30 continue C C Update the counter for steps that do not activate new constraints. C if(onwall.eq.0) offwall=offwall+1 fpp=fp C Check to see if line search came up w/ same point as last time i=1 do while ( dabs(y(i)-yold(i)).le.meps.and.i.lt.n) i=i+1 end do if(i.lt.n.or.dabs(y(n)-yold(n)).gt.meps ) then C Evaluate the function at the new point y. call ifcinfcn(fcn,y,xs1,xs2,fscal,n,func,fp) C Save point. do 40 i=1,n yold(i) = y(i) 40 continue end if C C Check if the new point satisfies the Armijo rule. C If so, return and report success in the line search. C Otherwise, continue. C if(fp.le.f-1.d-4*alp*gf2) return 100 continue C C Set the check for line search success to no. C 101 linf = 0 return end subroutine ifccubcalc(f,gf2,fp,fpp,alp,alpp,alp2p) real*8 f,gf2,fp,fpp,alp,alpp,alp2p,a1,a2,a3,a4,a5 real*8 a6,a7,a,b,talp CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCCUBCALC C C Subroutine ifccubcalc calculates the step length parameter for the cutback C line search algorithm ifclines, using a cubic model of the object function. C See the user's manual for details on the formation of the cubic model. C C On entry C C f -- A double-precision scalar. The value of the object function at x. C C gft -- A double-precision scalar. Gft equals the inner product C of the descent direction for f at x and the gradient. C C fp -- A double-precision scalar. The value of the object function at C the previous point in the line search. C C fpp -- A double-precision scalar. Fpp is the value of the object C function at the point prior to the previous point in the C line search. C C alpp -- A double-precision scalar. The previous step length parameter. C C alp2p -- A double-precision scalar. The step length parameter C prior to the previous step length parameter. C C On return C C alp -- A double-precision scalar. The new step length parameter. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Calculate the new step length parameter using a cubic model of the objective C function. C a1=1.d0/(alpp-alp2p) a2=1.d0/alpp**2 a3=-1.d0/alp2p**2 a4=-alp2p*a2 a5=-alpp*a3 a6=fp-f+gf2*alpp a7=fpp-f+gf2*alp2p a=(a2*a6+a3*a7)*a1 b=(a4*a6+a5*a7)*a1 talp=(-b+dsqrt(b**2+3.d0*a*gf2))/(3.d0*a) C C If the new step length parameter is less than a tenth of the previous C step length parameter, set the new step length parameter equal to a C tenth of the previous step length parameter. C If(talp.lt.1.d-1*alpp) talp=1.d-1*alpp C C If the new step length parameter is more than a half of the previous C step length parameter, set the new step length parameter equal to a C half of the previous step length parameter. C If(talp.gt.5.d-1*alpp) talp=5.d-1*alpp alp=talp return end subroutine ifcprint(x,f,gf,n,k,h,cutnum,conv) integer n,k,cutnum,conv,p real*8 x(n),gf(n),f,h,enx,engf,ifc2nor CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCPRINT C C Subroutine ifcprint writes a history of the optimization iterative C process to standard output. Ifcprint is only called if the user-supplied C integer writ is set to 2 in the call to the subroutine ifcpreprc. C At each iteration of the optimization run, ifcprint prints the C the number of steps taken while using the current scale, the two norm of the C vector x, the value of the object function f at x, the two norm of the C gradient, and the scale currently being used to calculate the gradient. C If the termination criteria are satisfied at the current point with C with the current scale, ifcprint also prints "Convergence". C If the line search at the current point with the current scale succeeds C in finding a new point, ifcprint prints the number of cutbacks needed C to obtain this point. If the line search fails to find a new point, C ifcprint prints "Line search failed". C C On entry C C x -- A double-precision vector of length n. The current point C in the iterative process. C C f -- A double-precision scalar. The function value at x. C C gf -- A double-precision vector of length n. The finite C difference gradient at the point x. C C n -- An integer. The dimension of the problem. C C k -- An integer. The number of steps taken while using the current C scale plus one. C C cutnum -- An integer. The number of cutbacks used in the line search to C obtain the next point in the iterative process. C Cutnum is only used if the line search succeeds. C C conv -- An integer. Conv reports the results of the iterative process C at the current point with the current scale. Conv reports C either convergence, or the failure or success of the line C search. C C FUNCTIONS C C ifc2nor C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Set p to the number of iterates at the current point, scale pair. C p=k-1 C C Calculate the two norms of the current point and gradient. C enx=ifc2nor(x,n) engf=ifc2nor(gf,n) if(conv.eq.1)then C C If the termination criteria was satisfied at the current point, scale pair, C write to standard output the data for the current point, scale pair and C the message "Convergence". C write(6,10)p,enx,f,engf,h 10 format(i3,2x,4(d10.4,2x),'Convergence') else if(conv.eq.0)then C C If the line search failed at the current point, scale pair, C write to standard output the data for the current point, scale pair and C the message "Line search failed". C write(6,15)p,enx,f,engf,h 15 format(i3,2x,4(d10.4,2x),'Line search failed') else C C If the line search succeeded at the current point, scale pair, C write to standard output the data for the current point, scale pair and C the number of cutbacks needed to obtain the new point. C write(6,20)p,enx,f,engf,h,cutnum 20 format(i3,2x,4(d10.4,2x),i3) end if end if return end subroutine ifcpstprc(x,xs1,xs2,f,fscale,n,func,writ) integer n,i,func,writ real*8 x(n),xs1(n),xs2(n),f,fscale CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE IFCPSTPRC C C Subroutine ifcpstprc is the post-processing subroutine for iffco. C Subroutine ifcpstprc unscales the final solution x and the function value at x. C C On entry C C x -- A double-precision vector of length n. The final solution C to the problem. X is still scaled for use by iffco. C C xs1 -- A double-precision vector of length n, used for scaling C and unscaling the vector x. C C xs2 -- A double-precision vector of length n, used for scaling C and unscaling the vector x. C C f -- A double-precision scalar. The function value at x. C F is still scaled for use by iffco. C C fscale -- A double-precision scalar used to unscale f. C C On return C C x -- The unscaled final solution. C C f -- The unscaled function value at x. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C If writ = 2, print the number of function evaluations. C if(writ.eq.2)then write(6,*)'Number of function evaluations: ',func end if C C Unscale the function value at the final point. C f=fscale*f C C Unscale the final iterate. C do 60 i=1,n x(i)=xs1(i)*x(i)+xs1(i)*xs2(i) 60 continue return end