MITgcm/lsopt/lsline.F

4cee17c1be Patr*0001
                0002       subroutine lsline(
                0003      &     simul
                0004      &     , nn, ifail, lphprint
                0005      &     , ifunc, nfunc
                0006      &     , ff, dotdg
                0007      &     , tmin, tmax, tact, epsx
                0008      &     , dd, gg, xx, xdiff
                0009      &     )
                0010
                0011 c     ==================================================================
                0012 c     SUBROUTINE lsline
                0013 c     ==================================================================
                0014 c
                0015 c     o line search algorithm for determining control vector update;
                0016 c       After computing updated control vector from given gradient,
                0017 c       a forward and adjoint model run are performed (simul.F)
                0018 c       using the updated control vector.
                0019 c       Tests are then applied to see whether solution has improved.
                0020 c
                0021 c     o Reference: J.C. Gilbert & C. Lemarechal
                0022 c                  Some numerical experiments with variable-storage
                0023 c                  quasi-Newton algorithms
                0024 c                  Mathematical Programming 45 (1989), pp. 407-435
                0025 c
                0026 c     o started: ??? not reproducible
                0027 c
                0028 c     o changed: Patrick Heimbach, MIT/EAPS
                0029 c
                0030 c     o Version: 2.0, 24-Feb-2000: Patrick Heimbach, MIT/EAPS
                0031 c                   - severe changes in structure including various
                0032 c                     shifts of variables which are only used in this
                0033 c                     routine
                0034 c                   - number of 3 control flags for error handling
                0035 c                     (indic, moderl, ifail) reduced to 1 (ifail)
                0036 c                     and homogenized with routine lsoptv
                0037 c
                0038 c     o Version: 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS
                0039 c                   - initial computation of tact and
                0040 c                     xdiff = xx + tact*dd
                0041 c                     moved to new routine lsupdxx
                0042 c                     tmin, tmax, tact needed as parameters
                0043 c
                0044 c     ==================================================================
                0045 c     SUBROUTINE lsline
                0046 c     ==================================================================
                0047
f50e6c1777 Patr*0048 #include "blas1.h"
4cee17c1be Patr*0049
                0050       implicit none
                0051
                0052 c----------------------------------
                0053 c declare arguments
                0054 c----------------------------------
                0055       integer nn, ifail, ifunc, nfunc
                0056       double precision ff, dotdg, tmin, tmax, tact, epsx
                0057       double precision xx(nn), dd(nn), gg(nn), xdiff(nn)
                0058
                0059       logical lphprint
                0060
                0061       external simul
                0062
                0063 c----------------------------------
                0064 c declare local variables
                0065 c----------------------------------
                0066
                0067       double precision xpara1, xpara2
                0068       parameter( xpara1 = 0.0001, xpara2=0.9 )
                0069
                0070       double precision    factor
                0071       parameter( factor = 0.2 )
                0072
                0073       double precision    barmax
                0074       parameter( barmax = 0.3 )
                0075       double precision    barmul
                0076       parameter( barmul = 5.0 )
                0077       double precision    barmin
                0078       parameter( barmin = 0.01 )
                0079
                0080       integer i, indic
                0081
                0082       double precision    tg, fg, td, ta
                0083       double precision    fa, fpa, fp
                0084       double precision    fnew, fdiff
                0085       double precision    z, test, barr
                0086       double precision    left, right, told
                0087
ae287463fd Patr*0088       external DDOT
                0089       double precision     DDOT
4cee17c1be Patr*0090
                0091 c----------------------------------
                0092 c check parameters
                0093 c----------------------------------
                0094
                0095       if (  (nn.le.0)
                0096      &     .or. (dotdg.ge.0.0)
                0097      &     .or. (xpara1.le.0.0) .or. (xpara1.ge.0.5)
                0098      &     .or. (xpara2.le.xpara1)  .or. (xpara2.ge.1.0) ) then
                0099          ifail = 9
                0100          go to 999
                0101       endif
                0102
                0103 c----------------------------------
                0104 c initialization
                0105 c----------------------------------
                0106       indic = 0
                0107
                0108       barr   = barmin
                0109       fg     = ff
                0110       fa     = ff
                0111       fpa    = dotdg
                0112
                0113       td     = 0.0
                0114       tg     = 0.0
                0115       ta     = 0.0
                0116
                0117 c=======================================================================
                0118 c begin of simulation iter.
                0119 c=======================================================================
                0120
                0121       do ifunc = 1, nfunc
                0122
                0123          if (lphprint)
                0124      &        print *, 'pathei-lsopt: ', ifunc, ' simul.'
                0125
                0126 c------------------------------------
                0127 c compute cost function and gradient
                0128 c------------------------------------
                0129          call simul( indic, nn, xdiff, fnew, gg )
                0130
