1 |
|
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 |
#include "blas1.h" |
49 |
|
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 |
external DDOT |
89 |
double precision DDOT |
90 |
|
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 |
fp = DDOT( nn, dd, 1, gg, 1 ) |
132 |
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 |