subroutine opt(N,x,f,g,fevals,iter,M,retcode)
      implicit none
      double precision x(N), f, g(N)
      integer N,fevals, iter, M, retcode
c
c     ----------------------------------------------------------------------
c     Subroutine: opt
c
c     Author: Tamara Gibson (Kolda), Fall 1995.
c     Copyright (c) 1996.
c     University of Maryland at College Park.
c     All rights reserved.
c
c     Implements L-BFGS with variations as described in:
c     Tamara Gibson (Kolda), Dianne P. O'Leary, and Larry Nazareth, 
c     "BFGS with Update Skipping and Varying Memory", 
c     Technical Report CS-TR-3663, Department of Computer Science, 
c     University of Maryland at College Park, 1996. 
c     Also listed as Technical Report UMIACS-TR-96-49.
c
c     Parameters: (i - on input, o - on output)
c
c     N:         INTEGER, (i) dimension of problem
c     x:         DOUBLE PRECISION N-vector, (i) starting point 
c                                           (o) minimizer of f
c     f:         DOUBLE PRECISION, (i/o) value of f at x 
c     g:         DOUBLE PRECISION N-vector, (i/o) gradient vector of f at x
c     fevals:    INTEGER, (o) number of function evaluations
c     iter:      INTEGER, (o) number of bfgs iterations
c     M:         INTEGER, (o) final value of M
c     retcode:   INTEGER, (o) 
c                  -1 - no Progress (line search failed on 1st iteration)
c                   0 - |g| < eps
c                   1 - failed line search
c                   2 - misc error
c                   3 - reached maximum number of iterations (MAXIT)
c
c     ----------------------------------------------------------------------
c
c     Possible Compile Line Options
c
c     VARYMG 
c      0 -> M varies but never gets smaller
c      1 -> M varies 
c
c     USEGOVERX (used in conjunction with VARYMG)
c
c     SHRINKM
c      0 -> If a step larger than 1 is taken, k is set to 1
c
c     BACKUP
c      0 -> Backup if the iteration number is odd
c      1 -> Backup if the iteration number is even
c      2 -> Backup on the NEXT iteration if the current step is
c           larger than 1.0
c      3 -> Backup on the NEXT iteration if the current norm of g is 
c           larger than the previous norm of g
c     
c     BUOPT (used in conjunction with BACKUP to specify that there
c            should not be two backups in a row)
c
c     MERGE
c      0 -> Merge the the 2nd and 3rd most recent (s,y) pairs every
c           2nd step
c      1 -> Merge the 2nd and 3rd most recent (s,y) pairs if the
c           corresponding step lengths were 1 AND neither pair
c           is the result of a previous merge
c      2 -> Merge the 2nd and 3rd most recent (s,y) pairs if the
c           corresponding step lengths were >1 AND neither pair
c           is the result of a previous merge
c      3 -> Merge the 2nd and 3rd most recent (s,y) pairs if either
c           of the corresponding step lengths was >1 AND neither pair
c           is the result of a previous merge
c
c     SKIP
c      0 -> Skip if the iteration number is odd
c      1 -> Skip if the iteration number is even
c      2 -> Skip if the current norm of g is larger than the previous
c
c     ----------------------------------------------------------------------
c

c
c EPS is the covergence tolerance.
c
#ifndef EPS
#define EPS (1.0d-5)
#endif

c
c The following are line search parameters.
c See the documentation in the line search routine to
c understand how they are used.
c
#ifndef FTOL
#define FTOL (1.0d-4)
#endif
#ifndef GTOL
#define GTOL (0.9d0)
#endif
#ifndef XTOL
#define XTOL (1.0d-15)
#endif
#ifndef STPMIN
#define STPMIN (1.0d-15)
#endif
#ifndef STPMAX
#define STPMAX (1.0d+15)
#endif
#ifndef MAXFEV
#define MAXFEV 20
#endif

c
c MVAL is the maximum allowable value of M
c
#ifndef MVAL
#define MVAL 20
#endif

c
c MAXN is the maximum allowable value of N (problem dimension)
c
#ifndef MAXN
#define MAXN 10000
#endif

