MITgcm/lsopt/hessupd.F

4cee17c1be Patr*0001
                0002       subroutine hessupd( nn, mupd, dd, jmin, jmax, xdiff, lphprint )
                0003
                0004 c     ==================================================================
                0005 c     SUBROUTINE hessupd
                0006 c     ==================================================================
                0007 c
                0008 c     o controls update of descent vector using available
                0009 c       approximation of Hessian Matrix based on gradients of
                0010 c       previous iterations
                0011 c
                0012 c     o Reference: J.C. Gilbert & C. Lemarechal
                0013 c                  Some numerical experiments with variable-storage
                0014 c                  quasi-Newton algorithms
                0015 c                  Mathematical Programming 45 (1989), pp. 407-435
                0016 c
                0017 c     o started: ??? not reproducible
                0018 c
                0019 c     o changed: Patrick Heimbach, MIT/EAPS
                0020 c                24-Feb-2000:
                0021 c                   - changed some variable names to be consistent
                0022 c                     with routine lsoptv, lsline;
                0023 c
                0024 c     o Version: 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS
                0025 c
                0026 c     ==================================================================
                0027 c     SUBROUTINE hessupd
                0028 c     ==================================================================
                0029
                0030       implicit none
                0031
f50e6c1777 Patr*0032 #include "blas1.h"
4cee17c1be Patr*0033
                0034 c------------------------------------
                0035 c declare arguments
                0036 c------------------------------------
                0037       integer nn, mupd, jmin, jmax
                0038       double precision dd(nn), alpha(100), xdiff(nn)
                0039       logical lphprint
                0040
                0041 c------------------------------------
                0042 c declare local variables
                0043 c------------------------------------
ae287463fd Patr*0044       external DDOT
                0045       double precision     DDOT
4cee17c1be Patr*0046
                0047       integer jfin, i, j, jp
                0048       double precision    r
                0049
                0050 c------------------------------------
                0051 c initialization
                0052 c------------------------------------
                0053       jfin = jmax
                0054
                0055       if (lphprint)
                0056      &     print *, 'pathei-lsopt: in hessupd; ',
                0057      &     'jmin, jmax, mupd:', jmin, jmax, mupd
                0058
                0059       if (jfin.lt.jmin) jfin = jmax+mupd
                0060
                0061 c------------------------------------
                0062 c compute right hand side
                0063 c------------------------------------
                0064       do j = jfin,jmin,-1
                0065
                0066          if (lphprint)
                0067      &        print *, 'pathei-lsopt: in hessupd; loop ',
                0068      &        'j,jfin,jmin = ', j,jfin,jmin
                0069
                0070          jp = j
                0071          if (jp.gt.mupd) jp = jp-mupd
                0072          call dostore( nn, xdiff, .false., 2*jp+3 )
ae287463fd Patr*0073          r = DDOT( nn, dd, 1, xdiff,1 )
4cee17c1be Patr*0074          call dostore( nn, xdiff, .false., 2*jp+2 )
                0075          alpha(jp) = r
                0076          do i = 1, nn
                0077             dd(i) = dd(i) - r*xdiff(i)
                0078          end do
                0079       end do
                0080
                0081 c------------------------------------
                0082 c multiply precondition matrix
                0083 c------------------------------------
                0084       if (mupd .ne. 0) then
                0085          call dostore( nn, xdiff, .false., 3 )
                0086          do i = 1, nn
                0087             dd(i) = dd(i)*xdiff(i)
                0088          end do
                0089       end if
                0090
                0091 c------------------------------------
                0092 c compute left hand side
                0093 c------------------------------------
                0094       do j = jmin,jfin
                0095          jp = j
                0096          if (jp .gt. mupd) jp = jp-mupd
                0097          call dostore( nn, xdiff, .false., 2*jp+2 )
ae287463fd Patr*0098          r = alpha(jp) - DDOT( nn, dd,1 , xdiff, 1 )
4cee17c1be Patr*0099          call dostore( nn, xdiff, .false., 2*jp+3 )
                0100          do i = 1, nn
                0101             dd(i) = dd(i) + r*xdiff(i)
                0102          end do
                0103       end do
                0104
                0105       return
                0106
                0107       end