--- MITgcm/lsopt/lsopt_top.F 2002/02/05 20:34:34 1.1 +++ MITgcm/lsopt/lsopt_top.F 2002/11/15 04:03:24 1.2 @@ -0,0 +1,699 @@ + + subroutine lsopt_top( nn, xx, ff, gg, simul, optline + $ , epsx, fmin, epsg + $ , iprint, itmax, nfunc, nupdate + $ , dd, gold, xdiff + $ , loffline + $ , ifail ) + +c ================================================================== +c SUBROUTINE lsopt_top +c ================================================================== +c +c o uses a set of control variables, their adjoint variables, +c and a cost function value +c to compute an improved set of controls with respect to the +c cost function via a +c variable-storage Quasi-Newton method +c +c o Reference: J.C. Gilbert & C. Lemarechal +c Some numerical experiments with variable-storage +c quasi-Newton algorithms +c Mathematical Programming 45 (1989), pp. 407-435 +c +c o started: ??? not reproducible +c +c o changed: Patrick Heimbach, MIT/EAPS +c +c o Version: 2.0, 24-Feb-2000: +c (Version 1.0 is considered as version +c starting from which changes were made). +c - severe changes in structure including various +c bug-fixes to enable multi-optimization runs; +c - routine lsoptw incorporated into lsoptv +c - optimization iteration loop restructured +c - complete restructuring of handling +c cold start cases +c - number of 3 control flags for error handling +c (indic, moderl, ifail) reduced to 1 (ifail) +c and homogenized with routine lsline +c +c o Version: 2.1, 29-Feb-2000: +c - handling of case ifail = 6 changed; +c leads to stop of optimization +c (similar to case ifail = 8) +c - logical lphprint included +c +c ================================================================== +c SUBROUTINE lsopt_top +c ================================================================== + + implicit none + +#include + +c----------------------------------------- +c declare arguments +c----------------------------------------- + integer nn, iprint, itmax, nfunc, nupdate, ifail + + double precision xx(nn), ff, gg(nn), epsx, fmin, epsg + double precision dd(nn), gold(nn), xdiff(nn) + +cph( + integer phniter0, phniterold + double precision phff + COMMON /PH_OPTI/ phniter0, phniterold, phff +cph) + + external simul, optline + +c----------------------------------------- +C declare local variables +c----------------------------------------- + logical cold, lphprint, loffline + parameter (lphprint = .true.) + + integer mm, mupd, jmin, jmax, indic, isize, REAL_BYTE + integer i, iiter, ifunc + + double precision fupd + double precision r1, tmin, tmax, tact, gnorm, gnorm0, eps1 + double precision fold, ys + double precision dotdg + + external SDOT, SNRM2, SSCAL + double precision SDOT, SNRM2 + +c----------------------------------------- +c parameters +c----------------------------------------- + + double precision rmin + parameter( rmin = 1.e-20 ) + + character*(*) iform + parameter(iform='(i5,2x,1pe8.1,1x,i5,4x,1pe11.4,3(2x,1pe8.1))' ) + +c ================================================================== +c +c----------------------------------------- +c initialization +c----------------------------------------- + cold = .true. + fupd = 1.0e+10 + indic = 0 + tmin = 0. + tmax = 1.0e+10 + tact = 1.0 +cph( + phniterold = 0 +cph) + iiter = 0 + eps1 = 1.0 + ifunc = 0 + ifail = 0 + gnorm0 = 1. + +c----------------------------------------- +c initialization for dd and dds +c----------------------------------------- + + jmin = 1 + jmax = 0 + + mm = nn + mupd = nupdate + + REAL_BYTE = 8 + isize = REAL_BYTE + +c----------------------------------------- +c print information +c----------------------------------------- + if (iprint .ge. 1) then + print '(2x,a)', + $ '===============================================' + print '(2x,a)', + $ '=== O P T I M I Z A T I O N ===' + print '(2x,a)', + $ '===============================================' + print '(a,i9)' + $ , ' number of control variables.......', nn + print '(a,e9.2)' + $ , ' precision on x....................', epsx + print '(a,e9.2)' + $ , ' precision on g....................', epsg + print '(a,e9.2)' + $ , ' expected optimal function value...', fmin + print '(a,i9)' + $ , ' maximal number of iterations......', itmax + print '(a,i9)' + $ , ' maximal number of simulations.....', nfunc + print '(a,i9)' + $ , ' information level.................', iprint + print '(a,i9)' + $ , ' number of updates.................', nupdate + print '(a,i9)' + $ , ' size of used memory...............', 3*nn + endif + +c----------------------------------------- +c check arguments +c----------------------------------------- + + if (nn .le. 0) then + if (iprint.ge.1) then + print '(a,i6)' , ' ERROR : n = ', nn + endif + ifail = 1 + goto 999 + endif + + if (itmax .lt. 0) then + if (iprint.ge.1) then + print '(a,i6)' , ' ERROR : itmax = ', itmax + endif + ifail = 1 + goto 999 + endif + + if (nfunc .le. 0) then + if (iprint.ge.10) then + print '(a,i6)' , ' ERROR : nfunc = ', nfunc + endif + ifail = 1 + goto 999 + endif + + if (epsx .le. 0.) then + if (iprint.ge.1) then + print '(a,e9.2)', ' ERROR : epsx = ', epsx + endif + ifail = 1 + goto 999 + endif + + if (epsg .le. 0.) then + if (iprint.ge.1) then + print '(a,e9.2)', ' ERROR : epsg = ', epsg + endif + ifail = 1 + goto 999 + endif + + if (epsg .gt. 1.) then + if (iprint.ge.1) then + print '(a,e9.2)', ' ERROR : epsg = ', epsg + endif + ifail = 1 + goto 999 + endif + +cph( + print *, 'pathei: vor instore ' +cph) + call instore( mm, fupd, gnorm0, isize, mupd, jmin, jmax, cold, + & ifail ) + +cph( + phff = fupd +cph) + +c----------------------------------------- +c check warm start parameters +c----------------------------------------- + if (ifail .ne. 0) then + if (iprint.ge.1) then + print '(a)', ' ERROR : reading restart file' + end if + ifail = 2 + goto 999 + end if + + if (isize .ne. REAL_BYTE) then + if (iprint.ge.1) then + print '(a)', ' ERROR : uncompatible floating point format' + end if + ifail = 2 + goto 999 + end if + + if (mupd .lt. 1) then + if (iprint .ge. 1) then + print '(a)', ' ERROR : m is set too small in instore' + endif + ifail = 2 + goto 999 + endif + +c----------------------------------------- +c cold start or warm restart ? +c----------------------------------------- + + if (cold) then +c--- start if cold start --- + if (lphprint) then + print '(a)', 'pathei-lsopt: cold start' + print * + end if + + call simul( indic, nn, xx, ff, gg ) +cph( + print *, 'pathei: nach simul: nn, ff = ', nn, ff + print *, 'pathei: nach simul: xx, gg = ', xx(1), gg(1) +cph) + + do i = 1, nn + xdiff(i) = 1. + end do + +cph( + print *, 'pathei: vor dostore ' +cph) + call dostore( nn, xx, .true., 1 ) + call dostore( nn, gg, .true., 2 ) + call dostore( nn, xdiff, .true., 3 ) + +cph( + print *, 'pathei: vor lswri ' +cph) + call lswri( iiter, nn, xx, gg, lphprint ) + +cph( + print *, 'pathei: vor gnorm0 ' +cph) + gnorm0 = SNRM2( nn, gg, 1 ) +cph( + print *, 'pathei: gnorm0 = ', gnorm0 +cph) + if (gnorm0 .lt. rmin) then + ifail = 3 + goto 1000 + endif + +cph( + phff = ff +cph) + + if (lphprint) + & print *, 'pathei-lsopt: cold; first call simul: ff = ', + & phff + +c--- end if cold start --- + else +c--- start if warm start --- + if (mm .ne. nn) then + if (iprint.ge.1) then + print '(a,i6)' + $ , ' ERROR : inconsistent nn = ', mm + endif + ifail = 2 + goto 999 + endif + if (mupd .ne. nupdate) then + if (iprint.ge.1) then + print '(a,i6)' + $ , ' ERROR : inconsistent nupdate = ', mupd + endif + ifail = 2 + goto 999 + endif + if (jmin .lt. 1) then + if (iprint.ge.1) then + print '(a,i6)' + $ , ' ERROR : inconsistent jmin = ', jmin + endif + ifail = 2 + goto 999 + endif + if (jmin .gt. mupd) then + if (iprint.ge.1) then + print '(a,i6)' + $ , ' ERROR : inconsistent jmin = ', jmin + endif + ifail = 2 + goto 999 + endif + if (jmax .gt. mupd) then + if (iprint.ge.1) then + print '(a,i6)' + $ , ' ERROR : inconsistent jmax = ', jmax + endif + ifail = 2 + goto 999 + endif + + if (lphprint) then + print *, 'pathei-lsopt: warm start; read via dostore' + print * + endif + + call dostore( nn, xx, .false., 1 ) + call dostore( nn, gg, .false., 2 ) + ff = fupd + +cph( + phff = ff +cph) + + if (lphprint) + & print *, 'pathei-lsopt: warm; first dostore read: ff = ', + & ff + +c--- end if warm start --- + endif + + if (iprint .ge. 1) then + print '(2a)', ' Itn Step Nfun Objective ' + $ , 'Norm G Norm X Norm (X(k-1)-X(k))' + end if + +c----------------------------------------- +c print information line +c----------------------------------------- + if (cold) then + print iform, iiter, tact, ifunc, ff, gnorm0 + $ , SNRM2( nn, xx, 1 ), 0. + + write(94,'(i5,2x,1pe11.4,4(2x,1pe8.1))') + & iiter, ff, gnorm0, tact, + & SNRM2( nn, xx, 1 ), 0. + + if ( itmax .EQ. 0 ) then + ifail = 10 + goto 1000 + end if + end if + +c======================================================================= +c begin of loop +c compute x(k+1) out of x(k) + t*d, where t > 0. +c======================================================================= + + do iiter = 1, itmax + + if (lphprint) then + print *, 'pathei-lsopt: ++++++++++++++++++++++++' + print *, 'pathei-lsopt: entering iiter =', iiter + end if + +c----------------------------------------- +c store old values +c----------------------------------------- + do i = 1, nn + gold(i) = gg(i) + end do + + fold = ff +cph( + phniter0 = iiter + phff = ff +cph) + +c----------------------------------------- +c compute new dd and xx +c----------------------------------------- +c + call lsupdxx( + & nn, ifail, lphprint + & , jmin, jmax, nupdate + & , ff, fmin, fold, gnorm0, dotdg + & , gg, dd, xx, xdiff + & , tmin, tmax, tact, epsx + & ) + +c----------------------------------------- +c check whether new direction is a descent one +c----------------------------------------- + if ( ifail .eq. 4) goto 1000 + +c----------------------------------------- +c optline returns new values of x, f, and g +c----------------------------------------- +c + call optline( + & simul + & , nn, ifail, lphprint + & , ifunc, nfunc + & , ff, dotdg + & , tmin, tmax, tact, epsx + & , dd, gg, xx, xdiff + & ) + + if (lphprint) print *, 'pathei-lsopt: ', + & 'back from optline; ifail = ', ifail + if (lphprint) print *, 'pathei-lsopt: ', + & 'dostore 1,2 at iiter ', iiter + + call dostore( nn, xx, .true., 1 ) + call dostore( nn, gg, .true., 2 ) +cph( + call lswri( iiter, nn, xx, gg, lphprint ) +cph) + + gnorm = SNRM2( nn, gg, 1 ) + +c----------------------------------------- +c print information line +c----------------------------------------- + print iform, iiter, tact, ifunc, ff, gnorm + $ , SNRM2( nn, xx, 1 ), tact*SNRM2( nn, dd, 1 ) + + write(94,'(i5,2x,1pe11.4,4(2x,1pe8.1))') + & iiter, ff, gnorm, tact, + & SNRM2( nn, xx, 1 ), tact*SNRM2( nn, dd, 1 ) + +c----------------------------------------- +c check output mode of ifail +c----------------------------------------- + if (ifail .eq. 7 + & .or. ifail .eq. 8 + & .or. ifail .eq. 9) goto 1000 + +c----------------------------------------- +c maximal number of simulation reached +c no goto in order to update Hessian +c----------------------------------------- + if (ifail .eq. 6) ifail = 0 + +c----------------------------------------- +c NOTE: stopping tests are now done +c after having updated the matrix, +c so that update information can be stored +c in case of a later warm restart +c----------------------------------------- +c----------------------------------------- +c compute new s, y +c----------------------------------------- + do i = 1, nn + xdiff(i) = tact*dd(i) + end do + +c----------------------------------------- +c compute ys +c----------------------------------------- + do i = 1, nn + gold(i) = gg(i)-gold(i) + end do + + ys = SDOT( nn, gold, 1, xdiff, 1 ) + if (ys .le. 0.) then + ifail = 4 + print *, 'pathei-lsopt: ys < 0; ifail = ', ifail + goto 1000 + endif +cph( +cph----------------------------------------- +cph at this point it is clear that xdiff +cph provides a true optimized solution; +cph i.e. take new gradient gg to compute new dd +cph----------------------------------------- +cph) + +c----------------------------------------- +c update pointers for hessupd +c----------------------------------------- + if (nupdate .gt. 0) then + + if (jmax .eq. 0) then + jmax = jmax+1 + if (lphprint) + & print *, 'pathei-lsopt: ', + & 'first pointer update after cold start; ', + & 'iiter, jmin, jmax = ', iiter, jmin, jmax + else + jmax = jmax+1 + if (jmax .gt. nupdate) jmax = jmax-nupdate + + if (jmin .eq. jmax) then + if (lphprint) + & print *, 'pathei-lsopt: pointers updated for ', + & ' iiter, jmin, jmax = ', iiter, jmin, jmax + jmin = jmin+1 + if (jmin .gt. nupdate) jmin = jmin-nupdate + end if + end if + +c----------------------------------------- +c compute sbar, ybar store rec = min 4,5 +c----------------------------------------- + r1 = sqrt(1./ys) + call SSCAL( nn, r1, xdiff, 1 ) + call SSCAL( nn, r1, gold, 1 ) + + if (lphprint) + & print *, 'pathei-lsopt: dostore at iiter, jmin, jmax ', + & iiter, jmin, jmax + + call dostore( nn, gold, .true., 2*jmax+2 ) + call dostore( nn, xdiff, .true., 2*jmax+3 ) + +c----------------------------------------- +c compute the diagonal preconditioner +c use dd as temporary array +c----------------------------------------- + call dgscale( nn, gold, xdiff, dd, rmin ) + + endif + +c----------------------------------------- +c test convergence and stopping criteria +c----------------------------------------- + eps1 = gnorm / gnorm0 + + if (eps1 .lt. epsg) then + + ifail = 11 + goto 1000 + endif + +c======================================================================= +c return of loop +c======================================================================= + + end do + + iiter = iiter - 1 + ifail = 5 + +c----------------------------------------- +c loop exit +c----------------------------------------- + 1000 continue + nfunc = ifunc + epsg = eps1 + +c----------------------------------------------------------------------- +c save data for warm start +c----------------------------------------------------------------------- + + call outstore( nn, ff, gnorm0, nupdate, jmin, jmax ) + +c----------------------------------------------------------------------- +c compute dd(i+1), xx(i+1) based on latest available gg(i), xx(i) +c for offline version +c----------------------------------------------------------------------- + + if (loffline) then + + call lsupdxx( + & nn, ifail, lphprint + & , jmin, jmax, nupdate + & , ff, fmin, fold, gnorm0, dotdg + & , gg, dd, xx, xdiff + & , tmin, tmax, tact, epsx + & ) + +c Save xx(i+1) to file for offline version. + call optim_write_control( nn, xdiff ) + + end if + +c----------------------------------------- +c print final information +c----------------------------------------- + if (iprint .ge. 5) then + print * + print '(a,i9)' + $ , ' number of iterations..............', iiter + print '(a,i9)' + $ , ' number of simultations............', nfunc + print '(a,e9.2)' + $ , ' relative precision on g...........', epsg + end if + + if (iprint.ge.10) then + print * + print '(a,e15.8)' + $ , ' cost function...............', ff + print '(a,e15.8)' + $ , ' norm of x...................', SNRM2( nn, xx, 1 ) + print '(a,e15.8)' + $ , ' norm of g...................', SNRM2( nn, gg, 1 ) + end if + +c----------------------------------------- +c print error message +c----------------------------------------- + 999 continue + + if (ifail .ne. 0) then + if (iprint .ge. 5) then + print * + print '(a)', ' optimization stopped because :' + if (ifail .lt. 0) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' user request during simul' + else if (ifail .eq. 0) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' optimal solution found' + else if (ifail .eq. 1) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' an input argument is wrong' + else if (ifail .eq. 2) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' warm start file is corrupted' + else if (ifail .eq. 3) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' the initial gradient is too small' + else if (ifail .eq. 4) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' the search direction is not a descent one' + else if (ifail .eq. 5) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' maximal number of iterations reached' + else if (ifail .eq. 6) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' maximal number of simulations reached' + else if (ifail .eq. 7) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' the linesearch failed' + else if (ifail .eq. 8) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' the function could not be improved' + else if (ifail .eq. 9) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' optline parameters wrong' + else if (ifail .eq. 10) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' cold start, no optimization done' + else if (ifail .eq. 11) then + print '(2x,a8,I3,a)', 'ifail = ', ifail + $ , ' convergence achieved within precision' + end if + print * + else if (iprint .ge. 1) then + print * + print '(a,i9)' + $ , ' after optimization ifail..........', ifail + print * + end if + end if + +c----------------------------------------- +c end of subroutine +c----------------------------------------- + return + + end