MITgcm/lsopt/hessupd.F


      subroutine hessupd( nn, mupd, dd, jmin, jmax, xdiff, lphprint )

c     ==================================================================
c     SUBROUTINE hessupd
c     ==================================================================
c
c     o controls update of descent vector using available
c       approximation of Hessian Matrix based on gradients of
c       previous iterations
c
c     o Reference: J.C. Gilbert & C. Lemarechal
c                  Some numerical experiments with variable-storage
c                  quasi-Newton algorithms
c                  Mathematical Programming 45 (1989), pp. 407-435
c
c     o started: ??? not reproducible
c
c     o changed: Patrick Heimbach, MIT/EAPS
c                24-Feb-2000: 
c                   - changed some variable names to be consistent
c                     with routine lsoptv, lsline;
c
c     o Version: 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS
c
c     ==================================================================
c     SUBROUTINE hessupd
c     ==================================================================

      implicit none

#include "blas1.h"

c------------------------------------
c declare arguments
c------------------------------------
      integer nn, mupd, jmin, jmax
      double precision dd(nn), alpha(100), xdiff(nn)
      logical lphprint

c------------------------------------
c declare local variables
c------------------------------------
      external DDOT
      double precision     DDOT

      integer jfin, i, j, jp
      double precision    r

c------------------------------------
c initialization
c------------------------------------
      jfin = jmax

      if (lphprint) 
     &     print *, 'pathei-lsopt: in hessupd; ', 
     &     'jmin, jmax, mupd:', jmin, jmax, mupd

      if (jfin.lt.jmin) jfin = jmax+mupd

c------------------------------------
c compute right hand side
c------------------------------------
      do j = jfin,jmin,-1

         if (lphprint) 
     &        print *, 'pathei-lsopt: in hessupd; loop ',
     &        'j,jfin,jmin = ', j,jfin,jmin

         jp = j
         if (jp.gt.mupd) jp = jp-mupd
         call dostore( nn, xdiff, .false., 2*jp+3 )
         r = DDOT( nn, dd, 1, xdiff,1 )
         call dostore( nn, xdiff, .false., 2*jp+2 )
         alpha(jp) = r
         do i = 1, nn
            dd(i) = dd(i) - r*xdiff(i)
         end do
      end do

c------------------------------------
c multiply precondition matrix
c------------------------------------
      if (mupd .ne. 0) then
         call dostore( nn, xdiff, .false., 3 )
         do i = 1, nn
            dd(i) = dd(i)*xdiff(i)
         end do
      end if

c------------------------------------
c compute left hand side
c------------------------------------
      do j = jmin,jfin
         jp = j
         if (jp .gt. mupd) jp = jp-mupd
         call dostore( nn, xdiff, .false., 2*jp+2 )
         r = alpha(jp) - DDOT( nn, dd,1 , xdiff, 1 )
         call dostore( nn, xdiff, .false., 2*jp+3 )
         do i = 1, nn
            dd(i) = dd(i) + r*xdiff(i)
         end do
      end do

      return

      end
1
2	subroutine hessupd( nn, mupd, dd, jmin, jmax, xdiff, lphprint )
3
4	c ==================================================================
5	c SUBROUTINE hessupd
6	c ==================================================================
7	c
8	c o controls update of descent vector using available
9	c approximation of Hessian Matrix based on gradients of
10	c previous iterations
11	c
12	c o Reference: J.C. Gilbert & C. Lemarechal
13	c Some numerical experiments with variable-storage
14	c quasi-Newton algorithms
15	c Mathematical Programming 45 (1989), pp. 407-435
16	c
17	c o started: ??? not reproducible
18	c
19	c o changed: Patrick Heimbach, MIT/EAPS
20	c 24-Feb-2000:
21	c - changed some variable names to be consistent
22	c with routine lsoptv, lsline;
23	c
24	c o Version: 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS
25	c
26	c ==================================================================
27	c SUBROUTINE hessupd
28	c ==================================================================
29
30	implicit none
31
32	#include "blas1.h"
33
34	c------------------------------------
35	c declare arguments
36	c------------------------------------
37	integer nn, mupd, jmin, jmax
38	double precision dd(nn), alpha(100), xdiff(nn)
39	logical lphprint
40
41	c------------------------------------
42	c declare local variables
43	c------------------------------------
44	external DDOT
45	double precision DDOT
46
47	integer jfin, i, j, jp
48	double precision r
49
50	c------------------------------------
51	c initialization
52	c------------------------------------
53	jfin = jmax
54
55	if (lphprint)
56	& print *, 'pathei-lsopt: in hessupd; ',
57	& 'jmin, jmax, mupd:', jmin, jmax, mupd
58
59	if (jfin.lt.jmin) jfin = jmax+mupd
60
61	c------------------------------------
62	c compute right hand side
63	c------------------------------------
64	do j = jfin,jmin,-1
65
66	if (lphprint)
67	& print *, 'pathei-lsopt: in hessupd; loop ',
68	& 'j,jfin,jmin = ', j,jfin,jmin
69
70	jp = j
71	if (jp.gt.mupd) jp = jp-mupd
72	call dostore( nn, xdiff, .false., 2*jp+3 )
73	r = DDOT( nn, dd, 1, xdiff,1 )
74	call dostore( nn, xdiff, .false., 2*jp+2 )
75	alpha(jp) = r
76	do i = 1, nn
77	dd(i) = dd(i) - r*xdiff(i)
78	end do
79	end do
80
81	c------------------------------------
82	c multiply precondition matrix
83	c------------------------------------
84	if (mupd .ne. 0) then
85	call dostore( nn, xdiff, .false., 3 )
86	do i = 1, nn
87	dd(i) = dd(i)*xdiff(i)
88	end do
89	end if
90
91	c------------------------------------
92	c compute left hand side
93	c------------------------------------
94	do j = jmin,jfin
95	jp = j
96	if (jp .gt. mupd) jp = jp-mupd
97	call dostore( nn, xdiff, .false., 2*jp+2 )
98	r = alpha(jp) - DDOT( nn, dd,1 , xdiff, 1 )
99	call dostore( nn, xdiff, .false., 2*jp+3 )
100	do i = 1, nn
101	dd(i) = dd(i) + r*xdiff(i)
102	end do
103	end do
104
105	return
106
107	end