C ================================================================= C File: c-grasp.f C ================================================================= program c_gencan C This program try to find some y such that f(y) <= f(x) by C using the C-GRASP algorithm, described in: C C M.J. HIRSCH, P.M. PARDALOS, M.G.C. RESENDE. SPEEDING UP C CONTINUOUS GRASP. AT&T Labs Research Technical Report TD-6U2P2H, C AT&T Shannon Laboratory, 180 Park Avenue, Florham Park, C NJ 07932 USA, September 2006. C C and C C M.J. HIRSCH, C.N. Meneses, P.M. PARDALOS, M.G.C. RESENDE. Global C optimization by continuous GRASP. AT&T Labs Research Technical C Report TD-6MPUV9. AT&T Labs Research, Florham Park, NJ 07932, C March 8, 2006. C PARAMETERS integer MaxIter,nmax,mmax double precision INFINITO,alpha parameter ( MaxIter = 2000 ) parameter ( INFINITO = 10000000000000000.0d0 ) parameter ( alpha = 0.3d0 ) parameter ( nmax = 500000 ) parameter ( mmax = 500000 ) C LOCAL SCALARS real etime,ntime,btime,ctime integer n,m,i,iter,itbest,bfound, + nsamples,bnsamples,totevalf,bevalf,totevalg,bevalg,fcnmb double precision seed,drand,h,BestSolValue,SolValue,difSol, + bdifSol,vdummy,h_s,h_e logical impr_c,impr_l C LOCAL ARRAYS real dum(2) double precision BestSolution(nmax),y(nmax), + l(nmax),u(nmax) C GENCAN COMPATIBILITY C GENCAN'S SCALARS logical checkder character * 6 precond integer gtype,hptype,iprint,maxoutit,maxtotfc,maxtotit,ncomp, + inform,outiter,totcgcnt,totfcnt,totgcnt,totiter double precision epsfeas,epsopt,nalpsupn,snorm real time C GENCAN'S ARRAYS integer wi1(nmax) logical equatn(mmax),linear(mmax) double precision lambda(mmax),rho(mmax),wd1(mmax), + wd2(mmax),wd3(nmax),wd4(mmax),wd5(mmax),wd6(nmax), + wd7(nmax),wd8(nmax),wd9(nmax),wd10(nmax),wd11(nmax), + wd12(nmax),wd13(nmax),wd14(nmax),wd15(nmax),wd16(nmax), + wd17(nmax),wd18(nmax),wd19(nmax) C EXTERNAL SUBROUTINES external solver C SETUP PROBLEM call inip(n,y,l,u,m,lambda,rho,equatn,linear) C SET LOCAL SCALARS h_s = max(1.0d0, (u(1)-l(1)) / 1000.0d0 ) h_e = 0.1d0 * h_s C write(*,*) '>>>> Inicio <<<< | hs,he = ',h_s,h_e write(*,*) " " write(*,*) " " write(*,0020) write(*,0021) write(*,0022) SolValue = INFINITO BestSolValue = INFINITO C SET LOCAL COUNTERS totevalf = 0 totevalg = 0 nsamples = 0 bfound = 0 C SET TIME CHECK COUNTER ctime = etime(dum) C SET SOME SOLVER ARGUMENTS gtype = 0 hptype = 6 precond = 'NONE' C precond = 'QNCGNA' checkder = .false. epsfeas = 1.0d-04 epsopt = 1.0d-04 maxoutit = 50 maxtotit = 1000000 maxtotfc = 5 * maxtotit iprint = 0 ncomp = 5 C START METHOD iter = 1 C TEST STOP CRITERIA do while ( iter .le. MaxIter ) C SET INITIAL RANDOM POINT seed = iter vdummy = drand(seed) do i = 1,n y(i) = l(i) + ( u(i) - l(i) ) * drand(seed) end do nsamples = nsamples + 1 h = h_s C TEST GRID DENSITY do while( ( h .ge. h_e ) .and. + ( iter .le. MaxIter ) ) impr_c = .false. impr_l = .false. C C-GRASP-GENCAN GLOBAL PHASE call ConstructGreedyRandomizedSolution(y,n,h,l,u,seed, + alpha,fcnmb,impr_c) call evalf(n,y,difSol,inform) C C-GRASP-GENCAN LOCAL PHASE - CALL OPTIMIZATION SOLVER call solver(n,y,l,u,m,lambda,rho,equatn,linear,gtype, + hptype,precond,checkder,epsfeas,epsopt,maxoutit,maxtotit, + maxtotfc,iprint,ncomp,SolValue,snorm,nalpsupn,outiter,totiter, + totfcnt,totgcnt,totcgcnt, + time,inform,wi1,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8,wd9,wd10,wd11, + wd12,wd13,wd14,wd15,wd16,wd17,wd18,wd19) totevalf = totevalf + fcnmb + totfcnt totevalg = totevalg + totgcnt C TEST SOLUTION IMPROVEMENT IN LOCAL PHASE difSol = SolValue - difSol if( abs(difSol) .gt. 0.0d0 ) then impr_l = .true. end if C UPDATES THE BEST SOLUTION ntime = etime(dum) if( BestSolValue .gt. + ( SolValue + 0.0001d0 * abs( SolValue ) + + 1.0E-6 ) ) then do i = 1,n BestSolution(i) = y(i) end do itbest = iter btime = ntime bevalf = totevalf bevalg = totevalg bnsamples = nsamples bdifSol = difSol bfound = bfound + 1 write(*,0023) bfound,iter,SolValue,btime,bdifSol BestSolValue = SolValue end if C CHECK TIME SPENT if( ntime - ctime .gt. 60.0d0 ) then ctime = ntime write(*,0010) ctime,iter,h end if C MAKES THE GRID MORE DENSE if( .not.( impr_c ) .and. .not.( impr_l ) ) then h = h / 2.0d0 end if C ITERATE iter = iter + 1 end do end do C WRITES THE BEST SOLUTION FOUND write(*,*) " " write(*,*) " " write(*,*) " " write(*,0030) write(*,0031) BestSolValue,BestSolValue write(*,0035) write(*,0032) btime,ntime write(*,0033) itbest,bnsamples write(*,0034) bevalf,bevalg C write(*,0035) C write(*,*) " " C do i = 1,n C write(*,0036) i,BestSolution(i) C end do write(*,0037) write(*,*) " " write(*,*) '* Original experiments executed with G77. This ', + 'result may reach a diferent soluction since G77 and ', + 'gfortran are not fully compatible. *' C NON-EXECUTABLE STATEMENTS 0010 format("---- TIME : ",F7.2," -- Iter : ",I4, + " -- h : ",F9.4," ----") 0020 format("-----------------< Best solution found >", + "----------------") 0021 format(" # Iter f* t*", + " GENCAN Improv") 0022 format("--- ----- ------------------ -----", + " -----------------") 0023 format(I3,2X,I5,2X,F18.10,2X,F4.2,2X,F18.10) 0030 format("******************** FINAL SOLUTION ", + "********************") 0031 format("f(x*) : ",D25.16," ( ",F18.10," )" ) 0032 format("time(x*) : ",F12.5," | total time : ",F12.5 ) 0033 format("iter(x*) : ",I4,8X," | #x0 : ",I7 ) 0034 format("eval_f : ",I7,5X," | eval_g : ",I7 ) 0035 format("------------------------------------", + "--------------------") 0036 format("x( ",I4," ) = ",D25.16) 0037 format("************************************", + "********************") end; C ****************************************************************** C ****************************************************************** subroutine ConstructGreedyRandomizedSolution(x,n,h,l,u,seed, + alpha,totevalf,impr_c) C This subroutine try to find some y such that f(y) <= f(x) by C using the C-GRASP algorithm, described in: C C M.J. Hirsch, P.M. Pardalos, and M.G.C. Resende. C Speeding up continuous GRASP. C Technical Report TD-6U2P2H, AT&T Labs Research, Florham Park, NJ 07932, C 2006. C C M.J. Hirsch, C.N. Meneses, P.M. Pardalos, and M.G.C. Resende. C Global optimization by continuous GRASP. C Optimization Letters, 1:201–212, 2007. C C C.N. MENESES, P.M. PARDALOS, AND M.G.C. RESENDE. GRASP FOR C NONLINEAR OPTIMIZATION. AT&T Labs Research Technical C Report TD-6DUTRG. C C Parameters of the subroutine: C C On Entry: C C x double precision x(n) : Initial point x^0 C C n integer : number of variables C C h double precision : grid step length C C l integer : intervals lower bound C C u integer : intervals upper bound C C seed double precision : randomic factor (used by drand) C C alpha double precision : heuristic RCL constant C C On Return: C C x double precision x(n) : min(x^0,y) C C seed double precision : randomic factor (used by drand) C C totevalf integer : total functions evaluations C C impr_c boolean : true if y < x^0; false, otherwise C C PARAMETERS double precision INFINITO parameter ( INFINITO = 10000000000000000.0d0 ) C SCALAR ARGUMENTS integer n,totevalf,fcnmb double precision seed,alpha,h logical impr_c,ReUse C ARRAY ARGUMENTS double precision x(n),l(n),u(n) C LOCAL SCALARS integer i,j,k,RandomlySelectElement double precision vmax,vmin,xibest C LOCAL ARRAYS integer S(n+1),RCL(n+1) double precision z(n),g(n),y(n) C SET LOCAL SCALARS totevalf = 0 fcnmb = 0 do i = 1,n S(i) = i RCL(i) = -1 y(i) = x(i) end do S(n+1) = n RCL(n+1) = 0 C TRY TO FIND SOME POINT z SUCH THAT f(x) < f(z) ReUse = .false. do while( S(n+1) .gt. 0 ) vmin = INFINITO vmax = ( -1.0d0 ) * INFINITO do i = 1,S(n+1) C MINIMIZE COMPONENT y( S(i) ) WITH STEP h if( .not.( ReUse ) ) then C Minimize(x,h,indx,n,l,u,xibest,fbest,m) call Minimize(y,h,S(i),n,l,u,xibest,g(S(i)),fcnmb) totevalf = totevalf + fcnmb z(S(i)) = xibest end if if( vmin .gt. g(S(i)) ) then vmin = g(S(i)) end if if( vmax .lt. g(S(i)) ) then vmax = g(S(i)) end if end do C SET SOLUTION CANDIDATES RCL(n+1) = 0 do i=1,S(n+1) if( g( S(i) ) .le. ( vmin + alpha * ( vmax - vmin ) ) ) + then RCL(n+1) = RCL(n+1) + 1 RCL( RCL(n+1) ) = i end if end do C RANDOMLY SELECT A RCL ELEMENT seed = ( S(n+1) + seed ) / 2.0d0 j = RandomlySelectElement(seed,1,RCL(n+1)) C UPDATES x,y if( y( RCL(j) ) .eq. z( RCL(j) ) ) then ReUse = .true. else x( RCL(j) ) = z( RCL(j) ) y( RCL(j) ) = z( RCL(j) ) ReUse = .false. impr_c = .true. end if C S = S \ {j} k = S( RCL(j) ) S( RCL(j) ) = S( S(n+1) ) S( S(n+1) ) = k S(n+1) = S(n+1) - 1 end do end C ****************************************************************** C ****************************************************************** integer function RandomlySelectElement(seed,l,u) C This subroutine returns a random integer between l and u C Parameters of the subroutine: C C On Entry: C C seed double precision C C l integer : intervals lower bound C C u integer : intervals upper bound C C On Return: C C function integer : returns i \in [l,u] C C seed double precision C C SCALAR ARGUMENTS integer l,u double precision seed C LOCAL SCALARS integer delta double precision drand,nran delta = u - l + 1 nran = drand(seed) nran = drand(seed) RandomlySelectElement = max(nint( nran * delta ),1) return end; C ****************************************************************** C ****************************************************************** subroutine Minimize(x,h,indx,n,l,u,xibest,fbest,totevalf) C C Parameters of the subroutine: C C On Entry: C C x double precision x(n), C initial point, C C h double precision, C grid density C C indx integer, C index to be moved, C C n integer, C number of variables, C C l double precision l(n), C lower bounds on x, C C u double precision u(n), C upper bounds on x, C C On Return C C xibest double precision, C best position found C C fbest double precision C best f's value (with xibest) C PARAMETERS double precision INFINITO parameter ( INFINITO = 10000000000000000.0d0 ) C SCALAR ARGUMENTS integer n,totevalf,indx double precision h,xibest,fbest C ARRAY ARGUMENTS double precision x(n),l(n),u(n) C LOCAL SCALARS integer i,err double precision vmin,vle C LOCAL ARRAYS double precision t(n) C SET LOCAL SCALARS totevalf = 0 do i = 1,n t(i) = x(i) end do C SET STARTUP MIN VALUE call evalf(n,t,vle,err) totevalf = totevalf + 1 vmin = vle xibest = t(indx) t(indx) = l(indx) C SEARCH FOR A BEST VALUE do while ( t(indx) .le. u(indx) ) call evalf(n,t,vle,err) totevalf = totevalf + 1 if( vmin .gt. vle ) then vmin = vle xibest = t(indx) end if t(indx) = t(indx) + h end do fbest = vmin end