*-------------------------------------------------------------- * March 14, 2005 * =============== * * Program associated to the paper: * * A New Gradient Descent Method with an Anticipative Scalar * Approximation of Hessian for Unconstrained Optimization. * * * The test problems are in FUNC.FOR * * This is a draft code ! * There is no warranty ! *-------------------------------------------------------------- parameter(ia=20000) real*8 x(ia), xnew(ia), grad(ia), gradnew(ia), d(ia) real*8 bs(ia), by(ia), ss,sy real*8 eps,epsf, alfa,s, fx, fxnew, ginf, timp, fxmin real*8 step, stepini, stepmax, stepmin, stepar(25010) real*8 ps, delta, vv, eta, gradn real*8 gama,gamab logical igama, igamab integer fgcnt, fgcntt integer*4 gh,gm,gs,gc, ght,gmt,gst,gct, timpexp character*35, numef, fnumef(100) * * Input file open(unit=7,file='funcdat.dat',status='old') * Output files open(unit=1,file='tempbb.out',status='unknown') open(unit=2,file='tempbb.rez',status='unknown') c open(unit=3,file='gama.rez',status='unknown') c open(unit=4,file='fxnew.rez',status='unknown') open(unit=5,file='step.rez',status='unknown') eps = 10.d0**(-6) epsf = 10.d0**(-20) *-------------------------------------- Parameters in backtracking alfa=0.0001d0 s=0.8d0 delta=10.d0 maxiter=5000 c-------------------------------------- Logic: c igama = .true. (NA algorithm) c igamab= .true. (BB algorithm) ig=2 c ig=1 igama = .true. NA (Anticipative) c ig=2 igamab = .true. BB (Barzilai-Borwein) if(ig.eq.1) then igama = .true. igamab = .false. end if if(ig.eq.2) then igama = .false. igamab = .true. end if c * c *** c ******* c ********* *---------------------------------------------------------------------- * if(igama)write(2,947) 947 format(2x,'March 14, 2005 *** NEW GRADIENT ALGORITHM ***') if(igamab)write(2,948) 948 format(2x,'March 14, 2005 *** BARZILAI-BORWEIN ALGORITHM ***') write(2,949) 949 format(2x,' n',2x,' iter',2x,' fgcnt',2x,'time(c)') do i=1,69 read(7,21) numef fnumef(i)=numef 21 format(a35) end do *--------------------------------------------- Here Start Experiments do nexp =1,1 numef=fnumef(nexp) if(igama) then write(1,501) nexp, numef write(5,501) nexp, numef write(*,501) nexp, numef 501 format(/,1x,i3,2x,'New Gradient Algorithm:', * 1x,a35,'Function',/) end if if(igamab) then write(1,502) nexp, numef write(5,502) nexp, numef write(*,502) nexp, numef 502 format(/,1x,i3,2x,'Barzilai-Borwein Algorithm:', * 1x,a35,'Function',/) end if if(igama) write(1,504) 504 format(6x,'n',6x,'iter',4x,'fgcnt'2x,'time(c)',11x,'fxnew', * 16x,'ginf',8x,'ngama',8x,'stepinimin',11x,'stepinimax') if(igamab) write(1,505) 505 format(6x,'n',6x,'iter',4x,'fgcnt'2x,'time(c)',11x,'fxnew', * 16x,'ginf',17x,'stepmin',14x,'stepmax') if(igama)write(1,546) 546 format(1x,126('-')) if(igamab)write(1,547) 547 format(1x,121('-')) c ************** itert = 0 fgcntt = 0 timp = 0.d0 *----------------------------------------------- Here start minimization do n = 1000, 10000, 1000 *----------------------------------------- Initial point call inipoint(n,x,nexp) *----------------------------------------- call gettim(gh,gm,gs,gc) * iter=1 fgcnt=0 ngama =0 *----------------------------------------------------------------------- * ===== Beginning the algorithms ===== call evalfg(n,x,fx,grad,nexp) fgcnt=fgcnt+1 fxmin=fx c write(4,971) fx c971 format(4x,e20.13,',') *------------------------------------- Direction do i=1,n d(i) = -grad(i) end do ps= 0.d0 do i=1,n ps = ps + d(i)*grad(i) end do call back(n,x,d,fx,ps,alfa,s,nexp, step,fgcnt) gradn=0.d0 do i=1,n gradn = gradn + grad(i)*grad(i) end do * ------------------------------------ Begin Iterations 1 continue do i=1,n xnew(i) = x(i) - step * grad(i) end do call evalfg(n,xnew,fxnew,gradnew,nexp) fgcnt=fgcnt+1 if(fxnew .lt. fxmin) fxmin=fxnew c write(4,971) fxnew *-------------------------------------------> gama Barzilai & Borwein if(igamab) then do i=1,n bs(i) = xnew(i)-x(i) by(i) = gradnew(i) - grad(i) end do sy=0.d0 ss=0.d0 do i=1,n sy = sy + bs(i)*by(i) ss = ss + bs(i)*bs(i) end do gamab = sy/ss write(5,555) gamab 555 format(5x,'gamab=',e20.13) stepini = 1.d0/gamab end if *-------------------------------------------> End gama Barzilai & Borwein *-------------------------------------------> gama given by Function values if(igama) then if(fxnew-fx+step*gradn .lt. 0.d0) then ngama=ngama+1 vv = (fx - fxnew - step*gradn)/gradn eta = delta + vv step = step + eta end if gama = 2.d0*(fxnew-fx+step*gradn)/(gradn*step*step) stepini = 1.d0/gama end if *-------------------------------------------> End gama given by Function values c write(3,2445) gama-gamab c2445 format(3x,e20.13,',') *------------------------------------------- *-------------------------------------------> Determine the initial stepsize if(igama) stepini = 1.d0/gama if(igamab) stepini = 1.d0/gamab *------------------------------------------- do i=1,n d(i) =-gradnew(i) end do ps=0.d0 do i=1,n ps = ps + d(i)*gradnew(i) end do *------------------------------------------------- Backtracking step=stepini stepar(iter) = stepini call backnew(n,xnew,d,fxmin,ps,alfa,s,stepini,nexp, step,fgcnt) *---------------------------------------------------------------------- * ------------------------------------------------ STOP criterion ginf=dabs(gradnew(1)) do i=2,n if(dabs(gradnew(i)) .gt. ginf) ginf=dabs(gradnew(i)) end do if(ginf .le. eps .or. * dabs(step*ps) .le. epsf*dabs(fxnew)) go to 100 iter=iter+1 if(iter .gt. maxiter) go to 100 *----------------------------------------------------------------- *---------------------------------------- Updating for a new iteration do i=1,n x(i) = xnew(i) grad(i) = gradnew(i) end do fx=fxnew gradn=0.d0 do i=1,n gradn = gradn + grad(i)*grad(i) end do go to 1 *---------------------------------------- go to 1 100 continue cc call gettim(ght,gmt,gst,gct) call exetim(gh,gm,gs,gc, ght,gmt,gst,gct) *-- timpexp = ght*360000 + gmt*6000 + gst*100 + gct stepmin = stepar(1) stepmax = stepar(1) do i=2,iter if(stepar(i) .le. stepmin) stepmin=stepar(i) if(stepar(i) .gt. stepmax) stepmax=stepar(i) end do *---------------------------------------- Write some information if(igama) then if(ngama .eq. 0) then write(1,901)n,iter,fgcnt,timpexp,fxnew,ginf,stepmin,stepmax write(*,901)n,iter,fgcnt,timpexp,fxnew,ginf,stepmin,stepmax 901 format(1x,i6,3x,i7,2x,i7,2x,i7,3x,e20.13,1x,e20.13,7x,e20.13, * 1x,e20.13) end if if(ngama .ne. 0) then write(1,905)n,iter,fgcnt,timpexp,fxnew,ginf,ngama,stepmin,stepmax write(*,905)n,iter,fgcnt,timpexp,fxnew,ginf,ngama,stepmin,stepmax 905 format(1x,i6,3x,i7,2x,i7,2x,i7,3x,e20.13,1x,e20.13,1x,i4, * 2x,e20.13,1x,e20.13) end if end if if(igamab) then write(1,904)n,iter,fgcnt,timpexp,fxnew,ginf,stepmin,stepmax write(*,904)n,iter,fgcnt,timpexp,fxnew,ginf,stepmin,stepmax 904 format(1x,i6,3x,i7,2x,i7,2x,i7,3x,e20.13,1x,e20.13,1x,e20.13, * 1x,e20.13) end if if(n.eq.1000) then write(2,903)nexp,n, iter, fgcnt, timpexp 903 format(i2,i6,2x,i6,2x,i6,2x,i6) else write(2,902) n, iter, fgcnt, timpexp 902 format(2x,i6,2x,i6,2x,i6,2x,i6) end if itert = itert + iter fgcntt = fgcntt + fgcnt timp = timp + float(timpexp) cc c write(1,230)(i,x(i),i=1,n) c230 format(5x,i6,4x,e20.13) c-------------------------------------- End DO n end do if(igama) write(1,546) if(igamab) write(1,547) write(1,506) itert, fgcntt, timp/100 506 format(2x,'TOTAL',2x,i8,1x,i8,2x,f7.2,'(seconds)') c------------------------------------------- End DO nexp end do stop end *-------------------------------------------------------------- * *--------------------------------------------------------------- * Subroutine for backtracking. *-------------------------------------------------------------- * subroutine back(n,x,d,fx,ps,alfa,s,nexp, step,fgcnt) * real*8 x(n),d(n) real*8 fx, fy, ps real*8 alfa, s, step integer fgcnt real*8 y(100000), zz(100000) * step=1.d0 nt=1 * 1 continue do i=1,n y(i)=x(i) + step*d(i) end do * call evalfg(n, y, fy, zz, nexp) fgcnt=fgcnt+1 if(fy .le. fx + alfa*step*ps) go to 10 if(fy .gt. fx + alfa*step*ps) then step = step * s nt=nt+1 c write(*,20) nt,step c20 format(5x,'BACK, nt=',i5,4x,'step =',e20.13) go to 1 end if * 10 continue return end *------------------------------- * *--------------------------------------------------------------- * Subroutine for backtracking for the New algorithm *-------------------------------------------------------------- * subroutine backnew(n,x,d,fx,ps,alfa,s,stepini,nexp, step,fgcnt) * real*8 x(n),d(n) real*8 fx, fy, ps real*8 alfa, s, step, stepini integer fgcnt real*8 y(100000), zz(100000) * step=stepini nt=1 * 1 continue do i=1,n y(i)=x(i) + step*d(i) end do * call evalfg(n, y, fy, zz, nexp) fgcnt=fgcnt+1 if(fy .le. fx + alfa*step*ps) go to 10 if(fy .gt. fx + alfa*step*ps) then step = step * s nt=nt+1 go to 1 end if * 10 continue c write(1,50) nt, step c50 format(/,2x,'BACK: nt=',i4,3x,' step=',e20.13,' ====') return end *------------------------------- *-------------------------------------------------------------- c * c *** c ******* c ********* *-------------------------------------------------------------- * Date created : May 30, 1995 * Date last modified : May 30, 1995 * * Subroutine for execution time computation. * *---------------------------------------------------------------- * subroutine exetim(tih,tim,tis,tic,tfh,tfm,tfs,tfc) * integer*4 tih,tim,tis,tic integer*4 tfh,tfm,tfs,tfc * integer*4 ti,tf integer*4 ch,cm,cs data ch,cm,cs/360000,6000,100/ * ti=tih*ch+tim*cm+tis*cs+tic tf=tfh*ch+tfm*cm+tfs*cs+tfc tf=tf-ti tfh=tf/ch tf=tf-tfh*ch tfm=tf/cm tf=tf-tfm*cm tfs=tf/cs tfc=tf-tfs*cs * return end *---------------------------------------------------------------- c * c *** c ******* c ********* ****************************************************************** subroutine inipoint(n,x, nexp) C This subroutine computes the initial point real*8 x(n) go to ( 1, 2, 3, 4, 5, 6, 7, 8, 9,10, * 11,12,13,14,15,16,17,18,19,20, * 21,22,23,24,25,26,27,28,29,30, * 31,32,33,34,35,36,37,38,39,40, * 41,42,43,44,45,46,47,48,49,50, * 51,52,53,54,55,56,57,58,59,60, * 61,62,63,64,65) nexp 1 continue c Freudenstein & Roth i=1 991 x(i) = 0.5d0 x(i+1)= -2.d0 i=i+2 if(i.le.n) go to 991 return 2 continue c Trigonometric do i=1,n x(i) = 0.2d0 end do return 3 continue c Extended Rosenbrock i=1 993 x(i) = -1.2d0 x(i+1)= 1.d0 i=i+2 if(i.le.n) go to 993 return 4 continue c Extended White & Holst i=1 994 x(i) = -1.2d0 x(i+1)= 1.d0 i=i+2 if(i.le.n) go to 994 return 5 continue c Extended Beale i=1 995 x(i) = 1.d0 x(i+1)= 0.8d0 i=i+2 if(i.le.n) go to 995 return 6 continue c Penalty do i=1,n x(i) = float(i) end do return 7 continue c Perturbed Quadratic do i=1,n x(i) = 0.5d0 end do return 8 continue c Raydan 1 do i=1,n x(i) = 1.d0 end do return 9 continue c Raydan 2 do i=1,n x(i) = 1.d0 end do return 10 continue c Diagonal 1 do i=1,n x(i) = 1.d0/float(n) end do return 11 continue c Diagonal 2 do i=1,n x(i) = 1.d0/float(i) end do return 12 continue c Diagonal 3 do i=1,n x(i) = 1.d0 end do return 13 continue c Hager do i=1,n x(i) = 1.d0 end do return 14 continue c Generalized Tridiagonal 1 do i=1,n x(i) = 2.d0 end do return 15 continue c Extended Tridiagonal 1 do i=1,n x(i) = 2.d0 end do return 16 continue c Extended Three Expo Terms do i=1,n x(i) = 0.1d0 end do return 17 continue c Generalized Tridiagonal 2 do i=1,n x(i) = -1.d0 end do return 18 continue c Diagonal 4 do i=1,n x(i) = 1.d0 end do return 19 continue c Diagonal 5 do i=1,n x(i) = 1.1d0 end do return 20 continue c Extended Himmelblau do i=1,n x(i) = 1.d0 end do return 21 continue c Generalized PSC1 i=1 9921 x(i) = 3.d0 x(i+1)= 0.1d0 i=i+2 if(i.le.n) go to 9921 return 22 continue c Extended PSC1 i=1 9922 x(i) = 3.d0 x(i+1)= 0.1d0 i=i+2 if(i.le.n) go to 9922 return 23 continue c Extended Powell i=1 9923 x(i) = 3.d0 x(i+1)= -1.d0 x(i+2)= 0.d0 x(i+3)= 1.d0 i=i+4 if(i.le.n) go to 9923 return 24 continue c Extended BD1 do i=1,n x(i) = 0.1d0 end do return 25 continue c Extended Maratos i=1 9925 x(i) = 1.1d0 x(i+1)= 0.1d0 i=i+2 if(i.le.n) go to 9925 return 26 continue c Extended Cliff i=1 9926 x(i) = 0.d0 x(i+1)= -1.d0 i=i+2 if(i.le.n) go to 9926 return 27 continue c Quadratic Diagonal Perturbed do i=1,n x(i) = 0.5d0 end do return 28 continue c Extended Wood i=1 9928 x(i) = -3.d0 x(i+1)= -1.d0 x(i+2)= -3.d0 x(i+3)= -1.d0 i=i+4 if(i.le.n) go to 9928 return 29 continue c Extended Hiebert do i=1,n x(i) = 0.0001d0 end do return 30 continue c Quadratic QF1 do i=1,n x(i) = 1.d0 end do return 31 continue c Extended Quadratic Penalty QP1 do i=1,n x(i) = 1.d0 end do return 32 continue c Extended Quadratic Penalty QP2 do i=1,n x(i) = 1.d0 end do return 33 continue c Quadratic QF2 do i=1,n x(i) = 0.5d0 end do return 34 continue c Extended EP1 do i=1,n x(i) = 1.5d0 end do return 35 continue c Extended Tridiagonal do i=1,n x(i) = 1.d0 end do return 36 continue c BDQRTIC (CUTE) do i=1,n x(i) = 1.d0 end do return 37 continue c TRIDIA (CUTE) do i=1,n x(i) = 1.d0 end do return 38 continue c ARWHEAD (CUTE) do i=1,n x(i) = 1.d0 end do return 39 continue c NONDIA (CUTE) do i=1,n x(i) = -1.d0 end do return 40 continue c NONDQUAR (CUTE) i=1 9940 x(i) = 1.d0 x(i+1)= -1.d0 i=i+2 if(i.le.n) go to 9940 return 41 continue c DQDRTIC (CUTE) do i=1,n x(i) = 3.d0 end do return 42 continue c EG2 (CUTE) do i=1,n x(i) = 1.d0 end do return 43 continue c DIXMAANA (CUTE) do i=1,n x(i) = 2.d0 end do return 44 continue c DIXMAANB (CUTE) do i=1,n x(i) = 2.d0 end do return 45 continue c DIXMAANC (CUTE) do i=1,n x(i) = 2.d0 end do return 46 continue c DIXMAANE (CUTE) do i=1,n x(i) = 2.d0 end do return 47 continue c Partial Perturbed Quadratic do i=1,n x(i) = 0.5d0 end do return 48 continue c Broyden Tridiagonal do i=1,n x(i) = -1.d0 end do return 49 continue c Almost Perturbed Quadratic do i=1,n x(i) = 0.5d0 end do return 50 continue c Tridiagonal Perturbed Quadratic do i=1,n x(i) = 0.5d0 end do return 51 continue c Generalized Rosenbrock n=100,200,...,1000 i=1 9951 x(i) = -1.2d0 x(i+1)= 1.d0 i=i+2 if(i.le.n) go to 9951 return 52 continue c Generalized White & Holst i=1 9952 x(i) = -1.2d0 x(i+1)= 1.d0 i=i+2 if(i.le.n) go to 9952 return 53 continue c Full Hessian FH1 do i=1,n x(i) = 0.01d0 end do return 54 continue c Full Hessian FH2 do i=1,n x(i) = 0.01d0 end do return 55 continue c ARGLINB (CUTE) do i=1,n x(i) = 1.0d0 end do return 56 continue c CURLY20 do i=1,n x(i) = 0.0001d0/float(n+1) end do return 57 continue c DIAGONAL 6 do i=1,n x(i) = 1.d0 end do return 58 continue c DIXON3DQ 6 do i=1,n x(i) = -1.d0 end do return 59 continue c DIXMAANF (CUTE) do i=1,n x(i) = 2.d0 end do return 60 continue c DIXMAANG (CUTE) do i=1,n x(i) = 2.d0 end do return 61 continue c DIXMAANH (CUTE) do i=1,n x(i) = 2.d0 end do return 62 continue c DIXMAANI (CUTE) do i=1,n x(i) = 2.d0 end do return 63 continue c DIXMAANJ (CUTE) do i=1,n x(i) = 2.d0 end do return 64 continue c DIXMAANK (CUTE) do i=1,n x(i) = 2.d0 end do return 65 continue c DIXMAANL (CUTE) do i=1,n x(i) = 2.d0 end do return 66 continue c ENGVAL1 (CUTE) do i=1,n x(i) = 2.d0 end do return 67 continue c HATFLDA (CUTE) do i=1,n x(i) = 0.1d0 end do return 68 continue c FLETCHCR (CUTE) do i=1,n x(i) = 0.d0 end do return 69 continue c COSINE (CUTE) do i=1,n x(i) = 1.d0 end do return end c------------------------------------------------ End INIPOINT