mlosch/optim_m1qn3/optim_sub.F

C $Header: /u/gcmpack/MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F,v 1.6 2012/06/04 12:23:34 mlosch Exp $
C $Name:  $

C     ECCO_CPPOPTIONS used to affect maxcvars and defined ALLOW_OBCS?_CONTROL
C#include "ECCO_CPPOPTIONS.h"
C     now:
C     CTRL_OPTIONS affects maxcvars and may define ALLOW_OBCS?_CONTROL
#include "CTRL_OPTIONS.h"

      subroutine optim_sub(
     I                 nn, ff
     &               )

c     ==================================================================
c     subroutine optim_sub
c     ==================================================================
c
c     o This is the main driver routine for the offline version of 
c       m1qn3 (m1qn3_offline). It
c       - sets all m1qn3 relevant parameters
c       - reads the model state (control vector, cost function, 
c         and gradient)
c       - reads the state of m1qn3_offline
c       - calls m1qn3_offline
c       - saves model state (control vector) and state of m1qn3_offline
c     The routine is somewhat lengthy and could be split into separate 
c     subroutines, but I felt that it was easier to write and test it
c     in this form
c     Martin Losch (Martin.Losch@awi.de), Apr, 2012
c
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"

#include "ctrl.h"
#include "optim.h"
#include "m1qn3_common.h"

c     == routine arguments ==

      integer nn
      _RL ff

c     == local variables ==

      _RL   objf

#ifdef DYNAMIC
      _RL, dimension(:), allocatable :: xx, adxx
#else
      integer nmax
      parameter( nmax = MAX_INDEPEND )
      _RL   xx(nmax)
      _RL   adxx(nmax)
#endif
      _RL  xxmean

CML      logical coldStart
c     formal parameters of m1qn3
      integer reverse
      integer impres,imode(3),omode,niter,nsim,iz(5),indic
      _RL dxmin,df1
      character*3 normtype
c     work arrays
      integer ndz
CML      _RL dz(ndz)
      double precision, dimension(:), allocatable :: dz
c     extra dummy variables
      integer izs(1)
      _RS rzs(1)
      _RL dzs(1)
      integer, parameter :: io = 60
      character*(*), parameter :: fname_m1qn3 = 'm1qn3_output.txt'
c     end of m1qn3 parameters

      integer i

c     == external ==

      external simul_rc,euclid,ctonbe,ctcabe

c     == end of interface ==

c--   Allocate memory for the control variables and the gradient vector.
#if defined(DYNAMIC)
      allocate(   xx(nn) )
      allocate( adxx(nn) )
#endif

#ifndef DYNAMIC
      if (nn .gt. nmax) then
        print*,' OPTIMUM: Not enough space.'
        print*,'          nmax = ',nmax
        print*,'            nn = ',nn
        print*
        print*,'          Set MAX_INDEPEND in Makefile .ge. ',nn
        print*
        stop   ' ... stopped in OPTIMUM.'
      endif
#endif

      print*, ' OPTIM_SUB: Calling m1qn3_optim for iteration: ',
     &     optimcycle
      print*, ' OPTIM_SUB: with nn, REAL_BYTE = ', nn, REAL_BYTE

c     can be 'two','sup','dfn', see m1qn3 documentation for details
      normtype='two'
c     after reading data.optim some of these parameter values can be guessed
c     impres=6, impres determines the amount of m1qn3-output see documentation
      impres=iprint
c     these should strictly be different (nsim>niter), but in practice
c     it does not matter
      niter   = numiter
      nsim    = nfunc*niter
c     epsg=1.d-8
      dxmin=epsx
c     will be set later
      df1=-UNSET_RL
c     
      imode=(/0,1,0/)
      omode=-1
c     initialise work array
      ndz = 3*nn+nupdate*(3*nn+1)
      do i=1,5
       iz(i)=0
      enddo
      allocate(dz(ndz))
      do i=1,ndz
       dz(i) = 0.
      enddo