c
c MAXIT is the maximum allowable number of iterations.
c
#ifndef MAXIT
#define MAXIT 3000
#endif

c
c     ----------------------------------------------------------------------
c     Variable Declarations  
c     ----------------------------------------------------------------------
c
      double precision oldg(MAXN), d(MAXN), s(MAXN), y(MAXN), stp
      double precision bigs(MAXN,MVAL), bigy(MAXN,MVAL)
      double precision bigstg(MVAL), bigytg(MVAL), bigsty(MVAL)
      double precision bigstoldg(MVAL), bigytoldg(MVAL)
      double precision invr(MVAL,MVAL), yty(MVAL,MVAL), bigd(MVAL)
      double precision gamma
      double precision p1(MVAL),p2(MVAL)
      double precision vec1(MAXN)
      double precision rho, delta
      double precision one, zero
      double precision normg, normx, testval
      double precision wa(MAXN)
      integer info
      integer lsfevals
      integer i,j,k,l

#ifdef MERGE
      logical mergeprev
      double precision smerge(MAXN), ymerge(MAXN)
      double precision deltamerge, rhomerge
      double precision oldoldstp, oldstp
#endif
#ifdef SKIP
      double precision normoldg
#endif
#ifdef BACKUP
      integer backupnext
      double precision normoldg
#endif
#ifdef VARYMG
      integer prevM 
      double precision slope,lgnormg0,lgnormx0,testval0
#endif
#ifdef SHRINKM
      double precision snormoldg
#endif

c
c     ----------------------------------------------------------------------
c     Function Declarations
c     ----------------------------------------------------------------------
c
      double precision DDOT
c
c     ----------------------------------------------------------------------
c     Initialization
c     ----------------------------------------------------------------------
c
      retcode = 2

      if (MVAL .gt. MAXN) then
         write(6,*) 'ERROR: MAXN must be at least as big as MVAL.'
         return
      endif
      one = 1.0d0
      zero = 0.0d0
      fevals = 1

#ifdef VARYMG
      M = 1
      lgnormg0 = log10(sqrt(ddot(N,g,1,g,1)))

#ifdef WAIT
      lgnormg0 = 0.0d0
#endif 

#ifdef USEGOVERX
      lgnormx0 = sqrt(ddot(N,x,1,x,1))
      if (lgnormx0 .lt. 1) then
         lgnormx0 = 1
      else
         lgnormx0 = log10(lgnormx0)
      endif
      testval0 = lgnormg0/lgnormx0
      slope = (MVAL-1.0)/(log10(EPS*100) - testval0)
#else
      slope = (MVAL-1.0)/(log10(EPS*100) - lgnormg0)
#endif

#else
      M = MVAL
#endif

#ifdef BACKUP
      backupnext = 0
#endif
#ifdef MERGE
      mergeprev = .false.
      oldstp = zero
      oldoldstp = zero
#endif

c
c     ----------------------------------------------------------------------
c     Main Loop
c     ----------------------------------------------------------------------
c
      retcode = 3
      do iter = 1,MAXIT

#ifdef DEBUG
         write(6,*) ' '
         write(6,*) 'ITERATION # ',iter
         write(6,*) 'Value of k: ',k
         write(6,1500) 'g=',(g(i),i=1,N)
         write(6,1500) 's=',(s(i),i=1,N)
         write(6,1500) 'y=',(y(i),i=1,N)
 1500    format(A3,' ',11E10.2)
#endif

c
c        -------------------------------------------------------------------
c        Compute d = -Hg
c        -------------------------------------------------------------------
c
         if (iter .eq. 1) then
            do i = 1,N 
               d(i) = -g(i)
            enddo
            k = 1
         else
c
c           ---------------------------------------
c           The following are various modifications 
c           of the standard LM-BFGS algorithm.
c           ---------------------------------------
c

#ifdef MERGE

#if MERGE == 1
            if (mergeprev) then
               mergeprev = .false.
            elseif (k .gt. 2) then
               if ((oldoldstp .eq. 1.0d0) .and. (oldstp .eq. 1.0d0)) then
                  mergeprev = .true.
               endif
            endif
