/[MITgcm]/MITgcm_contrib/ecco_utils/lbfgs_jpl_version/lsopt.2/hessupd.F
ViewVC logotype

Contents of /MITgcm_contrib/ecco_utils/lbfgs_jpl_version/lsopt.2/hessupd.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Wed Apr 3 23:36:07 2013 UTC (12 years, 4 months ago) by heimbach
Branch: MAIN
CVS Tags: HEAD
Add L-BFGS code adapted to ECCO Production by JPL

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 include 'mpif.h'
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 integer status(MPI_STATUS_SIZE),ierr
41 integer myid, nprocs,mystart,myend
42 common /mpi_parm/ myid,nprocs,mystart,myend
43
44 c------------------------------------
45 c declare local variables
46 c------------------------------------
47 external DDOT
48 double precision DDOT
49
50 integer jfin, i, j, jp
51 double precision r
52 double precision temp
53
54 c------------------------------------
55 c initialization
56 c------------------------------------
57 jfin = jmax
58
59 if (lphprint .and. myid .eq. 0)
60 & print *, 'pathei-lsopt: in hessupd; ',
61 & 'jmin, jmax, mupd:', jmin, jmax, mupd
62
63
64 if (jfin.lt.jmin) jfin = jmax+mupd
65
66 c------------------------------------
67 c compute right hand side
68 c------------------------------------
69 do j = jfin,jmin,-1
70
71 if (lphprint .and. myid .eq. 0)
72 & print *, 'pathei-lsopt: in hessupd; loop ',
73 & 'j,jfin,jmin = ', j,jfin,jmin
74
75 jp = j
76 if (jp.gt.mupd) jp = jp-mupd
77 call dostore( nn, xdiff, .false., 2*jp+3 )
78
79 temp = DDOT( nn, dd, 1, xdiff, 1 )
80 call MPI_ALLREDUCE(temp,r,1,MPI_DOUBLE_PRECISION,MPI_SUM,
81 & MPI_COMM_WORLD,ierr)
82 ! r = DDOT( nn, dd, 1, xdiff,1 )
83
84 call dostore( nn, xdiff, .false., 2*jp+2 )
85 alpha(jp) = r
86 do i = 1, nn
87 dd(i) = dd(i) - r*xdiff(i)
88 end do
89 end do
90
91 c------------------------------------
92 c multiply precondition matrix
93 c------------------------------------
94 if (mupd .ne. 0) then
95 call dostore( nn, xdiff, .false., 3 )
96 do i = 1, nn
97 dd(i) = dd(i)*xdiff(i)
98 end do
99 end if
100
101 c------------------------------------
102 c compute left hand side
103 c------------------------------------
104 do j = jmin,jfin
105 jp = j
106 if (jp .gt. mupd) jp = jp-mupd
107 call dostore( nn, xdiff, .false., 2*jp+2 )
108
109 temp = DDOT( nn, dd, 1, xdiff, 1 )
110 call MPI_ALLREDUCE(temp,r,1,MPI_DOUBLE_PRECISION,MPI_SUM,
111 & MPI_COMM_WORLD,ierr)
112 r = alpha(jp) - r
113 ! r = alpha(jp) - DDOT( nn, dd,1 , xdiff, 1 )
114
115 call dostore( nn, xdiff, .false., 2*jp+3 )
116 do i = 1, nn
117 dd(i) = dd(i) + r*xdiff(i)
118 end do
119 end do
120
121 return
122
123 end

  ViewVC Help
Powered by ViewVC 1.1.22