c     these  alway have to be set like this
      reverse=1
      indic=4
c     initialise the dummy arguments that are not used
      izs(1)=UNSET_I
      rzs(1)=UNSET_RS
      dzs(1)=UNSET_RL

c--   first read the model output into xx, adxx, and cost function 
c     value into objf
      do i = 1,nn
       xx(i)   = 0.
       adxx(i) = 0.
      enddo
c
      print *, ' OPTIM_SUB: read model state'
      call optim_readdata( nn, ctrlname, .false., objf,   xx )
      call optim_readdata( nn, costname, .false., objf, adxx )
      print *, ' OPTIM_SUB after reading ', 
     &           ctrlname, ' and ', costname, ':'
      print *, ' OPTIM_SUB      nn = ', nn
      print *, ' OPTIM_SUB    objf = ', objf
      print *, ' OPTIM_SUB   xx(1) = ', xx(1)
      print *, ' OPTIM_SUB adxx(1) = ', adxx(1)
      
      if ( coldStart ) then
c--   cold start
       print *, ' OPTIM_SUB: cold start, optimcycle =', optimcycle
       imode(2) = 0
c     this variable is the only one of the m1qn3-related common blocks
c     that needs to be initialized here to make sure that we have a 
c     clean start
       reentry  = 0
c     compute expected decrease of cost function from objf and fmin;
c     this value is only used for a cold start of m1qn3_offline, for a
c     warm start df1 is overwritten with data from a restart file
       df1=objf-fmin
       if ( df1 .le. 0. ) then
        print *, ' OPTIM_SUB: df1 = objf-fmin = ', df1,
     &       ' but should be > 0.'
        stop 'ABNORMAL in S/R OPTIM_SUB'
       endif

c     open output file for m1qn3
       open(io,file=fname_m1qn3,status='unknown')
      else
c--   warm restart
c     requires restoring the state of m1qn3 with pickup file
       print *, ' OPTIM_SUB: warm start, optimcycle =', optimcycle
       imode(2) = 1
       call optim_store_m1qn3(ndz,iz,dz,niter,nsim,epsg,df1,
     I      optimcycle,
     I      .false.)
c     re-open output file for m1qn3
       open(io,file=fname_m1qn3,status='old',position='append')
      endif

c--   call the minimizer, a slightly modified version of m1qn3 v3.3
c     (Gilbert & Lemarechal, 1989), downloaded in April 2012.
c     simul_rc,euclid,ctonbe,ctcabe are external subroutines that
c     are provided within m1qn3. euclid, ctonbe, ctcabe can be replaced 
c     by something more efficient, simul_rc is a dummy routine for
c     the reverse communication mode and should not be changed.
      print *, ' OPTIM_SUB: call m1qn3_offline ........'
      call m1qn3_offline (simul_rc,euclid,ctonbe,ctcabe,
     &     nn,xx,objf,adxx,dxmin,df1,
     &     epsg,normtype,impres,io,imode,omode,niter,nsim,
     &     iz,dz,ndz,reverse,indic,izs,rzs,dzs)
      close(io)
      print *, ' OPTIM_SUB: ...........................'
      print *, ' OPTIM_SUB: returned from m1qn3_offline'
      print *, ' OPTIM_SUB:      nn = ', nn
      print *, ' OPTIM_SUB:   xx(1) = ', xx(1), xx(2)
      print *, ' OPTIM_SUB: adxx(1) = ', adxx(1), adxx(2)
      print *, ' OPTIM_SUB: omode   = ', omode
      print *, ' OPTIM_SUB: niter   = ', niter
      print *, ' OPTIM_SUB: nsim    = ', nsim 
      print *, ' OPTIM_SUB: reverse = ', reverse 

