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

Annotation 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 - (hide 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 heimbach 1.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