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
