1 |
heimbach |
1.2 |
|
2 |
|
|
subroutine lsopt_top( nn, xx, ff, gg, simul, optline |
3 |
|
|
$ , epsx, fmin, epsg |
4 |
|
|
$ , iprint, itmax, nfunc, nupdate |
5 |
|
|
$ , dd, gold, xdiff |
6 |
|
|
$ , loffline |
7 |
|
|
$ , ifail ) |
8 |
|
|
|
9 |
|
|
c ================================================================== |
10 |
|
|
c SUBROUTINE lsopt_top |
11 |
|
|
c ================================================================== |
12 |
|
|
c |
13 |
|
|
c o uses a set of control variables, their adjoint variables, |
14 |
|
|
c and a cost function value |
15 |
|
|
c to compute an improved set of controls with respect to the |
16 |
|
|
c cost function via a |
17 |
|
|
c variable-storage Quasi-Newton method |
18 |
|
|
c |
19 |
|
|
c o Reference: J.C. Gilbert & C. Lemarechal |
20 |
|
|
c Some numerical experiments with variable-storage |
21 |
|
|
c quasi-Newton algorithms |
22 |
|
|
c Mathematical Programming 45 (1989), pp. 407-435 |
23 |
|
|
c |
24 |
|
|
c o started: ??? not reproducible |
25 |
|
|
c |
26 |
|
|
c o changed: Patrick Heimbach, MIT/EAPS |
27 |
|
|
c |
28 |
|
|
c o Version: 2.0, 24-Feb-2000: |
29 |
|
|
c (Version 1.0 is considered as version |
30 |
|
|
c starting from which changes were made). |
31 |
|
|
c - severe changes in structure including various |
32 |
|
|
c bug-fixes to enable multi-optimization runs; |
33 |
|
|
c - routine lsoptw incorporated into lsoptv |
34 |
|
|
c - optimization iteration loop restructured |
35 |
|
|
c - complete restructuring of handling |
36 |
|
|
c cold start cases |
37 |
|
|
c - number of 3 control flags for error handling |
38 |
|
|
c (indic, moderl, ifail) reduced to 1 (ifail) |
39 |
|
|
c and homogenized with routine lsline |
40 |
|
|
c |
41 |
|
|
c o Version: 2.1, 29-Feb-2000: |
42 |
|
|
c - handling of case ifail = 6 changed; |
43 |
|
|
c leads to stop of optimization |
44 |
|
|
c (similar to case ifail = 8) |
45 |
|
|
c - logical lphprint included |
46 |
|
|
c |
47 |
|
|
c ================================================================== |
48 |
|
|
c SUBROUTINE lsopt_top |
49 |
|
|
c ================================================================== |
50 |
|
|
|
51 |
|
|
implicit none |
52 |
|
|
|
53 |
heimbach |
1.3 |
ccc#include <blas1.h> |
54 |
heimbach |
1.2 |
|
55 |
|
|
c----------------------------------------- |
56 |
|
|
c declare arguments |
57 |
|
|
c----------------------------------------- |
58 |
|
|
integer nn, iprint, itmax, nfunc, nupdate, ifail |
59 |
|
|
|
60 |
|
|
double precision xx(nn), ff, gg(nn), epsx, fmin, epsg |
61 |
|
|
double precision dd(nn), gold(nn), xdiff(nn) |
62 |
|
|
|
63 |
|
|
cph( |
64 |
|
|
integer phniter0, phniterold |
65 |
|
|
double precision phff |
66 |
|
|
COMMON /PH_OPTI/ phniter0, phniterold, phff |
67 |
|
|
cph) |
68 |
|
|
|
69 |
|
|
external simul, optline |
70 |
|
|
|
71 |
|
|
c----------------------------------------- |
72 |
|
|
C declare local variables |
73 |
|
|
c----------------------------------------- |
74 |
|
|
logical cold, lphprint, loffline |
75 |
|
|
parameter (lphprint = .true.) |
76 |
|
|
|
77 |
|
|
integer mm, mupd, jmin, jmax, indic, isize, REAL_BYTE |
78 |
|
|
integer i, iiter, ifunc |
79 |
|
|
|
80 |
|
|
double precision fupd |
81 |
|
|
double precision r1, tmin, tmax, tact, gnorm, gnorm0, eps1 |
82 |
|
|
double precision fold, ys |
83 |
|
|
double precision dotdg |
84 |
|
|
|
85 |
heimbach |
1.3 |
external DDOT, DNRM2, DSCAL |
86 |
|
|
double precision DDOT, DNRM2 |
87 |
heimbach |
1.2 |
|
88 |
|
|
c----------------------------------------- |
89 |
|
|
c parameters |
90 |
|
|
c----------------------------------------- |
91 |
|
|
|
92 |
|
|
double precision rmin |
93 |
|
|
parameter( rmin = 1.e-20 ) |
94 |
|
|
|
95 |
|
|
character*(*) iform |
96 |
|
|
parameter(iform='(i5,2x,1pe8.1,1x,i5,4x,1pe11.4,3(2x,1pe8.1))' ) |
97 |
|
|
|
98 |
|
|
c ================================================================== |
99 |
|
|
c |
100 |
|
|
c----------------------------------------- |
101 |
|
|
c initialization |
102 |
|
|
c----------------------------------------- |
103 |
|
|
cold = .true. |
104 |
|
|
fupd = 1.0e+10 |
105 |
|
|
indic = 0 |
106 |
|
|
tmin = 0. |
107 |
|
|
tmax = 1.0e+10 |
108 |
|
|
tact = 1.0 |
109 |
|
|
cph( |
110 |
|
|
phniterold = 0 |
111 |
|
|
cph) |
112 |
|
|
iiter = 0 |
113 |
|
|
eps1 = 1.0 |
114 |
|
|
ifunc = 0 |
115 |
|
|
ifail = 0 |
116 |
|
|
gnorm0 = 1. |
117 |
|
|
|
118 |
|
|
c----------------------------------------- |
119 |
|
|
c initialization for dd and dds |
120 |
|
|
c----------------------------------------- |
121 |
|
|
|
122 |
|
|
jmin = 1 |
123 |
|
|
jmax = 0 |
124 |
|
|
|
125 |
|
|
mm = nn |
126 |
|
|
mupd = nupdate |
127 |
|
|
|
128 |
heimbach |
1.7 |
REAL_BYTE = 4 |
129 |
heimbach |
1.2 |
isize = REAL_BYTE |
130 |
|
|
|
131 |
|
|
c----------------------------------------- |
132 |
|
|
c print information |
133 |
|
|
c----------------------------------------- |
134 |
|
|
if (iprint .ge. 1) then |
135 |
|
|
print '(2x,a)', |
136 |
|
|
$ '===============================================' |
137 |
|
|
print '(2x,a)', |
138 |
|
|
$ '=== O P T I M I Z A T I O N ===' |
139 |
|
|
print '(2x,a)', |
140 |
|
|
$ '===============================================' |
141 |
|
|
print '(a,i9)' |
142 |
|
|
$ , ' number of control variables.......', nn |
143 |
|
|
print '(a,e9.2)' |
144 |
|
|
$ , ' precision on x....................', epsx |
145 |
|
|
print '(a,e9.2)' |
146 |
|
|
$ , ' precision on g....................', epsg |
147 |
|
|
print '(a,e9.2)' |
148 |
|
|
$ , ' expected optimal function value...', fmin |
149 |
|
|
print '(a,i9)' |
150 |
|
|
$ , ' maximal number of iterations......', itmax |
151 |
|
|
print '(a,i9)' |
152 |
|
|
$ , ' maximal number of simulations.....', nfunc |
153 |
|
|
print '(a,i9)' |
154 |
|
|
$ , ' information level.................', iprint |
155 |
|
|
print '(a,i9)' |
156 |
|
|
$ , ' number of updates.................', nupdate |
157 |
|
|
print '(a,i9)' |
158 |
|
|
$ , ' size of used memory...............', 3*nn |
159 |
|
|
endif |
160 |
|
|
|
161 |
|
|
c----------------------------------------- |
162 |
|
|
c check arguments |
163 |
|
|
c----------------------------------------- |
164 |
|
|
|
165 |
|
|
if (nn .le. 0) then |
166 |
|
|
if (iprint.ge.1) then |
167 |
|
|
print '(a,i6)' , ' ERROR : n = ', nn |
168 |
|
|
endif |
169 |
|
|
ifail = 1 |
170 |
|
|
goto 999 |
171 |
|
|
endif |
172 |
|
|
|
173 |
|
|
if (itmax .lt. 0) then |
174 |
|
|
if (iprint.ge.1) then |
175 |
|
|
print '(a,i6)' , ' ERROR : itmax = ', itmax |
176 |
|
|
endif |
177 |
|
|
ifail = 1 |
178 |
|
|
goto 999 |
179 |
|
|
endif |
180 |
|
|
|
181 |
|
|
if (nfunc .le. 0) then |
182 |
|
|
if (iprint.ge.10) then |
183 |
|
|
print '(a,i6)' , ' ERROR : nfunc = ', nfunc |
184 |
|
|
endif |
185 |
|
|
ifail = 1 |
186 |
|
|
goto 999 |
187 |
|
|
endif |
188 |
|
|
|
189 |
|
|
if (epsx .le. 0.) then |
190 |
|
|
if (iprint.ge.1) then |
191 |
|
|
print '(a,e9.2)', ' ERROR : epsx = ', epsx |
192 |
|
|
endif |
193 |
|
|
ifail = 1 |
194 |
|
|
goto 999 |
195 |
|
|
endif |
196 |
|
|
|
197 |
|
|
if (epsg .le. 0.) then |
198 |
|
|
if (iprint.ge.1) then |
199 |
|
|
print '(a,e9.2)', ' ERROR : epsg = ', epsg |
200 |
|
|
endif |
201 |
|
|
ifail = 1 |
202 |
|
|
goto 999 |
203 |
|
|
endif |
204 |
|
|
|
205 |
|
|
if (epsg .gt. 1.) then |
206 |
|
|
if (iprint.ge.1) then |
207 |
|
|
print '(a,e9.2)', ' ERROR : epsg = ', epsg |
208 |
|
|
endif |
209 |
|
|
ifail = 1 |
210 |
|
|
goto 999 |
211 |
|
|
endif |
212 |
|
|
|
213 |
|
|
cph( |
214 |
|
|
print *, 'pathei: vor instore ' |
215 |
|
|
cph) |
216 |
|
|
call instore( mm, fupd, gnorm0, isize, mupd, jmin, jmax, cold, |
217 |
|
|
& ifail ) |
218 |
|
|
|
219 |
|
|
cph( |
220 |
|
|
phff = fupd |
221 |
|
|
cph) |
222 |
|
|
|
223 |
|
|
c----------------------------------------- |
224 |
|
|
c check warm start parameters |
225 |
|
|
c----------------------------------------- |
226 |
|
|
if (ifail .ne. 0) then |
227 |
|
|
if (iprint.ge.1) then |
228 |
|
|
print '(a)', ' ERROR : reading restart file' |
229 |
|
|
end if |
230 |
|
|
ifail = 2 |
231 |
|
|
goto 999 |
232 |
|
|
end if |
233 |
|
|
|
234 |
|
|
if (isize .ne. REAL_BYTE) then |
235 |
|
|
if (iprint.ge.1) then |
236 |
|
|
print '(a)', ' ERROR : uncompatible floating point format' |
237 |
|
|
end if |
238 |
|
|
ifail = 2 |
239 |
|
|
goto 999 |
240 |
|
|
end if |
241 |
|
|
|
242 |
|
|
if (mupd .lt. 1) then |
243 |
|
|
if (iprint .ge. 1) then |
244 |
|
|
print '(a)', ' ERROR : m is set too small in instore' |
245 |
|
|
endif |
246 |
|
|
ifail = 2 |
247 |
|
|
goto 999 |
248 |
|
|
endif |
249 |
|
|
|
250 |
|
|
c----------------------------------------- |
251 |
|
|
c cold start or warm restart ? |
252 |
|
|
c----------------------------------------- |
253 |
|
|
|
254 |
|
|
if (cold) then |
255 |
|
|
c--- start if cold start --- |
256 |
|
|
if (lphprint) then |
257 |
|
|
print '(a)', 'pathei-lsopt: cold start' |
258 |
|
|
end if |
259 |
|
|
|
260 |
heimbach |
1.4 |
print *, 'pathei-lsopt vor simul', nn |
261 |
|
|
print *, 'pathei-lsopt xx(1), gg(1) ', xx(1), gg(1) |
262 |
|
|
|
263 |
heimbach |
1.2 |
call simul( indic, nn, xx, ff, gg ) |
264 |
heimbach |
1.4 |
|
265 |
heimbach |
1.2 |
print *, 'pathei: nach simul: nn, ff = ', nn, ff |
266 |
heimbach |
1.4 |
print *, 'pathei: nach simul: xx(1), gg(1) = ', xx(1), gg(1) |
267 |
heimbach |
1.2 |
|
268 |
|
|
do i = 1, nn |
269 |
|
|
xdiff(i) = 1. |
270 |
|
|
end do |
271 |
|
|
|
272 |
|
|
cph( |
273 |
|
|
print *, 'pathei: vor dostore ' |
274 |
|
|
cph) |
275 |
|
|
call dostore( nn, xx, .true., 1 ) |
276 |
|
|
call dostore( nn, gg, .true., 2 ) |
277 |
|
|
call dostore( nn, xdiff, .true., 3 ) |
278 |
|
|
|
279 |
|
|
cph( |
280 |
|
|
print *, 'pathei: vor lswri ' |
281 |
|
|
cph) |
282 |
heimbach |
1.7 |
cph call lswri( isize, iiter, nn, xx, gg, lphprint ) |
283 |
heimbach |
1.2 |
|
284 |
|
|
cph( |
285 |
|
|
print *, 'pathei: vor gnorm0 ' |
286 |
|
|
cph) |
287 |
heimbach |
1.3 |
gnorm0 = DNRM2( nn, gg, 1 ) |
288 |
heimbach |
1.2 |
cph( |
289 |
|
|
print *, 'pathei: gnorm0 = ', gnorm0 |
290 |
|
|
cph) |
291 |
|
|
if (gnorm0 .lt. rmin) then |
292 |
|
|
ifail = 3 |
293 |
|
|
goto 1000 |
294 |
|
|
endif |
295 |
|
|
|
296 |
|
|
cph( |
297 |
|
|
phff = ff |
298 |
|
|
cph) |
299 |
|
|
|
300 |
|
|
if (lphprint) |
301 |
|
|
& print *, 'pathei-lsopt: cold; first call simul: ff = ', |
302 |
|
|
& phff |
303 |
|
|
|
304 |
|
|
c--- end if cold start --- |
305 |
|
|
else |
306 |
|
|
c--- start if warm start --- |
307 |
|
|
if (mm .ne. nn) then |
308 |
|
|
if (iprint.ge.1) then |
309 |
|
|
print '(a,i6)' |
310 |
|
|
$ , ' ERROR : inconsistent nn = ', mm |
311 |
|
|
endif |
312 |
|
|
ifail = 2 |
313 |
|
|
goto 999 |
314 |
|
|
endif |
315 |
|
|
if (mupd .ne. nupdate) then |
316 |
|
|
if (iprint.ge.1) then |
317 |
|
|
print '(a,i6)' |
318 |
|
|
$ , ' ERROR : inconsistent nupdate = ', mupd |
319 |
|
|
endif |
320 |
|
|
ifail = 2 |
321 |
|
|
goto 999 |
322 |
|
|
endif |
323 |
|
|
if (jmin .lt. 1) then |
324 |
|
|
if (iprint.ge.1) then |
325 |
|
|
print '(a,i6)' |
326 |
|
|
$ , ' ERROR : inconsistent jmin = ', jmin |
327 |
|
|
endif |
328 |
|
|
ifail = 2 |
329 |
|
|
goto 999 |
330 |
|
|
endif |
331 |
|
|
if (jmin .gt. mupd) then |
332 |
|
|
if (iprint.ge.1) then |
333 |
|
|
print '(a,i6)' |
334 |
|
|
$ , ' ERROR : inconsistent jmin = ', jmin |
335 |
|
|
endif |
336 |
|
|
ifail = 2 |
337 |
|
|
goto 999 |
338 |
|
|
endif |
339 |
|
|
if (jmax .gt. mupd) then |
340 |
|
|
if (iprint.ge.1) then |
341 |
|
|
print '(a,i6)' |
342 |
|
|
$ , ' ERROR : inconsistent jmax = ', jmax |
343 |
|
|
endif |
344 |
|
|
ifail = 2 |
345 |
|
|
goto 999 |
346 |
|
|
endif |
347 |
|
|
|
348 |
|
|
if (lphprint) then |
349 |
|
|
print *, 'pathei-lsopt: warm start; read via dostore' |
350 |
|
|
print * |
351 |
|
|
endif |
352 |
|
|
|
353 |
|
|
call dostore( nn, xx, .false., 1 ) |
354 |
|
|
call dostore( nn, gg, .false., 2 ) |
355 |
|
|
ff = fupd |
356 |
|
|
|
357 |
|
|
cph( |
358 |
|
|
phff = ff |
359 |
|
|
cph) |
360 |
|
|
|
361 |
|
|
if (lphprint) |
362 |
|
|
& print *, 'pathei-lsopt: warm; first dostore read: ff = ', |
363 |
|
|
& ff |
364 |
|
|
|
365 |
|
|
c--- end if warm start --- |
366 |
|
|
endif |
367 |
|
|
|
368 |
|
|
if (iprint .ge. 1) then |
369 |
|
|
print '(2a)', ' Itn Step Nfun Objective ' |
370 |
|
|
$ , 'Norm G Norm X Norm (X(k-1)-X(k))' |
371 |
|
|
end if |
372 |
|
|
|
373 |
|
|
c----------------------------------------- |
374 |
|
|
c print information line |
375 |
|
|
c----------------------------------------- |
376 |
|
|
if (cold) then |
377 |
|
|
print iform, iiter, tact, ifunc, ff, gnorm0 |
378 |
heimbach |
1.3 |
$ , DNRM2( nn, xx, 1 ), 0. |
379 |
heimbach |
1.2 |
|
380 |
mlosch |
1.5 |
Cml write(94,'(i5,2x,1pe11.4,4(2x,1pe8.1))') |
381 |
|
|
Cml & iiter, ff, gnorm0, tact, |
382 |
|
|
Cml & DNRM2( nn, xx, 1 ), 0. |
383 |
heimbach |
1.2 |
|
384 |
|
|
if ( itmax .EQ. 0 ) then |
385 |
|
|
ifail = 10 |
386 |
|
|
goto 1000 |
387 |
|
|
end if |
388 |
|
|
end if |
389 |
|
|
|
390 |
|
|
c======================================================================= |
391 |
|
|
c begin of loop |
392 |
|
|
c compute x(k+1) out of x(k) + t*d, where t > 0. |
393 |
|
|
c======================================================================= |
394 |
|
|
|
395 |
|
|
do iiter = 1, itmax |
396 |
|
|
|
397 |
|
|
if (lphprint) then |
398 |
|
|
print *, 'pathei-lsopt: ++++++++++++++++++++++++' |
399 |
|
|
print *, 'pathei-lsopt: entering iiter =', iiter |
400 |
|
|
end if |
401 |
|
|
|
402 |
|
|
c----------------------------------------- |
403 |
|
|
c store old values |
404 |
|
|
c----------------------------------------- |
405 |
|
|
do i = 1, nn |
406 |
|
|
gold(i) = gg(i) |
407 |
|
|
end do |
408 |
|
|
|
409 |
|
|
fold = ff |
410 |
|
|
cph( |
411 |
|
|
phniter0 = iiter |
412 |
|
|
phff = ff |
413 |
|
|
cph) |
414 |
|
|
|
415 |
|
|
c----------------------------------------- |
416 |
|
|
c compute new dd and xx |
417 |
|
|
c----------------------------------------- |
418 |
|
|
c |
419 |
|
|
call lsupdxx( |
420 |
|
|
& nn, ifail, lphprint |
421 |
|
|
& , jmin, jmax, nupdate |
422 |
|
|
& , ff, fmin, fold, gnorm0, dotdg |
423 |
|
|
& , gg, dd, xx, xdiff |
424 |
|
|
& , tmin, tmax, tact, epsx |
425 |
|
|
& ) |
426 |
|
|
|
427 |
|
|
c----------------------------------------- |
428 |
|
|
c check whether new direction is a descent one |
429 |
|
|
c----------------------------------------- |
430 |
|
|
if ( ifail .eq. 4) goto 1000 |
431 |
|
|
|
432 |
|
|
c----------------------------------------- |
433 |
|
|
c optline returns new values of x, f, and g |
434 |
|
|
c----------------------------------------- |
435 |
|
|
c |
436 |
|
|
call optline( |
437 |
|
|
& simul |
438 |
|
|
& , nn, ifail, lphprint |
439 |
|
|
& , ifunc, nfunc |
440 |
|
|
& , ff, dotdg |
441 |
|
|
& , tmin, tmax, tact, epsx |
442 |
|
|
& , dd, gg, xx, xdiff |
443 |
|
|
& ) |
444 |
|
|
|
445 |
|
|
if (lphprint) print *, 'pathei-lsopt: ', |
446 |
|
|
& 'back from optline; ifail = ', ifail |
447 |
|
|
if (lphprint) print *, 'pathei-lsopt: ', |
448 |
|
|
& 'dostore 1,2 at iiter ', iiter |
449 |
|
|
|
450 |
|
|
call dostore( nn, xx, .true., 1 ) |
451 |
|
|
call dostore( nn, gg, .true., 2 ) |
452 |
|
|
cph( |
453 |
heimbach |
1.7 |
cph call lswri( isize, iiter, nn, xx, gg, lphprint ) |
454 |
heimbach |
1.2 |
cph) |
455 |
|
|
|
456 |
heimbach |
1.3 |
gnorm = DNRM2( nn, gg, 1 ) |
457 |
heimbach |
1.2 |
|
458 |
|
|
c----------------------------------------- |
459 |
|
|
c print information line |
460 |
|
|
c----------------------------------------- |
461 |
|
|
print iform, iiter, tact, ifunc, ff, gnorm |
462 |
heimbach |
1.3 |
$ , DNRM2( nn, xx, 1 ), tact*DNRM2( nn, dd, 1 ) |
463 |
heimbach |
1.2 |
|
464 |
mlosch |
1.5 |
Cml write(94,'(i5,2x,1pe11.4,4(2x,1pe8.1))') |
465 |
|
|
Cml & iiter, ff, gnorm, tact, |
466 |
|
|
Cml & DNRM2( nn, xx, 1 ), tact*DNRM2( nn, dd, 1 ) |
467 |
heimbach |
1.2 |
|
468 |
|
|
c----------------------------------------- |
469 |
|
|
c check output mode of ifail |
470 |
|
|
c----------------------------------------- |
471 |
|
|
if (ifail .eq. 7 |
472 |
|
|
& .or. ifail .eq. 8 |
473 |
|
|
& .or. ifail .eq. 9) goto 1000 |
474 |
|
|
|
475 |
|
|
c----------------------------------------- |
476 |
|
|
c maximal number of simulation reached |
477 |
|
|
c no goto in order to update Hessian |
478 |
|
|
c----------------------------------------- |
479 |
|
|
if (ifail .eq. 6) ifail = 0 |
480 |
|
|
|
481 |
|
|
c----------------------------------------- |
482 |
|
|
c NOTE: stopping tests are now done |
483 |
|
|
c after having updated the matrix, |
484 |
|
|
c so that update information can be stored |
485 |
|
|
c in case of a later warm restart |
486 |
|
|
c----------------------------------------- |
487 |
|
|
c----------------------------------------- |
488 |
|
|
c compute new s, y |
489 |
|
|
c----------------------------------------- |
490 |
|
|
do i = 1, nn |
491 |
|
|
xdiff(i) = tact*dd(i) |
492 |
|
|
end do |
493 |
|
|
|
494 |
|
|
c----------------------------------------- |
495 |
|
|
c compute ys |
496 |
|
|
c----------------------------------------- |
497 |
|
|
do i = 1, nn |
498 |
|
|
gold(i) = gg(i)-gold(i) |
499 |
|
|
end do |
500 |
|
|
|
501 |
heimbach |
1.3 |
ys = DDOT( nn, gold, 1, xdiff, 1 ) |
502 |
heimbach |
1.2 |
if (ys .le. 0.) then |
503 |
|
|
ifail = 4 |
504 |
|
|
print *, 'pathei-lsopt: ys < 0; ifail = ', ifail |
505 |
|
|
goto 1000 |
506 |
|
|
endif |
507 |
|
|
cph( |
508 |
|
|
cph----------------------------------------- |
509 |
|
|
cph at this point it is clear that xdiff |
510 |
|
|
cph provides a true optimized solution; |
511 |
|
|
cph i.e. take new gradient gg to compute new dd |
512 |
|
|
cph----------------------------------------- |
513 |
|
|
cph) |
514 |
|
|
|
515 |
|
|
c----------------------------------------- |
516 |
|
|
c update pointers for hessupd |
517 |
|
|
c----------------------------------------- |
518 |
|
|
if (nupdate .gt. 0) then |
519 |
|
|
|
520 |
|
|
if (jmax .eq. 0) then |
521 |
|
|
jmax = jmax+1 |
522 |
|
|
if (lphprint) |
523 |
|
|
& print *, 'pathei-lsopt: ', |
524 |
|
|
& 'first pointer update after cold start; ', |
525 |
|
|
& 'iiter, jmin, jmax = ', iiter, jmin, jmax |
526 |
|
|
else |
527 |
|
|
jmax = jmax+1 |
528 |
|
|
if (jmax .gt. nupdate) jmax = jmax-nupdate |
529 |
|
|
|
530 |
|
|
if (jmin .eq. jmax) then |
531 |
|
|
if (lphprint) |
532 |
|
|
& print *, 'pathei-lsopt: pointers updated for ', |
533 |
|
|
& ' iiter, jmin, jmax = ', iiter, jmin, jmax |
534 |
|
|
jmin = jmin+1 |
535 |
|
|
if (jmin .gt. nupdate) jmin = jmin-nupdate |
536 |
|
|
end if |
537 |
|
|
end if |
538 |
|
|
|
539 |
|
|
c----------------------------------------- |
540 |
|
|
c compute sbar, ybar store rec = min 4,5 |
541 |
|
|
c----------------------------------------- |
542 |
|
|
r1 = sqrt(1./ys) |
543 |
heimbach |
1.3 |
call DSCAL( nn, r1, xdiff, 1 ) |
544 |
|
|
call DSCAL( nn, r1, gold, 1 ) |
545 |
heimbach |
1.2 |
|
546 |
|
|
if (lphprint) |
547 |
|
|
& print *, 'pathei-lsopt: dostore at iiter, jmin, jmax ', |
548 |
|
|
& iiter, jmin, jmax |
549 |
|
|
|
550 |
|
|
call dostore( nn, gold, .true., 2*jmax+2 ) |
551 |
|
|
call dostore( nn, xdiff, .true., 2*jmax+3 ) |
552 |
|
|
|
553 |
|
|
c----------------------------------------- |
554 |
|
|
c compute the diagonal preconditioner |
555 |
|
|
c use dd as temporary array |
556 |
|
|
c----------------------------------------- |
557 |
|
|
call dgscale( nn, gold, xdiff, dd, rmin ) |
558 |
|
|
|
559 |
|
|
endif |
560 |
|
|
|
561 |
|
|
c----------------------------------------- |
562 |
|
|
c test convergence and stopping criteria |
563 |
|
|
c----------------------------------------- |
564 |
|
|
eps1 = gnorm / gnorm0 |
565 |
|
|
|
566 |
|
|
if (eps1 .lt. epsg) then |
567 |
|
|
|
568 |
|
|
ifail = 11 |
569 |
|
|
goto 1000 |
570 |
|
|
endif |
571 |
|
|
|
572 |
|
|
c======================================================================= |
573 |
|
|
c return of loop |
574 |
|
|
c======================================================================= |
575 |
|
|
|
576 |
|
|
end do |
577 |
|
|
|
578 |
|
|
iiter = iiter - 1 |
579 |
|
|
ifail = 5 |
580 |
|
|
|
581 |
|
|
c----------------------------------------- |
582 |
|
|
c loop exit |
583 |
|
|
c----------------------------------------- |
584 |
|
|
1000 continue |
585 |
|
|
nfunc = ifunc |
586 |
|
|
epsg = eps1 |
587 |
|
|
|
588 |
|
|
c----------------------------------------------------------------------- |
589 |
|
|
c save data for warm start |
590 |
|
|
c----------------------------------------------------------------------- |
591 |
|
|
|
592 |
|
|
call outstore( nn, ff, gnorm0, nupdate, jmin, jmax ) |
593 |
|
|
|
594 |
|
|
c----------------------------------------------------------------------- |
595 |
|
|
c compute dd(i+1), xx(i+1) based on latest available gg(i), xx(i) |
596 |
|
|
c for offline version |
597 |
|
|
c----------------------------------------------------------------------- |
598 |
|
|
|
599 |
|
|
if (loffline) then |
600 |
|
|
|
601 |
|
|
call lsupdxx( |
602 |
|
|
& nn, ifail, lphprint |
603 |
|
|
& , jmin, jmax, nupdate |
604 |
|
|
& , ff, fmin, fold, gnorm0, dotdg |
605 |
|
|
& , gg, dd, xx, xdiff |
606 |
|
|
& , tmin, tmax, tact, epsx |
607 |
|
|
& ) |
608 |
|
|
|
609 |
|
|
c Save xx(i+1) to file for offline version. |
610 |
|
|
call optim_write_control( nn, xdiff ) |
611 |
|
|
|
612 |
|
|
end if |
613 |
|
|
|
614 |
|
|
c----------------------------------------- |
615 |
|
|
c print final information |
616 |
|
|
c----------------------------------------- |
617 |
|
|
if (iprint .ge. 5) then |
618 |
|
|
print * |
619 |
|
|
print '(a,i9)' |
620 |
|
|
$ , ' number of iterations..............', iiter |
621 |
|
|
print '(a,i9)' |
622 |
|
|
$ , ' number of simultations............', nfunc |
623 |
|
|
print '(a,e9.2)' |
624 |
|
|
$ , ' relative precision on g...........', epsg |
625 |
|
|
end if |
626 |
|
|
|
627 |
|
|
if (iprint.ge.10) then |
628 |
|
|
print * |
629 |
|
|
print '(a,e15.8)' |
630 |
|
|
$ , ' cost function...............', ff |
631 |
|
|
print '(a,e15.8)' |
632 |
heimbach |
1.3 |
$ , ' norm of x...................', DNRM2( nn, xx, 1 ) |
633 |
heimbach |
1.2 |
print '(a,e15.8)' |
634 |
heimbach |
1.3 |
$ , ' norm of g...................', DNRM2( nn, gg, 1 ) |
635 |
heimbach |
1.2 |
end if |
636 |
|
|
|
637 |
|
|
c----------------------------------------- |
638 |
|
|
c print error message |
639 |
|
|
c----------------------------------------- |
640 |
|
|
999 continue |
641 |
|
|
|
642 |
|
|
if (ifail .ne. 0) then |
643 |
|
|
if (iprint .ge. 5) then |
644 |
|
|
print * |
645 |
|
|
print '(a)', ' optimization stopped because :' |
646 |
|
|
if (ifail .lt. 0) then |
647 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
648 |
|
|
$ , ' user request during simul' |
649 |
|
|
else if (ifail .eq. 0) then |
650 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
651 |
|
|
$ , ' optimal solution found' |
652 |
|
|
else if (ifail .eq. 1) then |
653 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
654 |
|
|
$ , ' an input argument is wrong' |
655 |
|
|
else if (ifail .eq. 2) then |
656 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
657 |
|
|
$ , ' warm start file is corrupted' |
658 |
|
|
else if (ifail .eq. 3) then |
659 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
660 |
|
|
$ , ' the initial gradient is too small' |
661 |
|
|
else if (ifail .eq. 4) then |
662 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
663 |
|
|
$ , ' the search direction is not a descent one' |
664 |
|
|
else if (ifail .eq. 5) then |
665 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
666 |
|
|
$ , ' maximal number of iterations reached' |
667 |
|
|
else if (ifail .eq. 6) then |
668 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
669 |
|
|
$ , ' maximal number of simulations reached' |
670 |
|
|
else if (ifail .eq. 7) then |
671 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
672 |
|
|
$ , ' the linesearch failed' |
673 |
|
|
else if (ifail .eq. 8) then |
674 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
675 |
|
|
$ , ' the function could not be improved' |
676 |
|
|
else if (ifail .eq. 9) then |
677 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
678 |
|
|
$ , ' optline parameters wrong' |
679 |
|
|
else if (ifail .eq. 10) then |
680 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
681 |
|
|
$ , ' cold start, no optimization done' |
682 |
|
|
else if (ifail .eq. 11) then |
683 |
|
|
print '(2x,a8,I3,a)', 'ifail = ', ifail |
684 |
|
|
$ , ' convergence achieved within precision' |
685 |
|
|
end if |
686 |
|
|
print * |
687 |
|
|
else if (iprint .ge. 1) then |
688 |
|
|
print * |
689 |
|
|
print '(a,i9)' |
690 |
|
|
$ , ' after optimization ifail..........', ifail |
691 |
|
|
print * |
692 |
|
|
end if |
693 |
|
|
end if |
694 |
|
|
|
695 |
|
|
c----------------------------------------- |
696 |
|
|
c end of subroutine |
697 |
|
|
c----------------------------------------- |
698 |
|
|
return |
699 |
|
|
|
700 |
|
|
end |