#elif MERGE == 2
            if (mergeprev) then
               mergeprev = .false.
            elseif (k .gt. 2) then
               if ((oldoldstp .gt. 1.0d0) .and. (oldstp .gt. 1.0d0)) then
                  write(6,*) oldoldstp, oldstp
                  mergeprev = .true.
               endif
            endif
#elif MERGE == 3
            if (mergeprev) then
               mergeprev = .false.
            elseif (k .gt. 2) then
               if ((oldoldstp .gt. 1.0d0) .or. (oldstp .gt. 1.0d0)) then
                  mergeprev = .true.
               endif
            endif
#elif MERGE == 0
            if (mergeprev) then
               mergeprev = .false.
            elseif (k .gt. 2) then
               mergeprev = .true.
            endif
#endif

            if (mergeprev) then
               k = k-1

c     Merge s_k and s_(k-1) into s_(k-1).  
c     Do the same for y_k and y_(k-1).

               do i = 1,N
                  smerge(i) = bigS(i,k) + bigS(i,k-1)
                  bigS(i,k-1) = smerge(i)
                  ymerge(i) = bigY(i,k) + bigY(i,k-1)
                  bigY(i,k-1) = ymerge(i)
               enddo

c     Compute deltamerge and rhomerge

               deltamerge = DDOT(N,ymerge,1,smerge,1)
               rhomerge = 1.0d0/deltamerge

c     Compute bigstoldg and bigytoldg

               bigstg(k-1) = DDOT(N,smerge,1,oldg,1)
               bigytg(k-1) = DDOT(N,ymerge,1,oldg,1)

c     Compute invR and YtY

               if (k .gt. 2) then
                  call DGEMV('T',N,k-2,one,bigs,MAXN,ymerge,1,zero,bigsty,1)
                  do i = 1,k-2
                     vec1(i) = -rhomerge * bigsty(i)
                  enddo
                  call DTRMV('U','N','N',k-2,invr,MVAL,vec1,1)
                  do i = 1,k-2
                     invr(i,k-1) = vec1(i)
                  enddo
                  do i = 1,k-2
                     yty(i,k-1) = yty(i,k-1) + yty(i,k)
                  enddo
               endif

               invr(k-1,k-1) = rhomerge
               bigd(k-1) = deltamerge
               yty(k-1,k-1) = DDOT(N,ymerge,1,ymerge,1)

            endif

#endif
#if BACKUP == 0
            if ((iter - (iter/2)*2) .eq. 1) then
               k = k-1
               if (k .lt. 1) k = 1
            endif
#elif BACKUP == 1
            if ((iter - (iter/2)*2) .eq. 0) then
               k = k-1
               if (k .lt. 1) k = 1
            endif
#elif BACKUP == 2
            if (backupnext .eq. 1) then
               k = k - 1
            endif
#ifdef BUOPT
            if (backupnext .eq. 1) then
               backupnext = 0
            elseif (stp .gt. 1.0d0) then
               backupnext = 1
            endif
#else
            if (stp .gt. 1.0d0) then
               backupnext = 1
            else
               backupnext = 0
            endif
#endif
#elif BACKUP == 3
            if (backupnext .eq. 1) then
               k = k - 1
            endif
#ifdef BUOPT
            if (backupnext .eq. 1) then
               backupnext = 0
            elseif (normoldg .lt. normg) then
               backupnext = 1
            endif
#else
            if (normoldg .lt. normg) then
               backupnext = 1
            else
               backupnext = 0
            endif
#endif
#endif


#ifdef SKIP

#if SKIP == 2
            if (normoldg .lt. normg) then
#elif SKIP == 0
            if ((iter - (iter/2)*2) .eq. 1) then
#elif SKIP == 1
            if ((iter - (iter/2)*2) .eq. 0) then
