/[MITgcm]/MITgcm/lsopt/lsopt_top.F
ViewVC logotype

Contents of /MITgcm/lsopt/lsopt_top.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (show annotations) (download)
Thu Sep 9 15:51:26 2004 UTC (18 years, 3 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint56b_post, checkpoint55, checkpoint56, checkpoint54f_post, checkpoint55i_post, checkpoint55c_post, checkpoint55g_post, checkpoint55d_post, checkpoint55d_pre, checkpoint55j_post, checkpoint55h_post, checkpoint55b_post, checkpoint55f_post, checkpoint56a_post, checkpoint55a_post, checkpoint55e_post
Changes since 1.3: +8 -7 lines
Small modifs and fixes
(mostly change to real*4 for large-scale ECCO runs)

1
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 ccc#include <blas1.h>
54
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 external DDOT, DNRM2, DSCAL
86 double precision DDOT, DNRM2
87
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 REAL_BYTE = 4
129 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 print *, 'pathei-lsopt vor simul', nn
261 print *, 'pathei-lsopt xx(1), gg(1) ', xx(1), gg(1)
262
263 call simul( indic, nn, xx, ff, gg )
264
265 print *, 'pathei: nach simul: nn, ff = ', nn, ff
266 print *, 'pathei: nach simul: xx(1), gg(1) = ', xx(1), gg(1)
267
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 call lswri( isize, iiter, nn, xx, gg, lphprint )
283
284 cph(
285 print *, 'pathei: vor gnorm0 '
286 cph)
287 gnorm0 = DNRM2( nn, gg, 1 )
288 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 $ , DNRM2( nn, xx, 1 ), 0.
379
380 write(94,'(i5,2x,1pe11.4,4(2x,1pe8.1))')
381 & iiter, ff, gnorm0, tact,
382 & DNRM2( nn, xx, 1 ), 0.
383
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 call lswri( isize, iiter, nn, xx, gg, lphprint )
454 cph)
455
456 gnorm = DNRM2( nn, gg, 1 )
457
458 c-----------------------------------------
459 c print information line
460 c-----------------------------------------
461 print iform, iiter, tact, ifunc, ff, gnorm
462 $ , DNRM2( nn, xx, 1 ), tact*DNRM2( nn, dd, 1 )
463
464 write(94,'(i5,2x,1pe11.4,4(2x,1pe8.1))')
465 & iiter, ff, gnorm, tact,
466 & DNRM2( nn, xx, 1 ), tact*DNRM2( nn, dd, 1 )
467
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 ys = DDOT( nn, gold, 1, xdiff, 1 )
502 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 call DSCAL( nn, r1, xdiff, 1 )
544 call DSCAL( nn, r1, gold, 1 )
545
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 $ , ' norm of x...................', DNRM2( nn, xx, 1 )
633 print '(a,e15.8)'
634 $ , ' norm of g...................', DNRM2( nn, gg, 1 )
635 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

  ViewVC Help
Powered by ViewVC 1.1.22