c     compute min/max/mean/std of output vector see if it is within 
c     reasonable bounds (prior to scaling)
      xxmean = sum(xx)/dble(nn)
      print *
      print *,' OPTIM_SUB: mean(xx) =', xxmean
      print *,' OPTIM_SUB:  max(xx) =', maxval(xx)
      print *,' OPTIM_SUB:  min(xx) =', minval(xx)
      print *,' OPTIM_SUB:  std(xx) =',sum((xx-xxmean)**2)/dble(nn)
      print *
c     write state of m1qn3 into pickup file for warm restart
      call optim_store_m1qn3(ndz,iz,dz,niter,nsim,epsg,df1,
     I     optimcycle,
     I     .true.)
c     write model control vector
      print *,' OPTIM_SUB: writing ', nn,' sized control to file ',
     &     ctrlname
c     give the cost function a funny value to make sure that nobody
c     mistakes it for the real one
      call optim_writedata( nn, ctrlname, .false., -9999., xx )

c     clean up
#ifdef DYNAMIC
      deallocate(xx, adxx)
#endif /* DYNAMIC */
      deallocate(dz)

c     stopping criterion
      if ( reverse .lt. 0 ) then
       print *, ' OPTIM_SUB: reverse = ', reverse
       print *, ' OPTIM_SUB: optimization terminated with omode = ', 
     &      omode
       stop 'ABNORMAL in S/R OPTIM_SUB'
      endif

      return
      end
