mlosch/optim_m1qn3/optim_sub.F

C $Header: /u/gcmpack/MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F,v 1.1 2012/04/26 11:10:06 mlosch Exp $
C $Name:  $

c     Include ECCO_CPPOPTIONS because the ecco_ctrl,cost files 
c     have headers with options for OBCS masks.
#include "ECCO_CPPOPTIONS.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

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

      if ( optimcycle .eq. 0 ) 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     ff has be read in optim_readparms, so we do not read it here again
       objf = ff
       df1  = objf-fmin
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--   read the model output into xx,adxx
      if ( indic .eq. 4 ) then
       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 nn, objf = ', nn, objf,
     &      xx(1), adxx(1)
      else
       print *, ' OPTIM_SUB: indic = ', indic, ' is not possible'
       stop 'ABNORMAL in S/R OPTIM_SUB'
      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: returned from m1qn3_offline'

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 *,' OPTIMS_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.2	C $Header: /u/gcmpack/MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F,v 1.1 2012/04/26 11:10:06 mlosch Exp $
2			C $Name: $
3	mlosch	1.1
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	mlosch	1.2	#ifdef DYNAMIC
52			_RL, dimension(:), allocatable :: xx, adxx
53	mlosch	1.1	#else
54			integer nmax
55			parameter( nmax = MAX_INDEPEND )
56			_RL xx(nmax)
57			_RL adxx(nmax)
58			#endif
59
60			c formal parameters of m1qn3
61			integer reverse
62			integer impres,imode(3),omode,niter,nsim,iz(5),indic
63			_RL dxmin,df1
64			character*3 normtype
65			c work arrays
66	mlosch	1.2	integer ndz
67	mlosch	1.1	CML _RL dz(ndz)
68			double precision, dimension(:), allocatable :: dz
69			c extra dummy variables
70			integer izs(1)
71			_RS rzs(1)
72			_RL dzs(1)
73			integer, parameter :: io = 60
74			character(), parameter :: fname_m1qn3 = 'm1qn3_output.txt'
75			c end of m1qn3 parameters
76
77			integer i
78
79			c == external ==
80
81			external simul_rc,euclid,ctonbe,ctcabe
82
83			c == end of interface ==
84
85			c-- Allocate memory for the control variables and the gradient vector.
86			#if defined(DYNAMIC)
87	mlosch	1.2	allocate( xx(nn) )
88			allocate( adxx(nn) )
89	mlosch	1.1	#endif
90
91	mlosch	1.2	#ifndef DYNAMIC
92	mlosch	1.1	if (nn .gt. nmax) then
93			print*,' OPTIMUM: Not enough space.'
94			print*,' nmax = ',nmax
95			print*,' nn = ',nn
96			print*
97			print*,' Set MAX_INDEPEND in Makefile .ge. ',nn
98			print*
99			stop ' ... stopped in OPTIMUM.'
100			endif
101			#endif
102
103			print*, ' OPTIM_SUB: Calling m1qn3_optim for iteration: ',
104			& optimcycle
105			print*, ' OPTIM_SUB: with nn, REAL_BYTE = ', nn, REAL_BYTE
106
107			c can be 'two','sup','dfn', see m1qn3 documentation for details
108			normtype='two'
109			c after reading data.optim some of these parameter values can be guessed
110			c impres=6, impres determines the amount of m1qn3-output see documentation
111			impres=iprint
112			c these should strictly be different (nsim>niter), but in practice
113			c it does not matter
114			niter = numiter
115			nsim = nfunc*niter
116			c epsg=1.d-8
117			dxmin=epsx
118			c will be set later
119			df1=-UNSET_RL
120			c
121			imode=(/0,1,0/)
122			omode=-1
123			c initialise work array
124			ndz = 3nn+nupdate(3*nn+1)
125			do i=1,5
126			iz(i)=0
127			enddo
128			allocate(dz(ndz))
129			do i=1,ndz
130			dz(i) = 0.
131			enddo
132			c these alway have to be set like this
133			reverse=1
134			indic=4
135			c initialise the dummy arguments that are not used
136			izs(1)=UNSET_I
137			rzs(1)=UNSET_RS
138			dzs(1)=UNSET_RL
139
140			if ( optimcycle .eq. 0 ) then
141			c-- cold start
142			print *, ' OPTIM_SUB: cold start, optimcycle =', optimcycle
143			imode(2) = 0
144			c this variable is the only one of the m1qn3-related common blocks
145			c that needs to be initialized here to make sure that we have a
146			c clean start
147			reentry = 0
148			c ff has be read in optim_readparms, so we do not read it here again
149			objf = ff
150			df1 = objf-fmin
151			c open output file for m1qn3
152			open(io,file=fname_m1qn3,status='unknown')
153			else
154			c-- warm restart
155			c requires restoring the state of m1qn3 with pickup file
156			print *, ' OPTIM_SUB: warm start, optimcycle =', optimcycle
157			imode(2) = 1
158			call optim_store_m1qn3(ndz,iz,dz,niter,nsim,epsg,df1,
159			I optimcycle,
160			I .false.)
161			c re-open output file for m1qn3
162			open(io,file=fname_m1qn3,status='old',position='append')
163			endif
164			c-- read the model output into xx,adxx
165			if ( indic .eq. 4 ) then
166			do i = 1,nn
167			xx(i) = 0.
168			adxx(i) = 0.
169			enddo
170			c
171			print *, ' OPTIM_SUB: read model state'
172			call optim_readdata( nn, ctrlname, .false., objf, xx )
173			call optim_readdata( nn, costname, .false., objf, adxx )
174			print *, ' OPTIM_SUB after reading nn, objf = ', nn, objf,
175			& xx(1), adxx(1)
176			else
177			print *, ' OPTIM_SUB: indic = ', indic, ' is not possible'
178			stop 'ABNORMAL in S/R OPTIM_SUB'
179			endif
180
181			c-- call the minimizer, a slightly modified version of m1qn3 v3.3
182			c (Gilbert & Lemarechal, 1989), downloaded in April 2012.
183			c simul_rc,euclid,ctonbe,ctcabe are external subroutines that
184			c are provided within m1qn3. euclid, ctonbe, ctcabe can be replaced
185			c by something more efficient, simul_rc is a dummy routine for
186			c the reverse communication mode and should not be changed.
187			print *, ' OPTIM_SUB: call m1qn3_offline'
188			call m1qn3_offline (simul_rc,euclid,ctonbe,ctcabe,
189			& nn,xx,objf,adxx,dxmin,df1,
190			& epsg,normtype,impres,io,imode,omode,niter,nsim,
191			& iz,dz,ndz,reverse,indic,izs,rzs,dzs)
192			close(io)
193			print *, ' OPTIM_SUB: returned from m1qn3_offline'
194
195			c write state of m1qn3 into pickup file for warm restart
196			call optim_store_m1qn3(ndz,iz,dz,niter,nsim,epsg,df1,
197			I optimcycle,
198			I .true.)
199			c write model control vector
200			print *,' OPTIMS_SUB: writing ', nn,' sized control to file ',
201			& ctrlname
202			c give the cost function a funny value to make sure that nobody
203			c mistakes it for the real one
204			call optim_writedata( nn, ctrlname, .false., -9999., xx )
205
206			c clean up
207	mlosch	1.2	#ifdef DYNAMIC
208			deallocate(xx, adxx)
209			#endif /* DYNAMIC */
210	mlosch	1.1	deallocate(dz)
211
212			c stopping criterion
213			if ( reverse .lt. 0 ) then
214			print *, ' OPTIM_SUB: reverse = ', reverse
215			print *, ' OPTIM_SUB: optimization terminated with omode = ',
216			& omode
217			stop 'ABNORMAL in S/R OPTIM_SUB'
218			endif
219
220			return
221			end