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	C $Header: /u/gcmpack/MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F,v 1.6 2012/06/04 12:23:34 mlosch Exp $
2	C $Name: $
3
4	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
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	#ifdef DYNAMIC
54	_RL, dimension(:), allocatable :: xx, adxx
55	#else
56	integer nmax
57	parameter( nmax = MAX_INDEPEND )
58	_RL xx(nmax)
59	_RL adxx(nmax)
60	#endif
61	_RL xxmean
62
63	CML logical coldStart
64	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	integer ndz
71	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	allocate( xx(nn) )
92	allocate( adxx(nn) )
93	#endif
94
95	#ifndef DYNAMIC
96	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	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	if ( coldStart ) then
162	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	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	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	print *, ' OPTIM_SUB: call m1qn3_offline ........'
200	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	print *, ' OPTIM_SUB: ...........................'
206	print *, ' OPTIM_SUB: returned from m1qn3_offline'
207	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
215	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	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	print *,' OPTIM_SUB: writing ', nn,' sized control to file ',
230	& 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	#ifdef DYNAMIC
237	deallocate(xx, adxx)
238	#endif /* DYNAMIC */
239	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