mlosch/optim_m1qn3/optim_sub.F

C $Header: $
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

#if defined (DYNAMIC)
      _RL   xx(nn)
      _RL   adxx(nn)
#elif defined (USE_POINTER) || (MAX_INDEPEND == 0)
      _RL   xx
      _RL   adxx
      pointer (pxx,xx(1))
      pointer (padxx,adxx(1))
#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, mupdate
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)
#elif defined(USE_POINTER) || (MAX_INDEPEND == 0)
      call myalloc( pxx  ,  nn*REAL_BYTE )
      call myalloc( padxx,  nn*REAL_BYTE )
#endif

#if defined (DYNAMIC)
#elif defined(USE_POINTER) || (MAX_INDEPEND == 0)
#else
      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
      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.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 = 3nn+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