ae287463fd Patr*0131          fp = DDOT( nn, dd, 1, gg, 1 )
4cee17c1be Patr*0132          fdiff = fnew - ff
                0133
                0134 c-----------------------------------------
                0135 c apply 1st Wolfe test
                0136 c-----------------------------------------
                0137          if (fdiff .gt. tact*xpara1*dotdg) then
                0138             td     = tact
                0139             ifail  = 0
                0140             go to 500
                0141          endif
                0142
                0143 c-----------------------------------------
                0144 c 1st Wolfe test 1 ok, apply 2nd Wolf test
                0145 c-----------------------------------------
                0146          if (fp .gt. xpara2*dotdg) then
                0147             ifail = 0
                0148             go to 320
                0149          endif
                0150
                0151          if (ifail.eq.0) go to 350
                0152
                0153 c-----------------------------------------
                0154 c 2nd Wolfe test 2 ok, donc pas serieux, on sort
                0155 c-----------------------------------------
                0156
                0157  320     continue
                0158
                0159          ff = fnew
                0160          do i = 1, nn
                0161             xx(i) = xdiff(i)
                0162          end do
                0163 cph(
                0164          if (lphprint)
                0165      &        print *, 'pathei-lsopt: no inter-/extrapolation in lsline'
                0166 cph)
                0167          go to 999
                0168
                0169
                0170 c-----------------------------------------
                0171 c extrapolation
                0172 c-----------------------------------------
                0173
                0174  350     continue
                0175
                0176          tg  = tact
                0177          fg  = fnew
                0178          if (td .ne. 0.0) go to 500
                0179
                0180          told    = tact
                0181          left = (1.0+barmin)*tact
                0182          right = 10.0*tact
                0183          call cubic( tact, fnew, fp, ta, fa, fpa, left, right )
                0184          ta     = told
                0185          if (tact.ge.tmax) then
                0186             ifail = 7
                0187             tact     = tmax
                0188          endif
                0189
                0190          if (lphprint)
                0191      &        print *, 'pathei-lsopt: extrapolation: ',
                0192      &        'td, tg, tact, ifail = ', td, tg, tact, ifail
                0193
                0194          go to 900
                0195
                0196 c-----------------------------------------
                0197 c interpolation
                0198 c-----------------------------------------
                0199
                0200  500     continue
                0201
                0202          test   = barr*(td-tg)
                0203          left = tg+test
                0204          right = td-test
                0205          told    = tact
                0206          call cubic( tact, fnew, fp, ta, fa, fpa, left, right )
                0207          ta = told
                0208          if (tact.gt.left .and. tact.lt.right) then
                0209             barr = dmax1( barmin, barr/barmul )
                0210          else
                0211             barr = dmin1( barmul*barr, barmax )
                0212          endif
                0213
                0214          if (lphprint)
                0215      &        print *, 'pathei-lsopt: interpolation: ',
                0216      &        'td, tg, tact, ifail = ', td, tg, tact, ifail
                0217
                0218 c-----------------------------------------
                0219 c end of iteration loop
                0220 c     - tact peut etre bloque sur tmax
                0221 c       (venant de lextrapolation avec ifail=7)
                0222 c-----------------------------------------
                0223
                0224  900     continue
                0225
                0226          fa  = fnew
                0227          fpa = fp
                0228
                0229 c-----------------------------------------
                0230 c --- faut-il continuer ?
                0231 c-----------------------------------------
                0232          if (td .eq. 0.0)     go to 950
                0233          if (td-tg .lt. tmin) go to 920
                0234
                0235 c-----------------------------------------
                0236 c limit due to machine precision
                0237 c-----------------------------------------
                0238          do i = 1, nn
                0239             z = xx(i) + tact*dd(i)
                0240             if ((z.ne.xx(i)) .and. (z.ne.xdiff(i))) go to 950
                0241          end do
                0242
                0243 c-----------------------------------------
                0244 c arret sur dxmin ou de secours
                0245 c-----------------------------------------
                0246  920     continue
                0247          ifail = 8
                0248
                0249 c-----------------------------------------
                0250 c     si tg=0, xx = xx_depart,
                0251 c     sinon on prend xx=x_left qui fait decroitre f
                0252 c-----------------------------------------
                0253          if (tg .ne. 0.0) then
                0254             ff = fg
                0255             do i = 1, nn
                0256                xx(i) = xx(i) + tg*dd(i)
                0257             end do
                0258          endif
                0259
                0260          go to 999
                0261
                0262 c-----------------------------------------
                0263 c update vector for new simulation
                0264 c-----------------------------------------
                0265  950     continue
                0266
                0267          do i = 1, nn
                0268             xdiff(i) = xx(i) + tact*dd(i)
                0269          end do
                0270
                0271 c=======================================================================
                0272 c end of simulation iter.
                0273 c=======================================================================
                0274       end do
                0275
                0276 c-----------------------------------------
                0277 c too many function evaluations
                0278 c-----------------------------------------
                0279       ifail = 6
                0280       ff    = fg
                0281       do i = 1, nn
                0282          xx(i) = xx(i) + tg*dd(i)
                0283       end do
                0284
                0285 c-----------------------------------------
                0286 c end of routine
                0287 c-----------------------------------------
                0288   999 continue
                0289
                0290       return
                0291
                0292       end