1 |
heimbach |
1.2 |
|
2 |
|
|
subroutine lsline( |
3 |
|
|
& simul |
4 |
|
|
& , nn, ifail, lphprint |
5 |
|
|
& , ifunc, nfunc |
6 |
|
|
& , ff, dotdg |
7 |
|
|
& , tmin, tmax, 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: Patrick Heimbach, MIT/EAPS |
29 |
|
|
c |
30 |
|
|
c o Version: 2.0, 24-Feb-2000: Patrick Heimbach, MIT/EAPS |
31 |
|
|
c - severe changes in structure including various |
32 |
|
|
c shifts of variables which are only used in this |
33 |
|
|
c routine |
34 |
|
|
c - number of 3 control flags for error handling |
35 |
|
|
c (indic, moderl, ifail) reduced to 1 (ifail) |
36 |
|
|
c and homogenized with routine lsoptv |
37 |
|
|
c |
38 |
|
|
c o Version: 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS |
39 |
|
|
c - initial computation of tact and |
40 |
|
|
c xdiff = xx + tact*dd |
41 |
|
|
c moved to new routine lsupdxx |
42 |
|
|
c tmin, tmax, tact needed as parameters |
43 |
|
|
c |
44 |
|
|
c ================================================================== |
45 |
|
|
c SUBROUTINE lsline |
46 |
|
|
c ================================================================== |
47 |
|
|
|
48 |
heimbach |
1.4 |
#include "blas1.h" |
49 |
heimbach |
1.2 |
|
50 |
|
|
implicit none |
51 |
|
|
|
52 |
|
|
c---------------------------------- |
53 |
|
|
c declare arguments |
54 |
|
|
c---------------------------------- |
55 |
|
|
integer nn, ifail, ifunc, nfunc |
56 |
|
|
double precision ff, dotdg, tmin, tmax, tact, epsx |
57 |
|
|
double precision xx(nn), dd(nn), gg(nn), xdiff(nn) |
58 |
|
|
|
59 |
|
|
logical lphprint |
60 |
|
|
|
61 |
|
|
external simul |
62 |
|
|
|
63 |
|
|
c---------------------------------- |
64 |
|
|
c declare local variables |
65 |
|
|
c---------------------------------- |
66 |
|
|
|
67 |
|
|
double precision xpara1, xpara2 |
68 |
|
|
parameter( xpara1 = 0.0001, xpara2=0.9 ) |
69 |
|
|
|
70 |
|
|
double precision factor |
71 |
|
|
parameter( factor = 0.2 ) |
72 |
|
|
|
73 |
|
|
double precision barmax |
74 |
|
|
parameter( barmax = 0.3 ) |
75 |
|
|
double precision barmul |
76 |
|
|
parameter( barmul = 5.0 ) |
77 |
|
|
double precision barmin |
78 |
|
|
parameter( barmin = 0.01 ) |
79 |
|
|
|
80 |
|
|
integer i, indic |
81 |
|
|
|
82 |
|
|
double precision tg, fg, td, ta |
83 |
|
|
double precision fa, fpa, fp |
84 |
|
|
double precision fnew, fdiff |
85 |
|
|
double precision z, test, barr |
86 |
|
|
double precision left, right, told |
87 |
|
|
|
88 |
heimbach |
1.3 |
external DDOT |
89 |
|
|
double precision DDOT |
90 |
heimbach |
1.2 |
|
91 |
|
|
c---------------------------------- |
92 |
|
|
c check parameters |
93 |
|
|
c---------------------------------- |
94 |
|
|
|
95 |
|
|
if ( (nn.le.0) |
96 |
|
|
& .or. (dotdg.ge.0.0) |
97 |
|
|
& .or. (xpara1.le.0.0) .or. (xpara1.ge.0.5) |
98 |
|
|
& .or. (xpara2.le.xpara1) .or. (xpara2.ge.1.0) ) then |
99 |
|
|
ifail = 9 |
100 |
|
|
go to 999 |
101 |
|
|
endif |
102 |
|
|
|
103 |
|
|
c---------------------------------- |
104 |
|
|
c initialization |
105 |
|
|
c---------------------------------- |
106 |
|
|
indic = 0 |
107 |
|
|
|
108 |
|
|
barr = barmin |
109 |
|
|
fg = ff |
110 |
|
|
fa = ff |
111 |
|
|
fpa = dotdg |
112 |
|
|
|
113 |
|
|
td = 0.0 |
114 |
|
|
tg = 0.0 |
115 |
|
|
ta = 0.0 |
116 |
|
|
|
117 |
|
|
c======================================================================= |
118 |
|
|
c begin of simulation iter. |
119 |
|
|
c======================================================================= |
120 |
|
|
|
121 |
|
|
do ifunc = 1, nfunc |
122 |
|
|
|
123 |
|
|
if (lphprint) |
124 |
|
|
& print *, 'pathei-lsopt: ', ifunc, ' simul.' |
125 |
|
|
|
126 |
|
|
c------------------------------------ |
127 |
|
|
c compute cost function and gradient |
128 |
|
|
c------------------------------------ |
129 |
|
|
call simul( indic, nn, xdiff, fnew, gg ) |
130 |
|
|
|
131 |
heimbach |
1.3 |
fp = DDOT( nn, dd, 1, gg, 1 ) |
132 |
heimbach |
1.2 |
fdiff = fnew - ff |
133 |
|
|
|
134 |
|
|
c----------------------------------------- |
135 |
|
|
c apply 1st Wolfe test |
136 |
|
|
c----------------------------------------- |
137 |
|
|
if (fdiff .gt. tact*xpara1*dotdg) then |
138 |
|
|
td = tact |
139 |
|
|
ifail = 0 |
140 |
|
|
go to 500 |
141 |
|
|
endif |
142 |
|
|
|
143 |
|
|
c----------------------------------------- |
144 |
|
|
c 1st Wolfe test 1 ok, apply 2nd Wolf test |
145 |
|
|
c----------------------------------------- |
146 |
|
|
if (fp .gt. xpara2*dotdg) then |
147 |
|
|
ifail = 0 |
148 |
|
|
go to 320 |
149 |
|
|
endif |
150 |
|
|
|
151 |
|
|
if (ifail.eq.0) go to 350 |
152 |
|
|
|
153 |
|
|
c----------------------------------------- |
154 |
|
|
c 2nd Wolfe test 2 ok, donc pas serieux, on sort |
155 |
|
|
c----------------------------------------- |
156 |
|
|
|
157 |
|
|
320 continue |
158 |
|
|
|
159 |
|
|
ff = fnew |
160 |
|
|
do i = 1, nn |
161 |
|
|
xx(i) = xdiff(i) |
162 |
|
|
end do |
163 |
|
|
cph( |
164 |
|
|
if (lphprint) |
165 |
|
|
& print *, 'pathei-lsopt: no inter-/extrapolation in lsline' |
166 |
|
|
cph) |
167 |
|
|
go to 999 |
168 |
|
|
|
169 |
|
|
|
170 |
|
|
c----------------------------------------- |
171 |
|
|
c extrapolation |
172 |
|
|
c----------------------------------------- |
173 |
|
|
|
174 |
|
|
350 continue |
175 |
|
|
|
176 |
|
|
tg = tact |
177 |
|
|
fg = fnew |
178 |
|
|
if (td .ne. 0.0) go to 500 |
179 |
|
|
|
180 |
|
|
told = tact |
181 |
|
|
left = (1.0+barmin)*tact |
182 |
|
|
right = 10.0*tact |
183 |
|
|
call cubic( tact, fnew, fp, ta, fa, fpa, left, right ) |
184 |
|
|
ta = told |
185 |
|
|
if (tact.ge.tmax) then |
186 |
|
|
ifail = 7 |
187 |
|
|
tact = tmax |
188 |
|
|
endif |
189 |
|
|
|
190 |
|
|
if (lphprint) |
191 |
|
|
& print *, 'pathei-lsopt: extrapolation: ', |
192 |
|
|
& 'td, tg, tact, ifail = ', td, tg, tact, ifail |
193 |
|
|
|
194 |
|
|
go to 900 |
195 |
|
|
|
196 |
|
|
c----------------------------------------- |
197 |
|
|
c interpolation |
198 |
|
|
c----------------------------------------- |
199 |
|
|
|
200 |
|
|
500 continue |
201 |
|
|
|
202 |
|
|
test = barr*(td-tg) |
203 |
|
|
left = tg+test |
204 |
|
|
right = td-test |
205 |
|
|
told = tact |
206 |
|
|
call cubic( tact, fnew, fp, ta, fa, fpa, left, right ) |
207 |
|
|
ta = told |
208 |
|
|
if (tact.gt.left .and. tact.lt.right) then |
209 |
|
|
barr = dmax1( barmin, barr/barmul ) |
210 |
|
|
else |
211 |
|
|
barr = dmin1( barmul*barr, barmax ) |
212 |
|
|
endif |
213 |
|
|
|
214 |
|
|
if (lphprint) |
215 |
|
|
& print *, 'pathei-lsopt: interpolation: ', |
216 |
|
|
& 'td, tg, tact, ifail = ', td, tg, tact, ifail |
217 |
|
|
|
218 |
|
|
c----------------------------------------- |
219 |
|
|
c end of iteration loop |
220 |
|
|
c - tact peut etre bloque sur tmax |
221 |
|
|
c (venant de lextrapolation avec ifail=7) |
222 |
|
|
c----------------------------------------- |
223 |
|
|
|
224 |
|
|
900 continue |
225 |
|
|
|
226 |
|
|
fa = fnew |
227 |
|
|
fpa = fp |
228 |
|
|
|
229 |
|
|
c----------------------------------------- |
230 |
|
|
c --- faut-il continuer ? |
231 |
|
|
c----------------------------------------- |
232 |
|
|
if (td .eq. 0.0) go to 950 |
233 |
|
|
if (td-tg .lt. tmin) go to 920 |
234 |
|
|
|
235 |
|
|
c----------------------------------------- |
236 |
|
|
c limit due to machine precision |
237 |
|
|
c----------------------------------------- |
238 |
|
|
do i = 1, nn |
239 |
|
|
z = xx(i) + tact*dd(i) |
240 |
|
|
if ((z.ne.xx(i)) .and. (z.ne.xdiff(i))) go to 950 |
241 |
|
|
end do |
242 |
|
|
|
243 |
|
|
c----------------------------------------- |
244 |
|
|
c arret sur dxmin ou de secours |
245 |
|
|
c----------------------------------------- |
246 |
|
|
920 continue |
247 |
|
|
ifail = 8 |
248 |
|
|
|
249 |
|
|
c----------------------------------------- |
250 |
|
|
c si tg=0, xx = xx_depart, |
251 |
|
|
c sinon on prend xx=x_left qui fait decroitre f |
252 |
|
|
c----------------------------------------- |
253 |
|
|
if (tg .ne. 0.0) then |
254 |
|
|
ff = fg |
255 |
|
|
do i = 1, nn |
256 |
|
|
xx(i) = xx(i) + tg*dd(i) |
257 |
|
|
end do |
258 |
|
|
endif |
259 |
|
|
|
260 |
|
|
go to 999 |
261 |
|
|
|
262 |
|
|
c----------------------------------------- |
263 |
|
|
c update vector for new simulation |
264 |
|
|
c----------------------------------------- |
265 |
|
|
950 continue |
266 |
|
|
|
267 |
|
|
do i = 1, nn |
268 |
|
|
xdiff(i) = xx(i) + tact*dd(i) |
269 |
|
|
end do |
270 |
|
|
|
271 |
|
|
c======================================================================= |
272 |
|
|
c end of simulation iter. |
273 |
|
|
c======================================================================= |
274 |
|
|
end do |
275 |
|
|
|
276 |
|
|
c----------------------------------------- |
277 |
|
|
c too many function evaluations |
278 |
|
|
c----------------------------------------- |
279 |
|
|
ifail = 6 |
280 |
|
|
ff = fg |
281 |
|
|
do i = 1, nn |
282 |
|
|
xx(i) = xx(i) + tg*dd(i) |
283 |
|
|
end do |
284 |
|
|
|
285 |
|
|
c----------------------------------------- |
286 |
|
|
c end of routine |
287 |
|
|
c----------------------------------------- |
288 |
|
|
999 continue |
289 |
|
|
|
290 |
|
|
return |
291 |
|
|
|
292 |
|
|
end |