C Spectral Conjugate Gradient method (SCG) C C Subroutine scg implements the Spectral Conjugate C Gradient method to find local minimizers of a given C function. The method is described in: C C E. G. Birgin and J. M. Martinez, "A spectral C conjugate gradient method for unconstrained C optimization", Applied Mathematics and C Optimization 43, pp. 117-128, 2001. C C Subroutines inipoint and evalfg must be supplied by C the user to set the initail guess and to evaluate the C function f and its gradient vector g, simultaneously, C respectively. C C This version: August 31, 2004. program scgma C LOCAL SCALARS integer betatype,n,iter,fgcnt,lscnt,maxiter double precision eps,f,gnorm C LOCAL ARRAYS double precision x(10000) C Dimension n = 100 C Initial guess call inipoint(n,x) C Solver parameters (betatype = 1, 2 and 3 means Perry, C Polak-Riviere and Fletcher-Reeves, respectively). eps = 1.0d-06 maxiter = 15000 betatype = 1 C Call the solver call scg(n,x,eps,maxiter,betatype,f,gnorm,iter,fgcnt,lscnt) C Write statistics write(*,100) iter,fgcnt,lscnt,f,gnorm stop C NON-EXECUTABLE STATEMENTS 100 format(/'NUMBER OF ITERATIONS = ',8X,I6, + /'NUMBER OF FUNCTION/GRADIENT EVALUATIONS = ',8X,I6, + /'TOTAL NUMBER OF LINE SEARCH ITERATIONS = ',8X,I6, + /'OBJECTIVE FUNCTION VALUE = ',1PD14.7, + /'GRADIENT EUCLIDIAN NORM = ',1PD14.7) end c ****************************************************************** c ****************************************************************** subroutine inipoint(n,x) C This subroutine computes the initial point C INPUT PARAMETERS C C integer n : problem dimension C C OUTPUT PARAMETERS C C double precision x(n) : initial point C SCALAR ARGUMENTS integer n C ARRAY ARGUMENTS double precision x(n) C LOCAL SCALARS integer i do i = 1,n x(i) = 1.0d0 end do end c ****************************************************************** c ****************************************************************** subroutine evalfg(n,x,f,g) C This subroutine computes the objective function and its C gradient simultaneously. C INPUT PARAMETERS C C integer n : problem dimension C C double precision x(n) : point at which the objective C funtion and its gradient will be computed C C OUTPUT PARAMETERS C C double precision f : objective function value at x C C double precision g(n) : gradient at x C SCALAR ARGUMENTS integer n double precision f C ARRAY ARGUMENTS double precision x(n),g(n) C LOCAL SCALARS integer i f = 0.0d0 do i = 1,n f = f + ( 1.0d0 + mod(i,10) ) * x(i) ** 2 g(i) = ( 1.0d0 + mod(i,10) ) * 2.0d0 * x(i) end do end c ****************************************************************** c ****************************************************************** subroutine scg(n,x,eps,maxiter,betatype,f,gnorm,iter,fgcnt,lscnt) C Spectral Conjugate Gradient C SCALAR ARGUMENTS integer betatype,n,iter,fgcnt,lscnt,maxiter double precision eps,f,gnorm C ARRAY ARGUMENTS double precision x(n) C LOCAL SCALARS integer i,lsflag double precision fnew,gtd,alpha,theta,thetaprev,beta,gtg, + gtgprev,dnorm,sts,sty,dnormprev,a,b C LOCAL ARRAYS double precision xnew(10000),g(10000),gnew(10000),d(10000), + y(10000),s(10000) C Initialization write(*,200) iter = 0 fgcnt = 0 lscnt = 0 theta = 1.0d0 call evalfg(n,x,f,g) fgcnt = fgcnt + 1 gtg = 0.0d0 do i = 1,n d(i) = - g(i) gtg = gtg + g(i) ** 2 end do gnorm = sqrt( gtg ) write(*,300) iter,f,gnorm gtd = -gtg dnorm = gnorm if ( gnorm .ne. 0.0d0 ) then alpha = 1.0d0 / gnorm end if C Main loop 10 if ( gnorm .gt. eps * max( 1.0d0, abs( f ) ) .and. + iter .lt. maxiter ) then iter = iter + 1 call LineSearch(n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew,fgcnt, + lscnt,lsflag) if ( lsflag .ne. 0 ) then write(*,*) 'Warning! Abnormal Line Search termination.' end if f = fnew gtgprev = gtg sts = 0.0d0 sty = 0.0d0 gtg = 0.0d0 do i = 1,n s(i) = xnew(i) - x(i) y(i) = gnew(i) - g(i) sts = sts + s(i) ** 2 sty = sty + s(i) * y(i) x(i) = xnew(i) g(i) = gnew(i) gtg = gtg + g(i) ** 2 end do gnorm= sqrt( gtg ) write(*,300) iter,f,gnorm thetaprev = theta if ( sty .le. 0.0d0 ) then theta = 1.0d+10 else theta = min( 1.0d+10, max( 1.0d-10, sts / sty ) ) end if if ( betatype .eq. 1 ) then C Perry a = 0.0d0 b = 0.0d0 do i = 1,n a = a + ( theta * y(i) - s(i) ) * g(i) b = b + s(i) * y(i) end do else if ( betatype .eq. 2 ) then C Polak-Ribiere a = 0.0d0 do i = 1,n a = a + g(i) * y(i) end do a = theta * a b = alpha * thetaprev * gtgprev else if ( betatype .eq. 3 ) then C Fletcher-Reeves a = theta * gtg b = alpha * thetaprev * gtgprev else write(*,*) 'Warning! Undefined beta type.' b = 0.0d0 end if if ( b .ne. 0.0d0 ) then beta = a / b else beta = 0.0d0 end if dnormprev = dnorm dnorm = 0.0d0 gtd = 0.0d0 do i = 1,n d(i) = - theta * g(i) + beta * s(i) dnorm = dnorm + d(i) ** 2 gtd = gtd + g(i) * d(i) end do dnorm = sqrt( dnorm ) if ( gtd .gt. -1.0d-03 * gnorm * dnorm ) then do i = 1,n d(i) = - theta * g(i) end do dnorm = theta * gnorm gtd = - theta * gtg end if alpha = alpha * dnormprev / dnorm goto 10 end if C End of main loop return C NON-EXECUTABLE STATEMENTS 200 format('SPECTRAL CONJUGATE GRADIENT METHOD') 300 format('ITER = ',I6,' F = ',1PD14.7,' GNORM = ',1PD14.7) end c ****************************************************************** c ****************************************************************** subroutine LineSearch (n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew, + fgcnt,lscnt,lsflag) C This is the one-dimensional line search used in CONMIN C SCALAR ARGUMENTS integer n,fgcnt,lscnt,lsflag double precision f,gtd,dnorm,alpha,fnew C ARRAY ARGUMENTS double precision x(n),d(n),xnew(n),gnew(n) C LOCAL SCALARS integer i,lsiter double precision alphap,alphatemp,fp,dp,gtdnew,a,b lsflag = 0 alphap = 0.0d0 fp = f dp = gtd do i = 1,n xnew(i) = x(i) + alpha * d(i) end do call evalfg(n,xnew,fnew,gnew) fgcnt = fgcnt + 1 gtdnew = 0.0d0 do i = 1,n gtdnew = gtdnew + gnew(i) * d(i) end do lsiter = 0 10 if ( alpha * dnorm .gt. 1.0d-30 .and. lsiter .lt. 50 .and. + .not. ( gtdnew .eq. 0.0d0 .and. fnew .lt. f ) .and. + ( ( fnew .gt. f + 1.0d-04 * alpha * gtd .or. + abs( gtdnew / gtd ) .gt. 0.9d0 ) .or. ( lsiter .eq. 0 .and. + abs( gtdnew / gtd ) .gt. 0.5d0 ) ) ) then 20 if ( alpha * dnorm .gt. 1.0d-30 .and. fnew .gt. f .and. + gtdnew .lt. 0.0d0 ) then alpha = alpha / 3.0d0 do i = 1,n xnew(i) = x(i) + alpha * d(i) end do call evalfg(n,xnew,fnew,gnew) fgcnt = fgcnt + 1 gtdnew = 0.0d0 do i = 1,n gtdnew = gtdnew + gnew(i) * d(i) end do alphap = 0.0d0 fp = f dp = gtd goto 20 end if a = dp + gtdnew - 3.0d0 * ( fp - fnew ) / ( alphap - alpha ) b = a ** 2 - dp * gtdnew if ( b .gt. 0.0d0 ) then b = sqrt( b ) else b = 0.0d0 end if alphatemp = alpha - ( alpha - alphap ) * ( gtdnew + b - a ) / + ( gtdnew - dp + 2.0d0 * b ) if ( gtdnew / dp .le. 0.0d0 ) then if ( 0.99d0 * max( alpha, alphap ) .lt. alphatemp .or. + alphatemp .lt. 1.01d0 * min( alpha, alphap ) ) then alphatemp = ( alpha + alphap ) / 2.0d0 end if else if ( gtdnew .lt. 0.0d0 .and. + alphatemp .lt. 1.01d0 * max( alpha, alphap ) ) then alphatemp = 2.0d0 * max( alpha, alphap ) end if if ( ( gtdnew .gt. 0.0d0 .and. + alphatemp .gt. 0.99d0 * min( alpha, alphap ) ) .or. + alphatemp .lt. 0.0d0 ) then alphatemp = min( alpha, alphap ) / 2.0d0 end if end if alphap = alpha fp = fnew dp = gtdnew alpha = alphatemp do i = 1,n xnew(i) = x(i) + alpha * d(i) end do call evalfg(n,xnew,fnew,gnew) fgcnt = fgcnt + 1 gtdnew = 0.0d0 do i = 1,n gtdnew = gtdnew + gnew(i) * d(i) end do lsiter = lsiter + 1 goto 10 end if if ( lsiter .ge. 50 ) then lsflag = 1 end if if ( lsiter .ne. 0 ) then lscnt = lscnt + 1 end if return end