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

Diff of /MITgcm/lsopt/lsopt_top.F

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

revision 1.1 by heimbach, Tue Feb 5 20:34:34 2002 UTC revision 1.2 by heimbach, Fri Nov 15 04:03:24 2002 UTC
# Line 0  Line 1 
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    #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 SDOT, SNRM2, SSCAL
86          double precision     SDOT, SNRM2
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 = 8
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                print *
259             end if
260    
261             call simul( indic, nn, xx, ff, gg )
262    cph(
263             print *, 'pathei: nach simul:   nn, ff = ', nn, ff
264             print *, 'pathei: nach simul:   xx, gg = ', xx(1), gg(1)
265    cph)
266    
267             do i = 1, nn
268                xdiff(i) = 1.
269             end do
270    
271    cph(
272          print *, 'pathei: vor dostore '
273    cph)
274             call dostore( nn, xx, .true., 1 )
275             call dostore( nn, gg, .true., 2 )
276             call dostore( nn, xdiff, .true., 3 )
277                  
278    cph(
279          print *, 'pathei: vor lswri '
280    cph)
281             call lswri( iiter, nn, xx, gg, lphprint )
282    
283    cph(
284          print *, 'pathei: vor gnorm0 '
285    cph)
286             gnorm0 = SNRM2( nn, gg, 1 )
287    cph(
288          print *, 'pathei: gnorm0 = ', gnorm0
289    cph)
290             if (gnorm0 .lt. rmin) then
291                ifail = 3
292               goto 1000
293             endif
294    
295    cph(
296             phff = ff
297    cph)
298    
299             if (lphprint)
300         &        print *, 'pathei-lsopt: cold; first call simul: ff = ',
301         &        phff
302    
303    c--- end if cold start ---
304          else
305    c--- start if warm start ---
306             if (mm .ne. nn) then
307                if (iprint.ge.1) then
308                   print '(a,i6)'
309         $              , '  ERROR : inconsistent nn = ', mm
310                endif
311                ifail = 2
312                goto 999
313             endif
314             if (mupd .ne. nupdate) then
315                if (iprint.ge.1) then
316                   print '(a,i6)'
317         $              , '  ERROR : inconsistent nupdate = ', mupd
318                endif
319                ifail = 2
320                goto 999
321             endif
322             if (jmin .lt. 1) then
323                if (iprint.ge.1) then
324                   print '(a,i6)'
325         $              , '  ERROR : inconsistent jmin = ', jmin
326                endif
327                ifail = 2
328                goto 999
329             endif
330             if (jmin .gt. mupd) then
331                if (iprint.ge.1) then
332                   print '(a,i6)'
333         $              , '  ERROR : inconsistent jmin = ', jmin
334                endif
335                ifail = 2
336                goto 999
337             endif
338             if (jmax .gt. mupd) then
339                if (iprint.ge.1) then
340                   print '(a,i6)'
341         $              , '  ERROR : inconsistent jmax = ', jmax
342                endif
343                ifail = 2
344                goto 999
345             endif
346            
347             if (lphprint) then
348                print *, 'pathei-lsopt: warm start; read via dostore'
349                print *
350             endif
351    
352             call dostore( nn, xx, .false., 1 )
353             call dostore( nn, gg, .false., 2 )
354             ff = fupd
355    
356    cph(
357             phff = ff
358    cph)
359    
360             if (lphprint)
361         &        print *, 'pathei-lsopt: warm; first dostore read: ff = ',
362         &        ff
363    
364    c--- end if warm start ---
365          endif
366    
367          if (iprint .ge. 1) then
368             print '(2a)', '  Itn      Step  Nfun      Objective    '
369         $        , 'Norm G    Norm X   Norm (X(k-1)-X(k))'
370          end if
371    
372    c-----------------------------------------
373    c print information line
374    c-----------------------------------------
375             if (cold) then
376                print iform, iiter, tact, ifunc, ff, gnorm0
377         $           , SNRM2( nn, xx, 1 ), 0.
378    
379                write(94,'(i5,2x,1pe11.4,4(2x,1pe8.1))')
380         &           iiter, ff, gnorm0, tact,
381         &           SNRM2( nn, xx, 1 ), 0.
382    
383                if ( itmax .EQ. 0 ) then
384                   ifail = 10
385                   goto 1000
386                end if
387             end if
388    
389    c=======================================================================
390    c begin of loop
391    c compute x(k+1) out of x(k) + t*d, where  t > 0.
392    c=======================================================================
393    
394          do iiter = 1, itmax
395    
396             if (lphprint) then
397                print *, 'pathei-lsopt: ++++++++++++++++++++++++'
398                print *, 'pathei-lsopt: entering iiter =', iiter
399             end if
400    
401    c-----------------------------------------
402    c store old values
403    c-----------------------------------------
404          do i = 1, nn
405             gold(i) = gg(i)
406          end do
407    
408          fold = ff
409    cph(
410          phniter0 = iiter
411          phff = ff
412    cph)
413    
414    c-----------------------------------------
415    c compute new dd and xx
416    c-----------------------------------------
417    c
418          call lsupdxx(
419         &     nn, ifail, lphprint
420         &     , jmin, jmax, nupdate
421         &     , ff, fmin, fold, gnorm0, dotdg
422         &     , gg, dd, xx, xdiff
423         &     , tmin, tmax, tact, epsx
424         &     )
425    
426    c-----------------------------------------
427    c check whether new direction is a descent one
428    c-----------------------------------------
429          if ( ifail .eq. 4) goto 1000
430    
431    c-----------------------------------------
432    c optline returns new values of x, f, and g
433    c-----------------------------------------
434    c
435          call optline(
436         &     simul
437         &     , nn, ifail, lphprint
438         &     , ifunc, nfunc
439         &     , ff, dotdg
440         &     , tmin, tmax, tact, epsx
441         &     , dd, gg, xx, xdiff
442         &     )
443    
444          if (lphprint) print *, 'pathei-lsopt: ',
445         &     'back from optline; ifail = ', ifail
446          if (lphprint) print *, 'pathei-lsopt: ',
447         &     'dostore 1,2 at iiter ', iiter
448    
449          call dostore( nn, xx, .true., 1 )
450          call dostore( nn, gg, .true., 2 )
451    cph(
452          call lswri( iiter, nn, xx, gg, lphprint )
453    cph)
454    
455          gnorm = SNRM2( nn, gg, 1 )
456    
457    c-----------------------------------------
458    c print information line
459    c-----------------------------------------
460          print iform, iiter, tact, ifunc, ff, gnorm
461         $     , SNRM2( nn, xx, 1 ), tact*SNRM2( nn, dd, 1 )
462    
463          write(94,'(i5,2x,1pe11.4,4(2x,1pe8.1))')
464         &     iiter, ff, gnorm, tact,
465         &     SNRM2( nn, xx, 1 ), tact*SNRM2( nn, dd, 1 )
466    
467    c-----------------------------------------
468    c check output mode of ifail
469    c-----------------------------------------
470          if (ifail .eq. 7
471         &     .or. ifail .eq. 8
472         &     .or. ifail .eq. 9) goto 1000
473    
474    c-----------------------------------------
475    c maximal number of simulation reached
476    c no goto in order to update Hessian
477    c-----------------------------------------
478             if (ifail .eq. 6) ifail = 0
479    
480    c-----------------------------------------
481    c NOTE: stopping tests are now done
482    c after having updated the matrix,
483    c so that update information can be stored
484    c in case of a later warm restart
485    c-----------------------------------------
486    c-----------------------------------------
487    c  compute new s, y
488    c-----------------------------------------
489             do i = 1, nn
490                xdiff(i) = tact*dd(i)
491             end do
492    
493    c-----------------------------------------
494    c compute ys
495    c-----------------------------------------
496             do i = 1, nn
497                gold(i) = gg(i)-gold(i)
498             end do
499    
500             ys = SDOT( nn, gold, 1, xdiff, 1 )
501             if (ys .le. 0.) then
502                ifail = 4
503                print *, 'pathei-lsopt: ys < 0; ifail = ', ifail
504                goto 1000
505             endif
506    cph(
507    cph-----------------------------------------
508    cph at this point it is clear that xdiff
509    cph provides a true optimized solution;
510    cph i.e. take new gradient gg to compute new dd
511    cph-----------------------------------------
512    cph)
513    
514    c-----------------------------------------
515    c update pointers for hessupd
516    c-----------------------------------------
517          if (nupdate .gt. 0) then
518    
519             if (jmax .eq. 0) then
520                jmax = jmax+1
521                if (lphprint)
522         &           print *, 'pathei-lsopt: ',
523         &           'first pointer update after cold start; ',
524         &           'iiter, jmin, jmax = ', iiter, jmin, jmax
525             else
526                jmax = jmax+1
527                if (jmax .gt. nupdate) jmax = jmax-nupdate
528    
529                if (jmin .eq. jmax) then
530                   if (lphprint)
531         &              print *, 'pathei-lsopt: pointers updated for ',
532         &              ' iiter, jmin, jmax = ', iiter, jmin, jmax
533                   jmin = jmin+1
534                   if (jmin .gt. nupdate) jmin = jmin-nupdate
535                end if
536             end if
537    
538    c-----------------------------------------
539    c compute sbar, ybar store rec = min 4,5
540    c-----------------------------------------
541             r1 = sqrt(1./ys)
542             call SSCAL( nn, r1, xdiff, 1 )
543             call SSCAL( nn, r1, gold, 1 )
544    
545          if (lphprint)
546         &        print *, 'pathei-lsopt: dostore at iiter, jmin, jmax ',
547         &        iiter, jmin, jmax
548    
549             call dostore( nn, gold, .true., 2*jmax+2 )
550             call dostore( nn, xdiff, .true., 2*jmax+3 )
551    
552    c-----------------------------------------
553    c compute the diagonal preconditioner
554    c use dd as temporary array
555    c-----------------------------------------
556             call dgscale( nn, gold, xdiff, dd, rmin )
557    
558          endif
559    
560    c-----------------------------------------
561    c test convergence and stopping criteria
562    c-----------------------------------------
563          eps1 = gnorm / gnorm0
564    
565          if (eps1 .lt. epsg) then
566    
567             ifail = 11
568             goto 1000
569          endif
570    
571    c=======================================================================
572    c return of loop
573    c=======================================================================
574          
575          end do
576    
577          iiter = iiter - 1
578          ifail = 5
579    
580    c-----------------------------------------
581    c loop exit
582    c-----------------------------------------
583     1000 continue
584          nfunc = ifunc
585          epsg = eps1
586    
587    c-----------------------------------------------------------------------
588    c save data for warm start
589    c-----------------------------------------------------------------------
590    
591          call outstore( nn, ff, gnorm0, nupdate, jmin, jmax )
592    
593    c-----------------------------------------------------------------------
594    c compute dd(i+1), xx(i+1) based on latest available gg(i), xx(i)
595    c for offline version
596    c-----------------------------------------------------------------------
597    
598          if (loffline) then
599    
600             call lsupdxx(
601         &        nn, ifail, lphprint
602         &        , jmin, jmax, nupdate
603         &        , ff, fmin, fold, gnorm0, dotdg
604         &        , gg, dd, xx, xdiff
605         &        , tmin, tmax, tact, epsx
606         &        )
607    
608    c        Save xx(i+1) to file for offline version.
609             call optim_write_control( nn, xdiff )
610    
611          end if
612    
613    c-----------------------------------------
614    c print final information
615    c-----------------------------------------
616          if (iprint .ge. 5) then
617             print *
618             print '(a,i9)'
619         $        , '  number of iterations..............', iiter
620             print '(a,i9)'
621         $        , '  number of simultations............', nfunc
622             print '(a,e9.2)'
623         $        , '  relative precision on g...........', epsg
624          end if
625    
626          if (iprint.ge.10) then
627             print *
628             print '(a,e15.8)'
629         $        , '  cost function...............', ff
630             print '(a,e15.8)'
631         $        , '  norm of x...................', SNRM2( nn, xx, 1 )
632             print '(a,e15.8)'
633         $        , '  norm of g...................', SNRM2( nn, gg, 1 )
634          end if
635    
636    c-----------------------------------------
637    c print error message
638    c-----------------------------------------
639     999  continue
640    
641          if (ifail .ne. 0) then
642             if (iprint .ge. 5) then
643                print *
644                print '(a)', '  optimization stopped because :'
645                if (ifail .lt. 0) then
646                   print '(2x,a8,I3,a)', 'ifail = ', ifail
647         $              , '    user request during simul'
648                else if (ifail .eq. 0) then
649                   print '(2x,a8,I3,a)', 'ifail = ', ifail
650         $              , '    optimal solution found'
651                else if (ifail .eq. 1) then
652                   print '(2x,a8,I3,a)', 'ifail = ', ifail
653         $              , '    an input argument is wrong'
654                else if (ifail .eq. 2) then
655                   print '(2x,a8,I3,a)', 'ifail = ', ifail
656         $              , '    warm start file is corrupted'
657                else if (ifail .eq. 3) then
658                   print '(2x,a8,I3,a)', 'ifail = ', ifail
659         $              , '    the initial gradient is too small'
660                else if (ifail .eq. 4) then
661                   print '(2x,a8,I3,a)', 'ifail = ', ifail
662         $              , '    the search direction is not a descent one'
663                else if (ifail .eq. 5) then
664                   print '(2x,a8,I3,a)', 'ifail = ', ifail
665         $              , '    maximal number of iterations reached'
666                else if (ifail .eq. 6) then
667                   print '(2x,a8,I3,a)', 'ifail = ', ifail
668         $              , '    maximal number of simulations reached'
669                else if (ifail .eq. 7) then
670                   print '(2x,a8,I3,a)', 'ifail = ', ifail
671         $              , '    the linesearch failed'
672                else if (ifail .eq. 8) then
673                   print '(2x,a8,I3,a)', 'ifail = ', ifail
674         $              , '    the function could not be improved'
675                else if (ifail .eq. 9) then
676                   print '(2x,a8,I3,a)', 'ifail = ', ifail
677         $              , '    optline parameters wrong'
678                else if (ifail .eq. 10) then
679                   print '(2x,a8,I3,a)', 'ifail = ', ifail
680         $              , '    cold start, no optimization done'
681                else if (ifail .eq. 11) then
682                   print '(2x,a8,I3,a)', 'ifail = ', ifail
683         $              , '    convergence achieved within precision'
684                end if
685                print *
686             else if (iprint .ge. 1) then
687                print *
688                print '(a,i9)'
689         $           , '  after optimization ifail..........', ifail
690                print *
691             end if
692          end if
693    
694    c-----------------------------------------
695    c end of subroutine
696    c-----------------------------------------
697          return
698    
699          end

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

  ViewVC Help
Powered by ViewVC 1.1.22