#endif
               k = k - 1
               if (k .eq. 0) then
                  do i = 1,N
                     d(i) = -g(i)
                  enddo
                  k = 1
                  goto 10
               else
                  call DGEMV('T',N,k,one,bigS,MAXN,g,1,zero,bigStg,1)
                  call DGEMV('T',N,k,one,bigY,MAXN,g,1,zero,bigYtg,1)
                  goto 6
               endif
            endif
#endif

#ifdef SHRINKM
#if SHRINKM == 0
            if (stp .gt. 1.0d0) then
#elif SHRINKM == 1
            if (normg .gt. snormoldg) then
#endif
               k = 1
            endif
#endif

#ifdef VARYMG

            prevM = M

#ifdef USEGOVERX
            M = (slope*(log10(testval)-testval0))/1 + 1
#else
            M = (slope*(log10(normg)-lgnormg0))/1 + 1
#endif

#if VARYMG == 0
            if (M .lt. prevM) M = prevM
#else
            if (M .lt. 1) M = 1
#endif

            if (M .gt. MVAL) M = MVAL

#endif

c
c           -----------------
c           End Modifications
c           -----------------
c


c
c     --- Step 0 ---
c
            if (k .eq. 1) then
               bigstoldg(1) = DDOT(N,s,1,oldg,1)
               bigytoldg(1) = DDOT(N,y,1,oldg,1)
            else
               do i = 1,k-1
                  bigstoldg(i) = bigstg(i)
                  bigytoldg(i) = bigytg(i)                  
               enddo
            endif

 1          continue

            if (k .gt. M) then 
               l = k-M
               do i = 1,M-1
                  bigstoldg(i) = bigstoldg(i+l)
                  bigytoldg(i) = bigytoldg(i+l)
                  bigd(i) = bigd(i+l)
                  do j = 1,N
                     bigs(j,i) = bigs(j,i+l)
                     bigy(j,i) = bigy(j,i+l)
                  enddo
                  do j = 1,M-1
                     invr(i,j) = invr(i+l,j+l)
                     yty(i,j) = yty(i+l,j+l)
                  enddo
               enddo
               k = M
            endif
c
c           --- Step 1 ---
c
            do i=1,N
               bigS(i,k) = s(i)
               bigY(i,k) = y(i)
            enddo
c
c           --- Step 2 ---
c
            call DGEMV('T',N,k,one,bigS,MAXN,g,1,zero,bigStg,1)
            call DGEMV('T',N,k,one,bigY,MAXN,g,1,zero,bigYtg,1)
c
c           --- Step 3 ---
c           --- Step 4 ---
c
            if (k .gt. 1) then
               do i = 1,k-1
                  bigSty(i) = bigStg(i) - bigStoldg(i)
               enddo

               do i = 1,k-1
                  vec1(i) = -rho*bigSty(i)
               enddo
               call DTRMV('U','N','N',k-1,invr,MVAL,vec1,1)
               do i = 1,k-1
                  invr(i,k) = vec1(i)
               enddo
               do i = 1,k-1
                  vec1(i) = bigytg(i) - bigytoldg(i)
               enddo
               do i = 1,k-1
                  yty(i,k) = vec1(i)
               enddo
            endif
            invr(k,k) = rho
            bigd(k) = delta
            yty(k,k) = DDOT(N,y,1,y,1)
c
c           --- Step 5 ---
c
            gamma = delta/yty(k,k)
c
c           --- Step 6 ---
c
 6          continue
            do i = 1,k
               p2(i) = -bigstg(i)
            enddo
            call DTRMV('U','N','N',k,invr,MVAL,p2,1) 
            call DSYMV('U',k,-gamma,yty,MVAL,p2,1,zero,p1,1)
            do i = 1,k
               p1(i) = p1(i) - bigd(i)*p2(i)
            enddo
            call DTRMV('U','T','N',k,invr,MVAL,p1,1)
            do i = 1,k
               vec1(i) = gamma*bigytg(i)
            enddo
            call DTRMV('U','T','N',k,invr,MVAL,vec1,1)
            call DAXPY(k,-one,vec1,1,p1,1)
