C ================================================================= C File: spg.f C ================================================================= C ================================================================= C Module: Spectral Projected Gradient Method C ================================================================= C Last update of any of the component of this module: C September 10, 2012. c C Users are encouraged to download periodically updated versions of C this code at the TANGO Project web page: C C www.ime.usp.br/~egbirgin/tango/ C Subroutine SPG implements the Spectral Projected Gradient Method C (Version 2: "Feasible continuous projected path") to find a C local minimizers of a given function with convex constraints, C described in C C E.G. Birgin, J.M. Martinez and M. Raydan, "Nonmonotone spectral C projected gradient methods for convex sets", SIAM Journal on C Optimization 10, pp. 1196-1211, 2000. C C The user must supply the external subroutines evalf, evalg and C proj (not necessarily with these names) to evaluate the objective C function and its gradient and to project an arbitrary point onto C the convex feasible region. C C subroutine evalf(n,x,f,flag): Given x, evaluates the function C value f. The user must set flag .ne. 0 if the evaluation is not C possible. C subroutine evalg(n,x,g,flag): Given x, evaluates the gradient g. C The user must set flag .ne. 0 if the evaluation is not possible. C subroutine proj(n,x,flag): Given x, computes its Euclidian C projection onto the feasible set (x is an input/output array in C this subrotuine). The user must set flag .ne. 0 if the projection C is not possible. C This version 10 SEP 2012 by E.G.Birgin, J.M.Martinez and M.Raydan. C Other parameters (i means input, o means output): C C n (i) number of variables C x (i/o) initial guess on input, solution on output C epsopt (i) tolerance for the convergence criterion c m (i) number of delays that defines the non-monotone strategy c (recommended 10, maximum 100) C maxit (i) maximum number of iterations C maxfc (i) maximum number of functional evaluations C iprint (i) controls output level (0 = no print) c work (i) working space (at least 8*n) C f (o) functional value at the solution C gpsupn (o) sup-norm of the projected gradient at the solution C iter (o) number of iterations C fcnt (o) number of functional evaluations C spginfo (o) indicates the reason for stopping C inform (o) indicates an error in an user supplied subroutine C spginfo: C C 0: Small continuous-projected-gradient norm C 1: Maximum number of iterations reached C 2: Maximum number of functional evaluations reached C 3: Too small step in the line search C C spginfo remains unset if inform is not equal to zero on output C inform: C C 0: ok C -90: error in the user supplied evalf subroutine C -91: error in the user supplied evalg subroutine C -92: error in the user supplied proj subroutine C ***************************************************************** C ***************************************************************** subroutine spg(n,x,epsopt,m,maxit,maxfc,iprint,f,gpsupn,iter,fcnt, +spginfo,inform,evalf,evalg,proj,work) implicit none C SCALAR ARGUMENTS integer fcnt,inform,iter,iprint,m,maxfc,maxit,n,spginfo double precision epsopt,f,gpsupn C ARRAY ARGUMENTS double precision x(n),work(3*n+m) C EXTERNAL SUBROUTINE ARGUMENTS external evalf,evalg,proj call spgb(n,x,epsopt,m,maxit,maxfc,iprint,f,gpsupn,iter,fcnt, +spginfo,inform,evalf,evalg,proj,work(1),work(n+1),work(2*n+1), +work(1),work(n+1),work(3*n+1)) C The call sequence above is o.k. In fact, the pair of array C parameters gp and d are not used at the same time and, in C consequence, can share the same memory location starting at C work(1). The same is true for the pair of array parameters C g and xnew and the memory location starting at work(n+1). end C ***************************************************************** C ***************************************************************** subroutine spgb(n,x,epsopt,m,maxit,maxfc,iprint,f,gpsupn,iter, +fcnt,spginfo,inform,evalf,evalg,proj,d,g,xbest,gp,xnew,lastfv) implicit none C PARAMETERS double precision lmax,lmin parameter ( lmin = 1.0d-30 ) parameter ( lmax = 1.0d+30 ) C SCALAR ARGUMENTS integer fcnt,inform,iter,iprint,m,maxfc,maxit,n,spginfo double precision epsopt,f,gpsupn C ARRAY ARGUMENTS double precision d(n),g(n),gp(n),lastfv(0:m-1),x(n),xbest(n), + xnew(n) C EXTERNAL SUBROUTINE ARGUMENTS external evalf,evalg,proj C LOCAL SCALARS integer i,lsinfo double precision alpha,dtd,dty,fbest,fnew,gtd,lambda C EXTERNAL SUBROUTINES external ls,sevalf,sevalg,sproj C INTRINSIC FUNCTIONS intrinsic abs,max,min,mod C ================================================================== C Initialization C ================================================================== C Print problem information if ( iprint .gt. 0 ) then write(* ,fmt=1000) write(10,fmt=1000) write(* ,fmt=1010) n write(10,fmt=1010) n end if C Set some initial values: C error tracker inform = 0 C for counting number of iterations as well as functional evaluations iter = 0 fcnt = 0 C for the non-monotone line search do i = 0,m - 1 lastfv(i) = - 1.0d+99 end do C Project initial guess call sproj(n,x,inform,proj) if ( inform .ne. 0 ) return C Compute function and gradient at the initial point call sevalf(n,x,f,inform,evalf) if ( inform .ne. 0 ) return fcnt = fcnt + 1 call sevalg(n,x,g,inform,evalg) if ( inform .ne. 0 ) return C Compute continuous-project-gradient sup-norm do i = 1,n gp(i) = x(i) - g(i) end do call sproj(n,gp,inform,proj) if (inform .ne. 0) return gpsupn = 0.0d0 do i = 1,n gpsupn = max( gpsupn, abs( gp(i) - x(i) ) ) end do C Compute spectral steplength if ( gpsupn .ne. 0.0d0) then lambda = min( lmax, max( lmin, 1.0d0 / gpsupn ) ) else lambda = 0.0d0 end if C Store functional value for the non-monotone line search lastfv(0) = f C Initiate best solution and functional value fbest = f do i = 1,n xbest(i) = x(i) end do C ================================================================== C Main loop C ================================================================== 100 continue C Print iteration information if ( iprint .gt. 0 ) then if ( mod(iter,10) .eq. 0 ) then write(* ,fmt=1020) write(10,fmt=1020) end if write(* ,fmt=1030) iter,f,gpsupn write(10,fmt=1030) iter,f,gpsupn end if open(20,file='spg-tabline.out') write(20,fmt=1040) n,iter,fcnt,f,gpsupn close(20) C ================================================================== C Test stopping criteria C ================================================================== C Test whether the continuous-projected-gradient sup-norm C is small enough to declare convergence if ( gpsupn .le. epsopt ) then spginfo = 0 if ( iprint .gt. 0 ) then write(* ,fmt=1030) iter,f,gpsupn write(10,fmt=1030) iter,f,gpsupn write(*, 1100) write(10,1100) end if go to 200 end if C Test whether the number of iterations is exhausted if ( iter .ge. maxit ) then spginfo = 1 if ( iprint .gt. 0 ) then write(*, 1110) write(10,1110) end if go to 200 end if C Test whether the number of functional evaluations is exhausted if ( fcnt .ge. maxfc ) then spginfo = 2 if ( iprint .gt. 0 ) then write(*, 1120) write(10,1120) end if go to 200 end if C ================================================================== C Iteration C ================================================================== iter = iter + 1 C Compute search direction do i = 1,n d(i) = x(i) - lambda * g(i) end do call sproj(n,d,inform,proj) if (inform .ne. 0) return dtd = 0.0d0 gtd = 0.0d0 do i = 1,n d(i) = d(i) - x(i) dtd = dtd + d(i) ** 2 gtd = gtd + g(i) * d(i) end do C Perform safeguarded quadratic interpolation along the spectral C continuous projected gradient call ls(n,x,f,d,gtd,m,lastfv,maxfc,iprint,fcnt,alpha,xnew,fnew, +lsinfo,inform,evalf) if ( inform .ne. 0 ) return if ( lsinfo .eq. 1 ) then spginfo = 3 if ( iprint .gt. 0 ) then write(*, 1130) write(10,1130) end if go to 200 end if if ( lsinfo .eq. 2 ) then spginfo = 2 if ( iprint .gt. 0 ) then write(*, 1120) write(10,1120) end if go to 200 end if C Set new functional value and save it for the non-monotone line C search, and set new iterate f = fnew lastfv(mod(iter,m)) = f do i = 1,n x(i) = xnew(i) end do C Best solution and functional value if ( f .lt. fbest ) then fbest = f do i = 1,n xbest(i) = x(i) end do end if C Compute gradient at the new iterate call sevalg(n,x,g,inform,evalg) if ( inform .ne. 0 ) return C Compute spectral steplength dty = - gtd do i = 1,n dty = dty + d(i) * g(i) end do if ( dty .le. 0.0d0 ) then lambda = lmax else lambda = max( lmin, min( alpha * dtd / dty, lmax ) ) end if C Compute continuous-project-gradient sup-norm do i = 1,n gp(i) = x(i) - g(i) end do call sproj(n,gp,inform,proj) if (inform .ne. 0) return gpsupn = 0.0d0 do i = 1,n gpsupn = max( gpsupn, abs( gp(i) - x(i) ) ) end do C ================================================================== C Iterate C ================================================================== go to 100 C ================================================================== C End of main loop C ================================================================== 200 continue C ================================================================== C Write statistics C ================================================================== if ( iprint .gt. 0 ) then write(* ,fmt=2000) iter,fcnt,f,gpsupn write(10,fmt=2000) iter,fcnt,f,gpsupn end if C ================================================================== C Finish returning the best point C ================================================================== f = fbest do i = 1,n x(i) = xbest(i) end do C ================================================================== C NON-EXECUTABLE STATEMENTS C ================================================================== 1000 format(/,1X,78('='), + /,1X,'This is the SPECTRAL PROJECTED GRADIENT (SPG) for ', + 'for convex-constrained',/,1X,'optimization. If you ', + 'use this code, please, cite:',/, + /,1X,'E. G. Birgin, J. M. Martinez and M. Raydan, ', + 'Nonmonotone spectral projected',/,1X,'gradient ', + 'methods on convex sets, SIAM Journal on ', + 'Optimization 10, pp.',/,1X,'1196-1211, 2000, and',/, + /,1X,'E. G. Birgin, J. M. Martinez and M. Raydan, ', + 'Algorithm 813: SPG - software',/,1X,'for ', + 'convex-constrained optimization, ACM Transactions ', + 'on Mathematical',/,1X,'Software 27, pp. 340-349, ', + '2001.',/,1X,78('=')) 1010 format(/,1X,'Entry to SPG.', + /,1X,'Number of variables: ',I7) 1020 format(/,4X,'ITER',10X,'F',8X,'GPSUPN') 1030 format( 1X,I7,1X,1P,D16.8,1X,1P,D7.1) 1040 format( 1X,I7,1X,I7,1X,I7,1X,1P,D16.8,1X,1P,D7.1,1X,'(Abnormal ', + 'termination. Probably killed by CPU time limit.)') 1100 format(/,1X,'Flag of SPG: Solution was found.') 1110 format(/,1X,'Flag of SPG: Maximum of iterations reached.') 1120 format(/,1X,'Flag of SPG: Maximum of functional evaluations ', + 'reached.') 1130 format(/,1X,'Flag of SPG: Too small step in the line search.', + /,1X,'Probably, an exaggerated small norm of the ', + 'continuous projected gradient', + /,1X,'is being required for declaring convergence.') 2000 format(/,1X,'Number of iterations : ',9X,I7, + /,1X,'Number of functional evaluations : ',9X,I7, + /,1X,'Objective function value : ',1P,D16.8, + /,1X,'Sup-norm of the projected gradient : ',9X,1P,D7.1) end C ***************************************************************** C ***************************************************************** subroutine ls(n,x,f,d,gtd,m,lastfv,maxfc,iprint,fcnt,alpha,xnew, +fnew,lsinfo,inform,evalf) implicit none C SCALAR ARGUMENTS integer inform,iprint,lsinfo,maxfc,fcnt,m,n double precision alpha,f,fnew,gtd C ARRAY ARGUMENTS double precision d(n),lastfv(0:m-1),x(n),xnew(n) C EXTERNAL SUBROUTINE ARGUMENTS external evalf C Nonmonotone line search with safeguarded quadratic interpolation C lsinfo: C C 0: Armijo-like criterion satisfied C 1: Too small step in the line search C 2: Maximum number of functional evaluations reached C PARAMETERS double precision gamma parameter ( gamma = 1.0d-04 ) C LOCAL SCALARS integer i double precision atemp,fmax C EXTERNAL SUBROUTINES external sevalf C INTRINSIC FUNCTIONS intrinsic abs,max C Initiate fmax = lastfv(0) do i = 1,m - 1 fmax = max( fmax, lastfv(i) ) end do alpha = 1.0d0 do i = 1,n xnew(i) = x(i) + alpha * d(i) end do call sevalf(n,xnew,fnew,inform,evalf) if ( inform .ne. 0 ) return fcnt = fcnt + 1 C Main loop 100 continue C Print iteration information if ( iprint .gt. 1 ) then write(* ,fmt=3000) fnew,alpha write(10,fmt=3000) fnew,alpha end if C Test stopping criteria if ( fnew .le. fmax + gamma * alpha * gtd ) then lsinfo = 0 return end if C set the "0.0d0" in the line below to a value that should be C considered a "too small step" in your application. if ( alpha .le. 0.0d0 ) then lsinfo = 1 return end if if (fcnt .ge. maxfc) then lsinfo = 2 return end if C Safeguarded quadratic interpolation if ( alpha .le. 0.1d0 ) then alpha = alpha / 2.0d0 else atemp = ( - gtd * alpha ** 2 ) / + ( 2.0d0 * ( fnew - f - alpha * gtd ) ) if ( atemp .lt. 0.1d0 .or. atemp .gt. 0.9d0 * alpha ) then atemp = alpha / 2.0d0 end if alpha = atemp end if C New trial do i = 1,n xnew(i) = x(i) + alpha * d(i) end do call sevalf(n,xnew,fnew,inform,evalf) if ( inform .ne. 0 ) return fcnt = fcnt + 1 C Iterate go to 100 C NON-EXECUTABLE STATEMENTS 3000 format(9X,1P,D16.8,9X,'(Step = ',1P,D7.1,')') end C ***************************************************************** C ***************************************************************** subroutine sevalf(n,x,f,inform,evalf) implicit none C SCALAR ARGUMENTS integer inform,n double precision f C ARRAY ARGUMENTS double precision x(n) C EXTERNAL SUBROUTINE ARGUMENTS external evalf C LOCAL SCALARS integer flag C EXTERNAL SUBROUTINES external reperr call evalf(n,x,f,flag) C This is true if f if Inf, - Inf or NaN if ( .not. f .gt. - 1.0d+99 .or. .not. f .lt. 1.0d+99 ) then f = 1.0d+99 end if if ( flag .ne. 0 ) then inform = - 90 call reperr(inform) return end if end C ***************************************************************** C ***************************************************************** subroutine sevalg(n,x,g,inform,evalg) implicit none C SCALAR ARGUMENTS integer inform,n C ARRAY ARGUMENTS double precision g(n),x(n) C EXTERNAL SUBROUTINE ARGUMENTS external evalg C LOCAL SCALARS integer flag C EXTERNAL SUBROUTINES external reperr call evalg(n,x,g,flag) if ( flag .ne. 0 ) then inform = - 91 call reperr(inform) return end if end C ***************************************************************** C ***************************************************************** subroutine sproj(n,x,inform,proj) implicit none C SCALAR ARGUMENTS integer inform,n C ARRAY ARGUMENTS double precision x(n) C EXTERNAL SUBROUTINE ARGUMENTS external proj C LOCAL SCALARS integer flag C EXTERNAL SUBROUTINES external reperr call proj(n,x,flag) if ( flag .ne. 0 ) then inform = - 92 call reperr(inform) return end if end C ****************************************************************** C ****************************************************************** subroutine reperr(inform) implicit none C SCALAR ARGUMENTS integer inform if ( inform .eq. -90 ) then write(* ,fmt=100) 'EVALF' write(10,fmt=100) 'EVALF' else if ( inform .eq. -91 ) then write(* ,fmt=100) 'EVALG' write(10,fmt=100) 'EVALG' else if ( inform .eq. -92 ) then write(* ,fmt=100) 'PROJ ' write(10,fmt=100) 'PROJ ' end if C NON-EXECUTABLE STATEMENTS 100 format(/,1X,'*** There was an error in the user supplied ', + 'subroutine ',A10,' ***',/) end