1	mlosch	1.7	C $Header: /u/gcmpack/MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F,v 1.6 2012/06/04 12:23:34 mlosch Exp $
2	mlosch	1.2	C $Name: $
3	mlosch	1.1
4	mlosch	1.7	C ECCO_CPPOPTIONS used to affect maxcvars and defined ALLOW_OBCS?_CONTROL
5			C#include "ECCO_CPPOPTIONS.h"
6			C now:
7			C CTRL_OPTIONS affects maxcvars and may define ALLOW_OBCS?_CONTROL
8			#include "CTRL_OPTIONS.h"
9	mlosch	1.1
10			subroutine optim_sub(
11			I nn, ff
12			& )
13
14			c ==================================================================
15			c subroutine optim_sub
16			c ==================================================================
17			c
18			c o This is the main driver routine for the offline version of
19			c m1qn3 (m1qn3_offline). It
20			c - sets all m1qn3 relevant parameters
21			c - reads the model state (control vector, cost function,
22			c and gradient)
23			c - reads the state of m1qn3_offline
24			c - calls m1qn3_offline
25			c - saves model state (control vector) and state of m1qn3_offline
26			c The routine is somewhat lengthy and could be split into separate
27			c subroutines, but I felt that it was easier to write and test it
28			c in this form
29			c Martin Losch (Martin.Losch@awi.de), Apr, 2012
30			c
31			c ==================================================================
32
33			implicit none
34
35			c == global variables ==
36
37			#include "EEPARAMS.h"
38			#include "SIZE.h"
39
40			#include "ctrl.h"
41			#include "optim.h"
42			#include "m1qn3_common.h"
43
44			c == routine arguments ==
45
46			integer nn
47			_RL ff
48
49			c == local variables ==
50
51			_RL objf
52
53	mlosch	1.2	#ifdef DYNAMIC
54			_RL, dimension(:), allocatable :: xx, adxx
55	mlosch	1.1	#else
56			integer nmax
57			parameter( nmax = MAX_INDEPEND )
58			_RL xx(nmax)
59			_RL adxx(nmax)
60			#endif
61	mlosch	1.7	_RL xxmean
62	mlosch	1.1
63	mlosch	1.4	CML logical coldStart
64	mlosch	1.1	c formal parameters of m1qn3
65			integer reverse
66			integer impres,imode(3),omode,niter,nsim,iz(5),indic
67			_RL dxmin,df1
68			character*3 normtype
69			c work arrays
70	mlosch	1.2	integer ndz
71	mlosch	1.1	CML _RL dz(ndz)
72			double precision, dimension(:), allocatable :: dz
73			c extra dummy variables
74			integer izs(1)
75			_RS rzs(1)
76			_RL dzs(1)
77			integer, parameter :: io = 60
78			character(), parameter :: fname_m1qn3 = 'm1qn3_output.txt'
79			c end of m1qn3 parameters
80
81			integer i
82
83			c == external ==
84
85			external simul_rc,euclid,ctonbe,ctcabe
86
87			c == end of interface ==
88
89			c-- Allocate memory for the control variables and the gradient vector.
90			#if defined(DYNAMIC)
91	mlosch	1.2	allocate( xx(nn) )
92			allocate( adxx(nn) )
93	mlosch	1.1	#endif
94
95	mlosch	1.2	#ifndef DYNAMIC
96	mlosch	1.1	if (nn .gt. nmax) then
97			print*,' OPTIMUM: Not enough space.'
98			print*,' nmax = ',nmax
99			print*,' nn = ',nn
100			print*
101			print*,' Set MAX_INDEPEND in Makefile .ge. ',nn
102			print*
103			stop ' ... stopped in OPTIMUM.'
104			endif
105			#endif
106
107			print*, ' OPTIM_SUB: Calling m1qn3_optim for iteration: ',
108			& optimcycle
109			print*, ' OPTIM_SUB: with nn, REAL_BYTE = ', nn, REAL_BYTE
110
111			c can be 'two','sup','dfn', see m1qn3 documentation for details
112			normtype='two'
113			c after reading data.optim some of these parameter values can be guessed
114			c impres=6, impres determines the amount of m1qn3-output see documentation
115			impres=iprint
116			c these should strictly be different (nsim>niter), but in practice
117			c it does not matter
118			niter = numiter
119			nsim = nfunc*niter
120			c epsg=1.d-8
121			dxmin=epsx
122			c will be set later
123			df1=-UNSET_RL
124			c
125			imode=(/0,1,0/)
126			omode=-1
127			c initialise work array
128			ndz = 3nn+nupdate(3*nn+1)
129			do i=1,5
130			iz(i)=0
131			enddo
132			allocate(dz(ndz))
133			do i=1,ndz
134			dz(i) = 0.
135			enddo
136			c these alway have to be set like this
137			reverse=1
138			indic=4
139			c initialise the dummy arguments that are not used
140			izs(1)=UNSET_I
141			rzs(1)=UNSET_RS
142			dzs(1)=UNSET_RL
143
144	mlosch	1.4	c-- first read the model output into xx, adxx, and cost function
145			c value into objf
146			do i = 1,nn
147			xx(i) = 0.
148			adxx(i) = 0.
149			enddo
150			c
151			print *, ' OPTIM_SUB: read model state'
152			call optim_readdata( nn, ctrlname, .false., objf, xx )
153			call optim_readdata( nn, costname, .false., objf, adxx )
154			print *, ' OPTIM_SUB after reading ',
155			& ctrlname, ' and ', costname, ':'
156			print *, ' OPTIM_SUB nn = ', nn
157			print *, ' OPTIM_SUB objf = ', objf
158			print *, ' OPTIM_SUB xx(1) = ', xx(1)
159			print *, ' OPTIM_SUB adxx(1) = ', adxx(1)
160
161	mlosch	1.3	if ( coldStart ) then
162	mlosch	1.1	c-- cold start
163			print *, ' OPTIM_SUB: cold start, optimcycle =', optimcycle
164			imode(2) = 0
165			c this variable is the only one of the m1qn3-related common blocks
166			c that needs to be initialized here to make sure that we have a
167			c clean start
168			reentry = 0
169	mlosch	1.5	c compute expected decrease of cost function from objf and fmin;
170			c this value is only used for a cold start of m1qn3_offline, for a
171			c warm start df1 is overwritten with data from a restart file
172			df1=objf-fmin
173			if ( df1 .le. 0. ) then
174			print *, ' OPTIM_SUB: df1 = objf-fmin = ', df1,
175			& ' but should be > 0.'
176			stop 'ABNORMAL in S/R OPTIM_SUB'
177			endif
178
179	mlosch	1.1	c open output file for m1qn3
180			open(io,file=fname_m1qn3,status='unknown')
181			else
182			c-- warm restart
183			c requires restoring the state of m1qn3 with pickup file
184			print *, ' OPTIM_SUB: warm start, optimcycle =', optimcycle
185			imode(2) = 1
186			call optim_store_m1qn3(ndz,iz,dz,niter,nsim,epsg,df1,
187			I optimcycle,
188			I .false.)
189			c re-open output file for m1qn3
190			open(io,file=fname_m1qn3,status='old',position='append')
191			endif
192
193			c-- call the minimizer, a slightly modified version of m1qn3 v3.3
194			c (Gilbert & Lemarechal, 1989), downloaded in April 2012.
195			c simul_rc,euclid,ctonbe,ctcabe are external subroutines that
196			c are provided within m1qn3. euclid, ctonbe, ctcabe can be replaced
197			c by something more efficient, simul_rc is a dummy routine for
198			c the reverse communication mode and should not be changed.
199	mlosch	1.6	print *, ' OPTIM_SUB: call m1qn3_offline ........'
200	mlosch	1.1	call m1qn3_offline (simul_rc,euclid,ctonbe,ctcabe,
201			& nn,xx,objf,adxx,dxmin,df1,
202			& epsg,normtype,impres,io,imode,omode,niter,nsim,
203			& iz,dz,ndz,reverse,indic,izs,rzs,dzs)
204			close(io)
205	mlosch	1.6	print *, ' OPTIM_SUB: ...........................'
206	mlosch	1.1	print *, ' OPTIM_SUB: returned from m1qn3_offline'
207	mlosch	1.6	print *, ' OPTIM_SUB: nn = ', nn
208			print *, ' OPTIM_SUB: xx(1) = ', xx(1), xx(2)
209			print *, ' OPTIM_SUB: adxx(1) = ', adxx(1), adxx(2)
210			print *, ' OPTIM_SUB: omode = ', omode
211			print *, ' OPTIM_SUB: niter = ', niter
212			print *, ' OPTIM_SUB: nsim = ', nsim
213			print *, ' OPTIM_SUB: reverse = ', reverse
214	mlosch	1.1
215	mlosch	1.7	c compute min/max/mean/std of output vector see if it is within
216			c reasonable bounds (prior to scaling)
217			xxmean = sum(xx)/dble(nn)
218			print *
219			print *,' OPTIM_SUB: mean(xx) =', xxmean
220			print *,' OPTIM_SUB: max(xx) =', maxval(xx)
221			print *,' OPTIM_SUB: min(xx) =', minval(xx)
222			print ,' OPTIM_SUB: std(xx) =',sum((xx-xxmean)*2)/dble(nn)
223			print *
224	mlosch	1.1	c write state of m1qn3 into pickup file for warm restart
225			call optim_store_m1qn3(ndz,iz,dz,niter,nsim,epsg,df1,
226			I optimcycle,
227			I .true.)
228			c write model control vector
229	mlosch	1.7	print *,' OPTIM_SUB: writing ', nn,' sized control to file ',
230	mlosch	1.1	& ctrlname
231			c give the cost function a funny value to make sure that nobody
232			c mistakes it for the real one
233			call optim_writedata( nn, ctrlname, .false., -9999., xx )
234
235			c clean up
236	mlosch	1.2	#ifdef DYNAMIC
237			deallocate(xx, adxx)
238			#endif /* DYNAMIC */
239	mlosch	1.1	deallocate(dz)
240
241			c stopping criterion
242			if ( reverse .lt. 0 ) then
243			print *, ' OPTIM_SUB: reverse = ', reverse
244			print *, ' OPTIM_SUB: optimization terminated with omode = ',
245			& omode
246			stop 'ABNORMAL in S/R OPTIM_SUB'
247			endif
248
249			return
250			end