MITgcm/lsopt/lsupdxx.F

4cee17c1be Patr*0001
                0002       subroutine lsupdxx(
                0003      &     nn, ifail, lphprint
                0004      &     , jmin, jmax, nupdate
                0005      &     , ff, fmin, fold, gnorm0, dotdg
                0006      &     , gg, dd, xx, xdiff
                0007      &     , tmin, tmax, tact, epsx
                0008      &     )
                0009
                0010 c     ==================================================================
                0011 c     SUBROUTINE lsupdxx
                0012 c     ==================================================================
                0013 c
                0014 c     o conceived for variable online/offline version
                0015 c       computes - new descent direction dd based on latest
                0016 c                  available gradient
                0017 c                - new tact based on new dd
                0018 c                - new control vector xx needed for offline run
                0019 c
                0020 c     o started: Patrick Heimbach, MIT/EAPS
                0021 c                29-Feb-2000:
                0022 c
                0023 c     o Version 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS
                0024 c
                0025 c     ==================================================================
                0026 c     SUBROUTINE lsupdxx
                0027 c     ==================================================================
                0028 c
                0029
f50e6c1777 Patr*0030 #include "blas1.h"
4cee17c1be Patr*0031
                0032       implicit none
                0033
                0034 c-----------------------------------------
                0035 c declare arguments
                0036 c-----------------------------------------
                0037       integer nn, jmin, jmax, nupdate, ifail
                0038       double precision    ff, fmin, fold, gnorm0, dotdg
                0039       double precision    gg(nn), dd(nn), xx(nn), xdiff(nn)
                0040       double precision    tmin, tmax, tact, epsx
                0041       logical lphprint
                0042
                0043 c-----------------------------------------
                0044 C declare local variables
                0045 c-----------------------------------------
                0046       integer i
                0047       double precision    fdiff, preco
                0048
ae287463fd Patr*0049       double precision     DDOT
                0050       external DDOT
4cee17c1be Patr*0051
                0052 c     ==================================================================
                0053
                0054 c-----------------------------------------
                0055 c use Fletchers scaling
                0056 c and initialize diagional to 1.
                0057 c-----------------------------------------
                0058 c
                0059       if ( ( jmax .eq. 0 ) .or. (nupdate .eq. 0 ) ) then
                0060
                0061          if (jmax .eq. 0) then
                0062             fold = fmin
                0063             if (lphprint)
                0064      &           print *, 'pathei-lsopt: using fold = fmin = ', fmin
                0065          end if
                0066          fdiff = fold - ff
                0067          if (jmax .eq. 0) fdiff = ABS(fdiff)
                0068
                0069          preco = 2. * fdiff / (gnorm0*gnorm0)
                0070          do i = 1, nn
                0071             dd(i)    = -gg(i)*preco
                0072          end do
                0073
                0074          if (lphprint)
                0075      &        print *, 'pathei-lsopt: first estimate of dd via ',
                0076      &        'fold - ff'
                0077
                0078 c-----------------------------------------
                0079 c use the matrix stored in [diag]
                0080 c and the (y,s) pairs
                0081 c-----------------------------------------
                0082
                0083          else
                0084
                0085             do i = 1, nn
                0086                dd(i) = -gg(i)
                0087             end do
                0088
                0089             if (jmax .gt. 0) then
                0090                call hessupd( nn, nupdate, dd, jmin, jmax, xdiff,
                0091      &              lphprint )
                0092             else
                0093                if (lphprint)
                0094      &              print *, 'pathei-lsopt: no hessupd for first optim.'
                0095             end if
                0096
                0097          endif
                0098
                0099 c-----------------------------------------
                0100 c check whether new direction is a descent one
                0101 c-----------------------------------------
ae287463fd Patr*0102          dotdg = DDOT( nn, dd, 1, gg, 1 )
4cee17c1be Patr*0103          if (dotdg .ge. 0.0) then
                0104             ifail = 4
                0105             goto 999
                0106          end if
                0107
                0108 c----------------------------------
                0109 c declare arguments
                0110 c----------------------------------
                0111
                0112       tmin = 0.
                0113       do i = 1, nn
                0114          tmin = max( tmin, abs(dd(i)) )
                0115       end do
                0116       tmin = epsx/tmin
                0117
                0118 c----------------------------------
                0119 c make sure that t is between
                0120 c tmin and tmax
                0121 c----------------------------------
                0122
                0123       tact  = 1.0
                0124       tmax = 1.0e+10
                0125       if (tact.le.tmin) then
                0126          tact = tmin
                0127          if (tact.gt.tmax) then
                0128             tmin = tmax
                0129          endif
                0130       endif
                0131
                0132       if (tact .gt. tmax) then
                0133           tact = tmax
                0134           ifail = 7
                0135       endif
                0136
                0137 c----------------------------------
                0138 c compute new x
                0139 c----------------------------------
                0140       do i = 1, nn
                0141          xdiff(i) = xx(i) + tact*dd(i)
                0142       end do
                0143
                0144 c----------------------------------
                0145 c save new x to file for offline version
                0146 c----------------------------------
                0147
                0148  999  continue
                0149
                0150       return
                0151
                0152       end