c
c           --- Step 7 ---
c
            call DGEMV('N',N,k,-gamma,bigY,MAXN,p2,1,zero,d,1)
            call DGEMV('N',N,k,-one,bigS,MAXN,p1,1,zero,vec1,1)
            call DAXPY(N,one,vec1,1,d,1)
            call DAXPY(N,-gamma,g,1,d,1)
            k = k + 1
         endif
 10      continue
c
c        -------------------------------------------------------------------
c        Perform Linesearch (which computes x & g)
c        -------------------------------------------------------------------
c
         do i = 1,N
            oldg(i) = g(i)
         enddo

c
c     Do Line Search
c

#if MERGE > 0
         oldoldstp = oldstp
         oldstp = stp
#endif

         if (iter .eq. 1) then
            normg = sqrt(DDOT(N,g,1,g,1))
            stp = 1.0d0/normg
         else

            stp = 1.0d0
         endif

         info = 0

         call cvsrch(N,x,f,g,d,stp,FTOL,GTOL,XTOL,
     .        STPMIN,STPMAX,MAXFEV,info,lsfevals,wa)
         fevals = fevals + lsfevals

#ifdef DUMP
c     Dump x to output
         write(6,999) (x(i), i=1,N)
 999     format(100E18.10)
#endif

c
c        Check exit conditions of line search
c
         if ((info .eq. 0) .or. (info .eq. 6)) then

#ifdef INFO
            if (info .eq. 0) then
               write(6,*) 'ERROR: Improper line search parameters'
            endif
#endif
#ifdef INFO
            write(6,1000) iter,f,zero,stp,info,lsfevals
#endif

            if (iter .eq. 1) then
               retcode = -1
            else
               retcode = 1
            endif

            return

         endif

#ifdef WARN 
         if (info .eq. 2) then
            write(6,*) 'WARNING: relative width < XTOL.' 
         elseif (info .eq. 3) then
            write(6,*) 'WARNING: Too many f-evals.'
         elseif (info .eq. 4) then
            write(6,*) 'WARNING: stp = stpmin'
         elseif (info .eq. 5) then
            write(6,*) 'WARNING: stp = stpmax'
         endif
#endif

         if (stp .eq. 0) then
            write(6,*) 'ERROR: Stp = 0.'
            retcode = 2
            return
         endif

c
c        -------------------------------------------------------------------
c        Compute s 
c        -------------------------------------------------------------------
c
         do i = 1,N
            s(i) = stp * d(i)
         enddo
c
c        -------------------------------------------------------------------
c        Check for Convergence
c        -------------------------------------------------------------------
c
#ifdef SKIP
         normoldg = normg
#endif
#ifdef SHRINKM
         snormoldg = normg
#endif
#ifdef BACKUP
         normoldg = normg
#endif
         normx = sqrt(DDOT(N,x,1,x,1))
         if (normx .lt. one) normx=1
         normg = sqrt(DDOT(N,g,1,g,1))
         testval = normg/normx

#ifdef INFO
            write(6,1000) iter,f,normg,stp,info,lsfevals
#endif

         if (testval .lt. EPS) then
            retcode = 0
            return
         endif
c
c        -------------------------------------------------------------------
c        Compute y
c        -------------------------------------------------------------------
c
         do i = 1,N
            y(i) = g(i) - oldg(i)
         enddo
c
c        -------------------------------------------------------------------
c        Compute delta and rho
c        -------------------------------------------------------------------
c
         delta = DDOT(N,y,1,s,1)

c
c        Check that delta is positive.
c

         if (delta .lt. 0) then
            write(6,*) 'ERROR: delta < 0'
         endif

#ifdef WARN
         if (delta .lt. XTOL) then
            write(6,*) 'WARNING: delta too small'
         endif
#endif
         rho = 1.0d0/delta

      enddo
c
c     ----------------------------------------------------------------------
c     Format Statements
c     ----------------------------------------------------------------------
c
 1000 format('Iter:',I4,' f= ',E10.4,' |g|= ',E10.4,' stp= ',E8.2,
     .     ' retcode:',I2,' fevals:',I2)




c
c     ----------------------------------------------------------------------
c     End of Subroutine
c     ----------------------------------------------------------------------
c
      iter = iter - 1
      end