/[MITgcm]/MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F
ViewVC logotype

Contents of /MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F

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


Revision 1.1 - (show annotations) (download)
Thu Apr 26 11:10:06 2012 UTC (11 years, 11 months ago) by mlosch
Branch: MAIN
First working version of a new optimization package that uses a slightly
modified version of m1qn3, v3.3
(https://who.rocq.inria.fr/Jean-Charles.Gilbert/modulopt/optimization-routines/m1qn3/m1qn3.html)
to work as an offline optimizer. The advantage of m1qn3_offline is, that
it is run in reverse communication control mode, so that it gives back
control to the call routine (here a script) to provide a new estimate of the
cost function and the gradient based on the control vector. This way we can
do complete line searches that are meaningful.

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

  ViewVC Help
Powered by ViewVC 1.1.22