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

Annotation of /MITgcm_contrib/ecco_utils/lbfgs_jpl_version/lsopt.2/lsline.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 lsline(
3     & simul
4     & , nn, ifail, lphprint
5     & , ifunc, nfunc
6     & , ff, dotdg
7     & , sflag, tflag, tact, epsx
8     & , dd, gg, xx, xdiff
9     & )
10    
11     c ==================================================================
12     c SUBROUTINE lsline
13     c ==================================================================
14     c
15     c o line search algorithm for determining control vector update;
16     c After computing updated control vector from given gradient,
17     c a forward and adjoint model run are performed (simul.F)
18     c using the updated control vector.
19     c Tests are then applied to see whether solution has improved.
20     c
21     c o Reference: J.C. Gilbert & C. Lemarechal
22     c Some numerical experiments with variable-storage
23     c quasi-Newton algorithms
24     c Mathematical Programming 45 (1989), pp. 407-435
25     c
26     c o started: ??? not reproducible
27     c
28     c o changed: Benny Cheng, NASA/JPL
29     c
30     c ==================================================================
31     c SUBROUTINE lsline
32     c ==================================================================
33    
34     #include "blas1.h"
35    
36     implicit none
37     include 'mpif.h'
38     c----------------------------------
39     c declare arguments
40     c----------------------------------
41     integer nmax
42     parameter( nmax = MAX_INDEPEND )
43     integer nn, ifail, ifunc, nfunc
44     double precision ff, dotdg, tact, epsx
45     double precision xx(nn), gg(nn), xdiff(nn),dd(nn)
46     logical tflag,sflag,isforward
47     logical lphprint
48    
49     external simul
50    
51     integer pidlen,myindx(2)
52     integer status(MPI_STATUS_SIZE),ierr
53     common /mpi_parm/myid,nprocs,mystart,myend
54     integer myid,nprocs,mystart,myend
55     c----------------------------------
56     c declare local variables
57     c----------------------------------
58    
59     double precision xpara1, xpara2
60     parameter( xpara1 = 0.001D0, xpara2=0.9D0 )
61    
62     double precision factor
63     parameter( factor = 0.2 )
64    
65     double precision barmax
66     parameter( barmax = 0.3 )
67     double precision barmul
68     parameter( barmul = 5.0 )
69     double precision barmin
70     parameter( barmin = 0.01 )
71     double precision tmax
72     parameter( tmax = 5. )
73    
74     integer i
75    
76     double precision tg, fg, td, ta
77     double precision fa, fpa, fp
78     double precision fnew, fdiff,ftemp
79     double precision z, test, barr
80     double precision left, right, told
81     double precision t,a,b
82     double precision xtemp(nn),gtemp(nn),gold(nn)
83     external DDOT,DNRM2
84     double precision DDOT, DNRM2
85    
86     double precision temp,gnorm0
87     real*8,allocatable:: vv(:)
88     integer, allocatable:: displs(:)
89     integer, allocatable:: counts(:)
90    
91     c----------------------------------
92     c check parameters
93     c----------------------------------
94    
95     if ( (nn.le.0)
96     & .or. (xpara1.le.0.0) .or. (xpara1.ge.0.5)
97     & .or. (xpara2.le.xpara1) .or. (xpara2.ge.1.0) ) then
98     ifail = 9
99     go to 999
100     endif
101    
102     c----------------------------------
103     c initialization
104     c----------------------------------
105    
106     c=======================================================================
107     c begin of simulation iter.
108     c=======================================================================
109    
110     if (lphprint .and. myid .eq. 0)
111     & print *, 'lsopt: simul.'
112    
113     c------------------------------------
114     c compute cost function and gradient
115     c------------------------------------
116    
117     gold = gg
118     temp = DNRM2(nn,gg,1)
119     call MPI_ALLREDUCE(temp*temp,gnorm0,1,MPI_DOUBLE_PRECISION,
120     & MPI_SUM,MPI_COMM_WORLD,ierr)
121     gnorm0 = dsqrt(gnorm0)
122    
123     if (sflag ) then
124    
125     !get estimated cost from hessian update control, this is the trial step
126     isforward = .true.
127     call simul( 0, isforward,nn, xdiff, fnew, gg )
128    
129     xtemp = xdiff
130     ! new direction vector dd
131    
132     temp = DDOT( nn, xdiff - xx, 1, gold, 1 )
133     call MPI_ALLREDUCE(temp,b,1,MPI_DOUBLE_PRECISION,MPI_SUM,
134     & MPI_COMM_WORLD,ierr)
135     dotdg = b
136    
137     c-----------------------------------------
138     c quadratic line search
139     c-----------------------------------------
140    
141     a = fnew - ff - b
142     t = -b/(2.*a)
143     t = dmax1(t,0.D0)
144    
145     if (t <= 1.d-2 ) then
146     tact = 1.
147     xdiff = xx - 2.*(.01*ff)*gold*tact/(gnorm0*gnorm0)
148     tflag = .false.
149    
150     if (lphprint .and. myid .eq. 0)
151     & print *, 'lsopt: using steepest descent, t value ',t,'< .01'
152    
153     else
154    
155     xdiff = (1. - t)*xx + t*xdiff
156     tact = t
157     tflag = .true.
158     endif
159    
160     if (lphprint .and. myid .eq. 0)
161     & print *, 'lsopt: trial stepsize =',t
162    
163     ! no hessian update , get new adjoint run
164     goto 999
165    
166     else
167    
168     !no line search
169     !read new ctrl vector and new cost gradient and new cost
170     isforward = .false.
171     call simul( 0, isforward, nn, xdiff, fnew, gg )
172    
173     if (tflag) then
174    
175     !read in trial step adjoint output
176     call simul( ifunc+1, isforward, nn, xtemp, ftemp, gtemp)
177    
178     dd = xtemp - xx
179     c-----------------------------------------
180     c apply 1st Wolfe test
181     c-----------------------------------------
182     fdiff = fnew - ff
183    
184     temp = DDOT( nn, dd, 1, gold, 1 )
185     call MPI_ALLREDUCE(temp,dotdg,1,MPI_DOUBLE_PRECISION,MPI_SUM,
186     & MPI_COMM_WORLD,ierr)
187    
188     if (fdiff .gt. tact*xpara1*dotdg) then
189     ifail = 7
190     if (lphprint .and. myid .eq. 0)
191     & print *, 'lsopt: fails wolfe test 1'
192    
193     endif
194    
195     c-----------------------------------------
196     c 1st Wolfe test 1 ok, apply 2nd Wolf test
197     c-----------------------------------------
198    
199     temp = DDOT( nn, dd, 1, gg, 1 )
200     call MPI_ALLREDUCE(temp,fp,1,MPI_DOUBLE_PRECISION,MPI_SUM,
201     & MPI_COMM_WORLD,ierr)
202    
203     if (fp .lt. xpara2*dotdg) then
204     ifail = 7
205     if (lphprint .and. myid .eq. 0)
206     & print *, 'lsopt: fails wolfe test 2'
207     endif
208    
209     c-----------------------------------------
210     c 2nd Wolfe test 2 ok, donc pas serieux, on sort
211     c-----------------------------------------
212    
213     if (ifail .eq. 7 ) then
214    
215     if (ifunc .lt. nfunc) then
216     t = dmin1(tact,tmax)
217     if (t < 5.D0 ) t = .5*t
218     xdiff = (1. - t)*xx + t*xtemp
219     tact = t
220     ifunc = ifunc + 1
221     goto 999
222     endif
223    
224     endif !end ifail = 7
225    
226     else
227    
228     if ( fnew > .999D0*ff ) then
229    
230     if (ifunc .lt. nfunc) then
231     tact = tact*.5
232     xdiff = xx - 2.*(.01*ff)*gold*tact/(gnorm0*gnorm0)
233     ifunc = ifunc + 1
234     goto 999
235     endif
236    
237     endif
238    
239     endif ! end of tflag
240    
241     dd = xdiff - xx
242     ff = fnew
243     xx =xdiff
244     ifunc = 0
245     tact = 1.
246     endif ! end of sflag
247     c-----------------------------------------
248     c end of routine
249     c-----------------------------------------
250     999 continue
251    
252     return
253    
254     end

  ViewVC Help
Powered by ViewVC 1.1.22