1module minimization_mod
2 ! MODULE minimization_mod (prefix='min' category='1. High-level functionality')
3 !
4 !:Purpose: Minimization for variational assimilation, including the
5 ! subroutine that evaluates the cost function and its gradient.
6 !
7 use codePrecision_mod
8 use timeCoord_mod
9 use obsTimeInterp_mod
10 use verticalCoord_mod
11 use columnData_mod
12 use obsSpaceData_mod
13 use controlVector_mod
14 use midasMpi_mod
15 use horizontalCoord_mod
16 use gridStateVector_mod
17 use gridStateVectorFileIO_mod
18 use bmatrix_mod
19 use bMatrix1DVar_mod
20 use stateToColumn_mod
21 use varqc_mod
22 use rmatrix_mod
23 use costFunction_mod
24 use residual_mod
25 use obsOperators_mod
26 use quasinewton_mod
27 use utilities_mod
28 use biasCorrectionSat_mod
29 use columnVariableTransforms_mod
30 implicit none
31 save
32 private
33
34 ! public variables
35 public :: min_niter, min_nsim
36
37 ! public procedures
38 public :: min_Setup, min_minimize, min_writeHessian
39
40 type(struct_obs) , pointer :: obsSpaceData_ptr => null()
41 type(struct_columnData), pointer :: columnAnlInc_ptr => null()
42 type(struct_columnData), pointer :: columnTrlOnAnlIncLev_ptr => null()
43 type(struct_hco) , pointer :: hco_anl => null()
44
45 logical :: initialized = .false.
46
47 integer :: nmtra,nwork,min_nsim
48 integer :: nvadim_mpilocal ! for mpi
49 integer :: min_niter
50 integer,external :: get_max_rss
51 logical :: preconFileExists
52 character(len=20) :: preconFileName = './preconin'
53 character(len=20) :: preconFileNameOut = './pm1q'
54 character(len=20) :: preconFileNameOut_pert = './pm1q_pert'
55 integer :: n1gc = 3
56
57 ! variables stored for later call to min_writeHessian
58 real(8), allocatable :: vatra(:)
59 real(8), pointer :: controlVectorIncrSum_ptr(:)
60 real(8), allocatable :: controlVectorIncrSumZero(:)
61 real(8) :: zeps0, zdf1
62 integer :: itertot, isimtot, iztrl(5), imode
63 integer :: outerLoopIndex
64 integer :: numIterMaxInnerLoopUsed
65 logical :: llvazx
66 logical :: initializeForOuterLoop
67 logical :: deallocHessian
68 logical :: isMinimizationFinalCall
69 logical :: oneDVarMode
70
71 ! namelist variables
72 real(8) :: REPSG ! relative gradient amplitude used as stopping criteria
73 real(8) :: rdf1fac ! factor applied to initial cost function value to approximate final value
74 integer :: NVAMAJ ! number of vector pairs to store in memory for Hessian approximation
75 integer :: NITERMAX ! maximum number of minimization iterations
76 integer :: NSIMMAX ! maximum number of cost function evaluations during minimization
77 integer :: nwoqcv ! number of iterations to initially perform without varQC
78 logical :: lxbar ! generally not used any longer
79 logical :: lwrthess ! choose to write the Hessian approximation
80 logical :: lgrtest ! choose to perform the 'gradient test" before and after minimization
81 logical :: lvazx ! generally not used any longer
82 logical :: lvarqc ! choose to activate varQC
83
84 NAMELIST /NAMMIN/ NVAMAJ, NITERMAX, NSIMMAX
85 NAMELIST /NAMMIN/ LGRTEST
86 NAMELIST /NAMMIN/ lxbar, lwrthess, lvazx
87 NAMELIST /NAMMIN/ REPSG, rdf1fac
88 NAMELIST /NAMMIN/ LVARQC, NWOQCV
89
90CONTAINS
91
92 subroutine min_setup( nvadim_mpilocal_in, hco_anl_in, oneDVarMode_opt, &
93 varqc_opt, nwoqcv_opt )
94 !
95 !:Purpose: Reading NAMMIN namelist to setup variables for minimization.
96 !
97 implicit none
98
99 ! Arguments:
100 integer, intent(in) :: nvadim_mpilocal_in
101 type(struct_hco), pointer, intent(in) :: hco_anl_in
102 logical, optional, intent(in) :: oneDVarMode_opt
103 logical, optional, intent(out) :: varqc_opt
104 integer, optional, intent(out) :: nwoqcv_opt
105
106 ! Locals:
107 integer :: ierr,nulnam
108 integer,external :: fnom,fclos
109
110 call utl_tmg_start(90,'--Minimization')
111
112 if ( nvadim_mpilocal_in /= cvm_nvadim ) then
113 write(*,*) 'nvadim_mpilocal_in,cvm_nvadim=',nvadim_mpilocal_in,cvm_nvadim
114 call utl_abort('min_setup: control vector dimension not consistent')
115 endif
116
117 nvadim_mpilocal = nvadim_mpilocal_in
118
119 if ( present(oneDVarMode_opt) ) then
120 oneDVarMode = oneDVarMode_opt
121 else
122 oneDVarMode = .false.
123 end if
124
125 hco_anl => hco_anl_in
126
127 ! set default values for namelist variables
128 nvamaj = 6
129 nitermax = 0
130 rdf1fac = 0.25d0
131 nsimmax = 500
132 lgrtest = .false.
133 lwrthess = .true.
134 lxbar = .false.
135 lvazx = .false.
136 repsg = 1.0d-5
137 lvarqc = .false.
138 nwoqcv = 5
139
140 ! read in the namelist NAMMIN
141 nulnam=0
142 ierr=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
143 read(nulnam,nml=nammin,iostat=ierr)
144 if(ierr.ne.0) call utl_abort('min_setup: Error reading namelist')
145 write(*,nml=nammin)
146 ierr=fclos(nulnam)
147
148 IF(N1GC == 3)THEN
149 NMTRA = (4 + 2*NVAMAJ)*nvadim_mpilocal
150 ELSE
151 call utl_abort('min_setup: only N1GC=3 currently supported!')
152 END IF
153 WRITE(*,9401)N1GC,NVAMAJ,NMTRA
154 9401 FORMAT(4X,'N1GC = ',I2,4X,'NVAMAJ = ',I3,/5X,"NMTRA =",1X,I14)
155
156 if ( present(varqc_opt) ) varqc_opt = lvarqc
157 if ( present(nwoqcv_opt) ) nwoqcv_opt = nwoqcv
158
159 if(LVARQC .and. mmpi_myid == 0) write(*,*) 'VARIATIONAL QUALITY CONTROL ACTIVATED.'
160
161 initialized=.true.
162
163 call utl_tmg_stop(90)
164
165 end subroutine min_setup
166
167
168 subroutine min_minimize( outerLoopIndex_in, columnTrlOnAnlIncLev, obsSpaceData, controlVectorIncrSum, &
169 vazx, numIterMaxInnerLoop, deallocHessian_opt, &
170 isMinimizationFinalCall_opt, numIterMaxInnerLoopUsed_opt )
171 !
172 !:Purpose: Minimizing cost function to get the increments.
173 ! The maximum number of inner-loop iterations is set to nitermax if the
174 ! namelist variable nitermax is provided. Otherwise, it is set to the
175 ! numIterMaxInnerLoop supplied by the calling subroutine/program.
176 ! numIterMaxInnerLoopUsed_opt is passing the maximum number of inner-loop
177 ! iterations to the calling subroutine/program.
178 !
179 implicit none
180
181 ! Arguments:
182 integer, intent(in) :: outerLoopIndex_in
183 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev
184 type(struct_obs), intent(inout) :: obsSpaceData
185 real(8), target, intent(inout) :: controlVectorIncrSum(:)
186 real(8), intent(inout) :: vazx(:)
187 integer, intent(in) :: numIterMaxInnerLoop
188 integer, optional, intent(out) :: numIterMaxInnerLoopUsed_opt
189 logical, optional, intent(in) :: deallocHessian_opt
190 logical, optional, intent(in) :: isMinimizationFinalCall_opt
191
192 ! Locals:
193 type(struct_columnData) :: columnAnlInc
194 integer :: get_max_rss
195
196 call utl_tmg_start(90,'--Minimization')
197
198 write(*,*) '--------------------------------'
199 write(*,*) '--Starting subroutine minimize--'
200 write(*,*) '--------------------------------'
201 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
202
203 initializeForOuterLoop = .true.
204 outerLoopIndex = outerLoopIndex_in
205
206 if ( present(deallocHessian_opt) ) then
207 deallocHessian = deallocHessian_opt
208 else
209 deallocHessian = .true.
210 end if
211 if ( present(isMinimizationFinalCall_opt) ) then
212 isMinimizationFinalCall = isMinimizationFinalCall_opt
213 else
214 isMinimizationFinalCall = .true.
215 end if
216
217 if ( (nitermax > 0 .and. numIterMaxInnerLoop > 0) .or. &
218 (nitermax == 0 .and. numIterMaxInnerLoop == 0) ) then
219 call utl_abort('min_minimize: one of nitermax or numIterMaxInnerLoop should be zero and the other positive')
220 end if
221
222 if ( nitermax > 0 ) then
223 numIterMaxInnerLoopUsed = nitermax
224 else if ( numIterMaxInnerLoop > 0 ) then
225 numIterMaxInnerLoopUsed = numIterMaxInnerLoop
226 else
227 call utl_abort('min_minimize: one of the variables nIterMax and numIterMaxInnerLoop must be positive')
228 end if
229 if ( present(numIterMaxInnerLoopUsed_opt) ) numIterMaxInnerLoopUsed_opt = numIterMaxInnerLoopUsed
230
231 controlVectorIncrSum_ptr => controlVectorIncrSum
232
233 ! zero array for writing to hessian
234 allocate(controlVectorIncrSumZero(cvm_nvadim))
235 controlVectorIncrSumZero(:) = 0.0d0
236
237 call col_setVco(columnAnlInc,col_getVco(columnTrlOnAnlIncLev))
238 call col_allocate(columnAnlInc,col_getNumCol(columnTrlOnAnlIncLev),mpiLocal_opt=.true.)
239
240 write(*,*) 'oti_timeBinning: For 4D increment'
241 call oti_timeBinning(obsSpaceData,tim_nstepobsinc)
242
243 call quasiNewtonMinimization( columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData, vazx )
244
245 call col_deallocate(columnAnlInc)
246
247 deallocate(controlVectorIncrSumZero)
248
249 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
250 write(*,*) '--Done subroutine minimize--'
251
252 call utl_tmg_stop(90)
253
254 end subroutine min_minimize
255
256
257 subroutine quasiNewtonMinimization( columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData, vazx )
258 !
259 !:Purpose: 3D/En-VAR minimization
260 !
261 implicit none
262
263 ! Arguments:
264 type(struct_columnData), target, intent(in) :: columnAnlInc
265 type(struct_columnData), target, intent(in) :: columnTrlOnAnlIncLev
266 type(struct_obs), target, intent(inout) :: obsSpaceData
267 real(8) , intent(out) :: vazx(:)
268
269 ! Locals:
270 integer :: nulout = 6
271 integer :: impres
272 INTEGER :: NGRANGE = 10 ! range of powers of 10 used for gradient test
273 real :: zzsunused(1)
274 integer :: intUnused(1)
275 real(8),allocatable :: vazg(:)
276 real(8) :: dlds(1)
277 logical :: llvarqc, lrdvatra, llxbar
278 integer :: itermax, iterdone, itermaxtodo, isimmax, indic, iitnovqc
279 integer :: ierr, isimdone, jdata, isimnovqc
280 integer :: ibrpstamp, isim3d
281 real(8) :: zjsp, zxmin, zeps1
282 real(8) :: dlgnorm, dlxnorm
283 real(8) :: zeps0_000,zdf1_000
284 integer :: iterdone_000,isimdone_000
285
286 if (lvarqc .and. outerLoopIndex==1) call vqc_setup(obsSpaceData)
287
288 min_nsim=0
289
290 if(mmpi_myid == 0) then
291 impres=5
292 else
293 impres=0
294 endif
295
296 ! Check for preconditioning file
297 inquire(file=preconFileName,exist=preconFileExists)
298 if(preconFileExists) then
299 if(mmpi_myid == 0) write(*,*) 'PRECONDITIONING FILE FOUND:',preconFileName
300 else
301 if(mmpi_myid == 0) write(*,*) 'NO PRECONDITIONING FILE FOUND:',preconFileName
302 endif
303
304 ! Initialize Hessian only at first outerLoop (mpilocal)
305 if ( outerLoopIndex == 1 ) then
306 allocate(vatra(nmtra))
307 vatra(:) = 0.0d0
308 end if
309
310 allocate(vazg(nvadim_mpilocal),stat=ierr)
311 if(ierr.ne.0) then
312 write(*,*) 'minimization: Problem allocating memory! id=2',ierr
313 call utl_abort('min quasiNewtonMinimization')
314 endif
315
316 ! set module variable pointers for obsspacedata and the two column objects
317 obsSpaceData_ptr => obsSpaceData
318 columnAnlInc_ptr => columnAnlInc
319 columnTrlOnAnlIncLev_ptr => columnTrlOnAnlIncLev
320
321 ! Set-up the minimization
322
323 ! initialize iteration/simulation counters to zero
324 ITERTOT = 0
325 isimtot = 0
326
327 ! initialize control vector related arrays to zero
328 vazx(:) = 0.0d0
329 vazg(:) = 0.0d0
330
331
332 ! If minimization start without qcvar : turn off varqc to compute
333 ! innovations and test the gradients
334
335 if (outerLoopIndex == 1) zeps0 = repsg
336
337 if (preconFileExists) then
338 if ( mmpi_myid == 0 ) write(*,*) 'quasiNewtonMinimization : Preconditioning mode'
339 lrdvatra = .true.
340 imode = 2
341 llvazx = lvazx ! from namelist (default is .false.)
342 llxbar = lxbar ! from namelist (default is .false.)
343 else
344 lrdvatra = .false.
345 imode = 0
346 endif
347
348 ! read the hessian from preconin file at first outer-loop iteration
349 if ( lrdvatra .and. outerLoopIndex == 1 ) then
350 ibrpstamp = tim_getDatestamp() ! ibrpstamp is a I/O argument of hessianIO
351
352 call hessianIO (preconFileName,0, &
353 isim3d,ibrpstamp,zeps0_000,zdf1_000,iterdone_000, &
354 isimdone_000,iztrl,vatra,controlVectorIncrSumZero, &
355 vazx,llxbar,llvazx,n1gc,imode)
356 endif
357
358 iterdone = 0
359 isimdone = 0
360 itermax = numIterMaxInnerLoopUsed
361 itermaxtodo = itermax
362 isimmax = nsimmax
363
364 ! do the gradient test for the starting point of minimization
365 if ( lgrtest .and. outerLoopIndex == 1 ) then
366 ! save user-requested varqc switch
367 llvarqc = lvarqc
368
369 if ( nwoqcv > 0 ) lvarqc = .false.
370 call utl_tmg_start(91,'----QuasiNewton')
371 call grtest2(simvar,nvadim_mpilocal,vazx,ngrange)
372 call utl_tmg_stop(91)
373
374 lvarqc = llvarqc
375 endif
376
377 itertot = iterdone
378 isimtot = isimdone
379
380 llvarqc = lvarqc
381 if ( nwoqcv > 0 .and. outerLoopIndex == 1 ) lvarqc = .false.
382 INDIC =2
383 call utl_tmg_start(91,'----QuasiNewton')
384 call simvar(indic,nvadim_mpilocal,vazx,zjsp,vazg)
385 call utl_tmg_stop(91)
386 lvarqc = llvarqc
387
388 if ( outerLoopIndex == 1 ) zdf1 = rdf1fac * ABS(zjsp)
389
390 CALL PRSCAL(nvadim_mpilocal,VAZG,VAZG,DLGNORM)
391 DLGNORM = DSQRT(DLGNORM)
392 CALL PRSCAL(nvadim_mpilocal,VAZX,VAZX,DLXNORM)
393 DLXNORM = DSQRT(DLXNORM)
394 WRITE(*,*)' |X| = ', DLXNORM
395 WRITE(*,FMT=9220) ZJSP, DLGNORM
396 9220 FORMAT(/4X,'J(X) = ',G23.16,4X,'|Grad J(X)| = ',G23.16)
397
398 ! Iterations of the minimization algorithm
399
400 ZXMIN = epsilon(ZXMIN)
401 WRITE(*,FMT=9320)ZXMIN,ZDF1,ZEPS0,IMPRES,ITERMAX,NSIMMAX
402
403 9320 FORMAT(//,10X,' Minimization QNA_N1QN3 starts ...',/ &
404 10x,'DXMIN =',G23.16,2X,'DF1 =',G23.16,2X,'EPSG =',G23.16 &
405 /,10X,'IMPRES =',I3,2X,'NITER = ',I3,2X,'NSIM = ',I3)
406
407 ! Begin the minimization
408 if ( numIterMaxInnerLoopUsed > 0 ) then
409
410 ! First do iterations without var-QC only at the beginning of first outer-loop.
411 if (lvarqc .and. nwoqcv > 0 .and. iterdone < nwoqcv .and. outerLoopIndex == 1 ) then
412 iitnovqc = min(nwoqcv - iterdone,itermax)
413 isimnovqc = isimmax
414 lvarqc = .false.
415 call utl_tmg_start(91,'----QuasiNewton')
416
417 zeps1 = zeps0
418
419 call qna_n1qn3(simvar, dscalqn, dcanonb, dcanab, nvadim_mpilocal, vazx, &
420 zjsp,vazg, zxmin, zdf1, zeps1, impres, nulout, imode, &
421 iitnovqc, isimnovqc ,iztrl, vatra, nmtra, intUnused, &
422 zzsunused, dlds)
423 call utl_tmg_stop(91)
424 call fool_optimizer(obsSpaceData)
425
426 isimnovqc = isimnovqc - 1
427 itermaxtodo = itermaxtodo - iitnovqc + 1
428 isimmax = isimmax - isimnovqc + 1
429
430 itertot = itertot + iitnovqc
431 isimtot = isimtot + isimnovqc
432
433 zeps0 = zeps0/zeps1
434 lvarqc = .true.
435
436 if ((imode == 4 .or. imode == 1) .and. itertot < itermax) then
437 imode = 2
438 INDIC = 2
439 call utl_tmg_start(91,'----QuasiNewton')
440 call simvar(indic,nvadim_mpilocal,vazx,zjsp,vazg)
441 call utl_tmg_stop(91)
442 else
443 write(*,*) 'minimization_mod: qna_n1qn3 imode = ', imode
444 call utl_abort('minimization_mod: qna_n1qn3 mode not equal to 1 or 4')
445 endif
446 endif
447
448 ! Now do main minimization with var-QC
449 call utl_tmg_start(91,'----QuasiNewton')
450
451 if ( outerLoopIndex > 1 ) imode = 2
452
453 zeps1 = zeps0
454
455 call qna_n1qn3(simvar, dscalqn, dcanonb, dcanab, nvadim_mpilocal, vazx, &
456 zjsp,vazg, zxmin, zdf1, zeps1, impres, nulout, imode, &
457 itermaxtodo,isimmax, iztrl, vatra, nmtra, intUnused, zzsunused, &
458 dlds)
459 call utl_tmg_stop(91)
460 call fool_optimizer(obsSpaceData)
461
462 itertot = itertot + itermaxtodo
463 isimtot = isimtot + isimmax
464
465 zeps0 = zeps0/zeps1
466
467
468 WRITE(*,FMT=9500) imode,iterdone,itertot-iterdone,itertot,isimdone,isimtot-isimdone,isimtot
469 9500 FORMAT(//,20X,20('*'),2X &
470 ,/,20X,' Minimization ended with MODE:',I4 &
471 ,/,20X,' Number of iterations done in previous job:',I4 &
472 ,/,20X,' Number of iterations in this job:',I4 &
473 ,/,20X,' Total number of iterations:',I4 &
474 ,/,20X,'Number of simulations done in previous job:',I4 &
475 ,/,20X,' Number of simulations in this job:',I4 &
476 ,/,20X,' Total number of simulations:',I4)
477
478 min_niter = itertot
479
480 ! Test the gradient at the final point
481 if ( lgrtest .and. isMinimizationFinalCall ) then
482 WRITE(*,FMT=9400)
483 9400 FORMAT(//,12X,40('**'),/,12X,'TESTING THE GRADIENT AT THE FINAL POINT',/,40('**'))
484 call utl_tmg_start(91,'----QuasiNewton')
485 call grtest2(simvar,nvadim_mpilocal,vazx,ngrange)
486 call utl_tmg_stop(91)
487 end if
488
489 ! Print some contents of obsSpaceData to the listing
490 if(mmpi_myid == 0) then
491 do jdata = 1, min(1,obs_numHeader(obsSpaceData))
492 call obs_prnthdr(obsSpaceData,jdata)
493 call obs_prntbdy(obsSpaceData,jdata)
494 end do
495 endif
496
497 else
498
499 ! no analysis performed for ensemble mean, ensure mean increment is zero
500 vazx(:) = 0.0d0
501 min_niter = 0
502
503 endif ! if numIterMaxInnerLoopUsed > 0
504
505 ! deallocate the gradient
506 deallocate(vazg)
507 if ( deallocHessian .and. .not. lwrthess ) deallocate(vatra)
508
509 end subroutine quasiNewtonMinimization
510
511
512 subroutine min_writeHessian(vazx)
513 implicit none
514
515 ! Arguments:
516 real(8), intent(inout) :: vazx(:)
517
518 ! Locals:
519 integer :: dateStamp
520
521 call utl_tmg_start(90,'--Minimization')
522
523 if ( lwrthess ) then
524 ! Write out the Hessian to file
525 if ( mmpi_myid == 0 ) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
526
527 ! zero array for writing to hessian
528 if ( .not. allocated(controlVectorIncrSumZero) ) then
529 allocate(controlVectorIncrSumZero(cvm_nvadim))
530 end if
531 controlVectorIncrSumZero(:) = 0.0d0
532
533 dateStamp = tim_getDatestamp()
534 call hessianIO (preconFileNameOut,1, &
535 min_nsim,datestamp,zeps0,zdf1,itertot,isimtot, &
536 iztrl,vatra,controlVectorIncrSumZero,vazx,.true.,llvazx,n1gc,imode)
537
538 deallocate(controlVectorIncrSumZero)
539
540 if ( mmpi_myid == 0 ) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
541 endif
542
543 if ( lwrthess ) then
544 deallocate(vatra)
545 end if
546
547 call utl_tmg_stop(90)
548
549 end subroutine min_writeHessian
550
551
552 subroutine simvar(na_indic,na_dim,da_v,da_J,da_gradJ)
553 !
554 !:Purpose: To implement the Variational solver as described in
555 ! Courtier, 1997, Dual Formulation of Four-Dimentional Variational
556 ! Assimilation, Q.J.R., pp2449-2461.
557 !
558 !:Arguments:
559 ! :na_indic:
560 ! =1 No action taken
561 !
562 ! =2 Same as 4 (compute J and gradJ) but do not interrupt
563 ! timer of the minimizer.
564 !
565 ! =3 Compute Jo and gradJo only.
566 !
567 ! =4 Both J(u) and its gradient are computed.
568 !
569 ! .. Note:: 1 and 4 are reserved values for call back from
570 ! m1qn3. For direct calls, use a value other than 1
571 ! and 4.
572 implicit none
573
574 ! Arguments:
575 integer, intent(in) :: na_indic
576 integer, intent(in) :: na_dim ! Control-vector dimension, forecast-error covariance space
577 real(8), intent(in) :: da_v(na_dim) ! Control variable, forecast-error covariance space
578 real(8), intent(out) :: da_J ! Cost function of the Variational algorithm
579 real(8), intent(out) :: da_gradJ(na_dim) ! Gradient of the Variational Cost funtion
580
581 ! Locals:
582 real*8, dimension(na_dim) :: dl_v
583 real*8 :: dl_Jb, dl_Jo
584 type(struct_gsv), save :: statevector
585 type(struct_vco), pointer :: vco_anl
586
587 call utl_tmg_stop(91)
588 call utl_tmg_stop(90)
589
590 if (na_indic .ne. 1) then ! No action taken if na_indic == 1
591 min_nsim = min_nsim + 1
592
593 if(mmpi_myid == 0) then
594 write(*,*) 'Entering simvar for simulation ',min_nsim
595 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
596 endif
597
598 ! note: controlVectorIncrSum_ptr is sum of previous outer-loops
599 dl_v(1:nvadim_mpilocal) = da_v(1:nvadim_mpilocal) + controlVectorIncrSum_ptr(1:nvadim_mpilocal)
600
601 ! Computation of background term of cost function:
602 dl_Jb = dot_product(dl_v(1:nvadim_mpilocal),dl_v(1:nvadim_mpilocal))/2.d0
603 call mmpi_allreduce_sumreal8scalar(dl_Jb,"GRID")
604
605 if (oneDVarMode) then
606 call bmat1D_sqrtB(da_v, nvadim_mpilocal, columnAnlInc_ptr, obsSpaceData_ptr)
607 call cvt_transform(columnAnlInc_ptr, 'ZandP_tl', columnTrlOnAnlIncLev_ptr)
608 else
609 if (.not.gsv_isAllocated(statevector)) then
610 write(*,*) 'min-simvar: allocating increment stateVector'
611 vco_anl => col_getVco(columnTrlOnAnlIncLev_ptr)
612 call gsv_allocate(statevector, tim_nstepobsinc, hco_anl, vco_anl, &
613 dataKind_opt=pre_incrReal, mpi_local_opt=.true.)
614 call gio_readMaskFromFile(statevector,'./analysisgrid')
615 end if
616
617 call bmat_sqrtB(da_v,nvadim_mpilocal,statevector)
618
619 ! put in columnAnlInc H_horiz dx
620 call s2c_tl(statevector,columnAnlInc_ptr,columnTrlOnAnlIncLev_ptr,obsSpaceData_ptr)
621 end if
622
623 ! Save as OBS_WORK: H_vert H_horiz dx = Hdx
624 call utl_tmg_start(10,'--Observations')
625 call utl_tmg_start(18,'----ObsOper_TL')
626 call oop_Htl(columnAnlInc_ptr,columnTrlOnAnlIncLev_ptr,obsSpaceData_ptr,min_nsim, &
627 initializeLinearization_opt=initializeForOuterLoop)
628 call utl_tmg_stop(18)
629 call utl_tmg_stop(10)
630
631 ! Calculate OBS_OMA from OBS_WORK : d-Hdx
632 call res_compute(obsSpaceData_ptr)
633
634 call bcs_calcbias_tl(da_v,OBS_OMA,obsSpaceData_ptr,columnTrlOnAnlIncLev_ptr)
635
636 ! Save as OBS_WORK : R**-1/2 (d-Hdx)
637 call utl_tmg_start(10,'--Observations')
638 call rmat_RsqrtInverseAllObs(obsSpaceData_ptr,OBS_WORK,OBS_OMA)
639 call utl_tmg_stop(10)
640
641 ! Store J-obs in OBS_JOBS : 1/2 * R**-1 (d-Hdx)**2
642 call cfn_calcJo(obsSpaceData_ptr)
643
644 ! Store modified J_obs in OBS_JOBS : -ln((gamma-exp(J))/(gamma+1))
645 IF ( LVARQC ) THEN
646 call vqc_NlTl(obsSpaceData_ptr)
647 endif
648
649 dl_Jo = 0.d0
650 call utl_tmg_start(90,'--Minimization')
651 call utl_tmg_start(92,'----SumCostFunction')
652 call cfn_sumJo(obsSpaceData_ptr,dl_Jo)
653 da_J = dl_Jb + dl_Jo
654 if (na_indic == 3) then
655 da_J = dl_Jo
656 IF(mmpi_myid == 0) write(*,FMT='(6X,"SIMVAR: JO = ",G23.16,6X)') dl_Jo
657 else
658 da_J = dl_Jb + dl_Jo
659 IF(mmpi_myid == 0) write(*,FMT='(6X,"SIMVAR: Jb = ",G23.16,6X,"JO = ",G23.16,6X,"Jt = ",G23.16)') dl_Jb,dl_Jo,da_J
660 endif
661 call utl_tmg_stop(92)
662 call utl_tmg_stop(90)
663
664 ! Modify OBS_WORK : R**-1 (d-Hdx)
665 call utl_tmg_start(10,'--Observations')
666 call rmat_RsqrtInverseAllObs(obsSpaceData_ptr,OBS_WORK,OBS_WORK)
667 call utl_tmg_stop(10)
668
669 IF ( LVARQC ) THEN
670 call vqc_ad(obsSpaceData_ptr)
671 endif
672
673 ! Calculate adjoint of d-Hdx (mult OBS_WORK by -1)
674 call res_computeAd(obsSpaceData_ptr)
675
676 call col_zero(columnAnlInc_ptr)
677
678 ! Put in column : -H_vert**T R**-1 (d-Hdx)
679 call utl_tmg_start(10,'--Observations')
680 call utl_tmg_start(19,'----ObsOper_AD')
681 call oop_Had(columnAnlInc_ptr,columnTrlOnAnlIncLev_ptr,obsSpaceData_ptr, &
682 initializeLinearization_opt=initializeForOuterLoop)
683 call utl_tmg_stop(19)
684 call utl_tmg_stop(10)
685
686 ! Put in statevector -H_horiz**T H_vert**T R**-1 (d-Hdx)
687 if (oneDVarMode) then
688 ! no interpolation needed for 1Dvar case
689 else
690 call s2c_ad(statevector,columnAnlInc_ptr,columnTrlOnAnlIncLev_ptr,obsSpaceData_ptr)
691 end if
692
693 da_gradJ(:) = 0.d0
694 call bcs_calcbias_ad(da_gradJ,OBS_WORK,obsSpaceData_ptr)
695
696 if (oneDVarMOde) then
697 call cvt_transform( columnAnlInc_ptr, 'ZandP_ad', columnTrlOnAnlIncLev_ptr) ! IN
698 call bmat1D_sqrtBT(da_gradJ, nvadim_mpilocal, columnAnlInc_ptr, obsSpaceData_ptr)
699 else
700 call bmat_sqrtBT(da_gradJ,nvadim_mpilocal,statevector)
701 end if
702
703 if (na_indic .ne. 3) then
704 da_gradJ(1:nvadim_mpilocal) = dl_v(1:nvadim_mpilocal) + da_gradJ(1:nvadim_mpilocal)
705 endif
706 endif
707
708 call utl_tmg_start(90,'--Minimization')
709 call utl_tmg_start(91,'----QuasiNewton')
710
711 initializeForOuterLoop = .false.
712
713 if(mmpi_myid == 0) write(*,*) 'end of simvar'
714
715 end subroutine simvar
716
717
718 subroutine dscalqn(kdim,px,py,ddsc)
719 !:Purpose: Interface for the inner product to be used by the minimization
720 ! subroutines QNA_N1QN3.
721 !
722 !:Arguments:
723 ! i : kdim
724 ! i : px,py
725 ! o : ddsc
726 ! i : kzs
727 ! i : pzs
728 ! i : ddzs
729 implicit none
730
731 ! Arguments:
732 integer, intent(in) :: kdim ! dimension of the vectors
733 real(8), intent(in) :: px(kdim) ! vector for which <PX,PY> is being calculated
734 real(8), intent(in) :: py(kdim) ! vector for which <PX,PY> is being calculated
735 real(8), intent(out) :: ddsc ! result of the inner product
736
737 CALL PRSCAL(KDIM,PX,PY,DDSC)
738
739 end subroutine dscalqn
740
741
742 subroutine prscal(kdim,px,py,ddsc)
743 !
744 !:Purpose: To evaluate the inner product used in the minimization
745 !
746 implicit none
747
748 ! Arguments:
749 integer, intent(in) :: kdim ! dimension of the vectors
750 real(8), intent(in) :: px(kdim) ! vector for which <PX,PY> is being calculated
751 real(8), intent(in) :: py(kdim) ! vector for which <PX,PY> is being calculated
752 real(8), intent(out) :: ddsc ! result of the inner product
753
754 ! Locals:
755 integer :: j
756
757 ddsc = 0.D0
758
759 do j=1,nvadim_mpilocal
760 ddsc = ddsc + px(j)*py(j)
761 end do
762
763 call mmpi_allreduce_sumreal8scalar(ddsc,"GRID")
764
765 end subroutine prscal
766
767
768 subroutine dcanab(kdim,py,px)
769 !
770 !:Purpose: Change of variable associated with the canonical inner product:
771 !
772 ! * PX = L^-1 * Py with L related to the inner product
773 ! * <PX,PY> = PX^t L^t L PY
774 ! * (see the modulopt documentation aboutn DTCAB)
775 ! * NOTE: L is assumed to be the identity!
776 !
777 implicit none
778
779 ! Arguments:
780 integer, intent(in) :: kdim
781 real(8), intent(out) :: px(kdim)
782 real(8), intent(in) :: py(kdim)
783
784 ! Locals:
785 integer :: jdim
786
787 do jdim = 1, kdim
788 px(jdim) = py(jdim)
789 end do
790
791 end subroutine DCANAB
792
793
794 subroutine dcanonb(kdim,px,py)
795 !
796 !:Purpose: Change of variable associated with the canonical inner product:
797 !
798 ! * PY = L * PX with L related to the inner product
799 ! * <PX,PY> = PX^t L^t L PY
800 ! * (see the modulopt documentation about DTONB)
801 !
802 implicit none
803
804 ! Arguments:
805 integer, intent(in) :: kdim
806 real(8), intent(in) :: px(kdim)
807 real(8), intent(out) :: py(kdim)
808
809 ! Locals:
810 INTEGER :: jdim
811
812 do jdim = 1, kdim
813 py(jdim) = px(jdim)
814 end do
815
816 end subroutine DCANONB
817
818
819 subroutine hessianIO (cfname,status, &
820 nsim,kbrpstamp,zeps1,zdf1,itertot,isimtot, &
821 iztrl,vatra,vazxbar,vazx,llxbar,llvazx,k1gc,imode)
822 !
823 !:Purpose: Read-Write Hessian and increment (possibly for outer loop) on a
824 ! file
825 !
826 !:Arguments:
827 ! * i cfname
828 ! * i status
829 ! * i nsim
830 ! * io kbrpstamp
831 ! * i zeps1
832 ! * i zdf1
833 ! * i itertot
834 ! * i isimtot
835 ! * i iztrl
836 ! * i vatra
837 ! * i vazxbar
838 ! * i vazx
839 ! * i llxbar
840 ! * i llvazx
841 ! * i k1gc
842 ! * o imode
843 !
844 implicit none
845
846 ! Arguments:
847 character(len=*), intent(in) :: cfname ! precon file
848 integer , intent(in) :: status ! = 0 if READ, = 1 if WRITE
849 integer , intent(inout) :: nsim ! Number of simulations in QNA_N1QN3
850 integer , intent(inout) :: kbrpstamp ! Date
851 real(8) , intent(inout) :: zeps1 ! Parameter in QNA_N1QN3
852 real(8) , intent(inout) :: zdf1 ! Parameter in QNA_N1QN3
853 integer , intent(inout) :: itertot ! Parameter in QNA_N1QN3
854 integer , intent(inout) :: isimtot ! Parameter in QNA_N1QN3
855 integer, target , intent(inout) :: iztrl(5) ! Localisation parameters for Hessian
856 real(8), target , intent(inout) :: vatra(nmtra) ! Hessian
857 real(8), target , intent(inout) :: vazxbar(nvadim_mpilocal) ! Vazx of previous loop
858 real(8), target , intent(inout) :: vazx(nvadim_mpilocal) ! Current state of the minimization
859 logical , intent(in) :: llxbar ! read in vaxzbar if dates are compatible
860 logical , intent(in) :: llvazx ! Logical to read vazx
861 integer , intent(in) :: k1gc ! Minimizer ID (2: m1qn2, 3: m1qn3)
862 integer , intent(out) :: imode ! If status=0, set imode=0 (no prec) or 2 (prec)
863
864 ! Locals:
865 real(4), allocatable :: vatravec_r4_mpiglobal(:)
866 real(4), allocatable :: vatra_r4(:)
867 real(8), allocatable :: vazxbar_mpiglobal(:),vazx_mpiglobal(:)
868 integer :: ibrpstamp,ireslun, ierr, fnom, fclos
869 integer :: nvadim_mpiglobal,nmtra_mpiglobal
870 integer :: ivadim, itrunc
871 integer :: ivamaj
872 integer :: jvec, i1gc,ictrlvec,ii
873 integer, dimension(10), target, save :: iztrl_io
874 character(len=3) :: cl_version
875
876 if (status == 0) then
877 if (mmpi_myid == 0) write(*,*) 'Read Hessian in min_hessianIO from file ', cfname
878 call utl_tmg_start(93,'----ReadHess')
879 elseif (status == 1) then
880 if (mmpi_myid == 0) write(*,*) 'Write Hessian in min_hessianIO to file ', cfname
881 call utl_tmg_start(94,'----WriteHess')
882 else
883 call utl_abort('min_hessianIO: status not valid ')
884 endif
885
886 call rpn_comm_allreduce(nvadim_mpilocal,nvadim_mpiglobal,1,"mpi_integer","mpi_sum","GRID",ierr)
887 call rpn_comm_allreduce(nmtra ,nmtra_mpiglobal, 1,"mpi_integer","mpi_sum","GRID",ierr)
888
889 ireslun=0
890
891 if (status == 0) then
892 !
893 !- 1. Read Hessian on processor 0 and distribute the data to the other processors
894 !
895
896 if (mmpi_myid == 0) then
897 !- Open the Hessian matrix file
898 ierr = fnom(ireslun,cfname,'FTN+SEQ+UNF+OLD+R/O',0)
899
900 !- Checking version number
901 read (ireslun) cl_version, i1gc
902 if (trim(cl_version) == 'V5') then
903 write(*,*) 'min_hessianIO: using single precision V5 Hessian'
904 else if (trim(cl_version) == 'V4') then
905 write(*,*) 'min_hessianIO: using DEPRECIATED single precision V5 Hessian'
906 else
907 write(*,*)
908 write(*,*) 'min_hessianIO : Only V5 Hessian are supported, found ', trim(cl_version)
909 call utl_abort('min_hessianIO')
910 endif
911
912 if (i1gc == 3 .and. i1gc == k1gc) then
913 write(*,*) trim(cl_version),' M1QN3'
914 else
915 write(*,*) 'Version, n1gc =',trim(cl_version),i1gc
916 call utl_abort('min_hessianIO: Inconsistant input hessian')
917 endif
918
919 rewind (ireslun)
920
921 read(ireslun) cl_version,i1gc,nsim,ibrpstamp,zeps1,zdf1,itertot,isimtot,ivadim,itrunc
922 read(ireslun) ivamaj,iztrl_io
923 if ((ivamaj.ne.nvamaj).or.(nvadim_mpiglobal.ne.ivadim)) then
924 write(*,*) nvamaj,ivamaj,nvadim_mpiglobal,ivadim
925 call utl_abort('min_hessianIO : ERROR, size of V5 Hessian not consistent')
926 endif
927
928 end if
929
930 ! ibrpstamp and iztrl_io must be broadcasted
931 call rpn_comm_bcast(ibrpstamp, 1, "MPI_INTEGER", 0, "GRID", ierr)
932 call rpn_comm_bcast(iztrl_io , 10, "MPI_INTEGER", 0, "GRID", ierr)
933
934 !- Read the Hessian
935 if(mmpi_myid == 0) then
936 write(*,*) 'min_hessianIO : reading Hessian'
937 allocate(vatravec_r4_mpiglobal(nvadim_mpiglobal))
938 else
939 allocate(vatravec_r4_mpiglobal(1))
940 end if
941 allocate(vatra_r4(nvadim_mpilocal))
942
943 if (k1gc == 3) ictrlvec = 2*nvamaj+1
944 do jvec = 1, ictrlvec
945
946 if (mmpi_myid == 0) then
947 read(ireslun) vatravec_r4_mpiglobal
948 end if
949
950 call bmat_reduceToMPILocal_r4( vatra_r4, & ! OUT
951 vatravec_r4_mpiglobal) ! IN (contains data only on proc 0)
952 !$OMP PARALLEL DO PRIVATE(ii)
953 do ii = 1, nvadim_mpilocal
954 vatra((jvec-1)*nvadim_mpilocal+ii) = real(vatra_r4(ii),8)
955 enddo
956 !$OMP END PARALLEL DO
957
958 end do
959
960 deallocate(vatravec_r4_mpiglobal)
961 deallocate(vatra_r4)
962
963 imode = 2
964 iztrl(:) = iztrl_io(1:5)
965 if(k1gc == 3) then
966 iztrl(1) = nvadim_mpilocal
967 iztrl(2) = 0
968 iztrl(3) = nvamaj
969 iztrl(4) = iztrl_io(4)
970 iztrl(5) = iztrl_io(5)
971 end if
972
973 !- Read VAZXBAR
974 if (ibrpstamp == kbrpstamp .and. llxbar) then
975 if (mmpi_myid == 0) then
976 write(*,*) 'min_hessianIO : reading vazxbar'
977 allocate(vazxbar_mpiglobal(nvadim_mpiglobal))
978 read(ireslun) vazxbar_mpiglobal
979 else
980 allocate(vazxbar_mpiglobal(1))
981 end if
982 call bmat_reduceToMPILocal( vazxbar, & ! OUT
983 vazxbar_mpiglobal ) ! IN (contains data only on proc 0)
984 deallocate(vazxbar_mpiglobal)
985 end if
986
987 !- Read VAZX
988 if (ibrpstamp == kbrpstamp .and. llvazx) then
989 if (mmpi_myid == 0) then
990 write(*,*) 'min_hessianIO : reading vazx'
991 allocate(vazx_mpiglobal(nvadim_mpiglobal))
992 read(ireslun) vazx_mpiglobal
993 else
994 allocate(vazx_mpiglobal(1))
995 end if
996 call bmat_reduceToMPILocal( vazx, & ! OUT
997 vazx_mpiglobal ) ! IN (contains data only on proc 0)
998 deallocate(vazx_mpiglobal)
999 end if
1000
1001 if (ibrpstamp /= kbrpstamp) then
1002 kbrpstamp = ibrpstamp
1003 end if
1004
1005 if (mmpi_myid == 0) ierr = fclos(ireslun)
1006
1007 elseif (status == 1) then
1008 !
1009 !- 2. Write Hessian
1010 !
1011
1012 !- Open the Hessian matrix file on processor 0 and write some metadata
1013 if (mmpi_myid == 0) ierr = fnom(ireslun,cfname, 'FTN+SEQ+UNF' , 0)
1014 cl_version = 'V5'
1015 itrunc=0
1016 iztrl_io(1: 5) = iztrl(:)
1017 iztrl_io(6:10) = 0 ! dummy values for hessian size legacy purposes
1018 if(mmpi_myid == 0) write(ireslun) cl_version,k1gc,nsim,kbrpstamp,zeps1,zdf1,itertot,isimtot,nvadim_mpiglobal,itrunc
1019 if(mmpi_myid == 0) write(ireslun) nvamaj,iztrl_io
1020
1021 !- Write the Hessian
1022 if (mmpi_myid == 0) then
1023 allocate(vatravec_r4_mpiglobal(nvadim_mpiglobal))
1024 else
1025 allocate(vatravec_r4_mpiglobal(1))
1026 endif
1027 allocate(vatra_r4(nvadim_mpilocal))
1028
1029 if (k1gc == 3) ictrlvec = 2*nvamaj+1
1030 do jvec = 1, ictrlvec
1031 !$OMP PARALLEL DO PRIVATE(ii)
1032 do ii = 1, nvadim_mpilocal
1033 vatra_r4(ii) = vatra((jvec-1)*nvadim_mpilocal+ii)
1034 enddo
1035 !$OMP END PARALLEL DO
1036 call bmat_expandToMPIGlobal_r4( vatra_r4, & ! IN
1037 vatravec_r4_mpiglobal ) ! OUT
1038 if (mmpi_myid == 0) write(ireslun) vatravec_r4_mpiglobal
1039 enddo
1040 deallocate(vatravec_r4_mpiglobal)
1041 deallocate(vatra_r4)
1042
1043 !- Write VAZXBAR
1044 if (mmpi_myid == 0) then
1045 allocate(vazxbar_mpiglobal(nvadim_mpiglobal))
1046 else
1047 allocate(vazxbar_mpiglobal(1)) ! dummy
1048 end if
1049 call bmat_expandToMPIGlobal( vazxbar, & ! IN
1050 vazxbar_mpiglobal ) ! OUT
1051 if (mmpi_myid == 0) write(ireslun) vazxbar_mpiglobal(1:nvadim_mpiglobal)
1052 deallocate(vazxbar_mpiglobal)
1053
1054 !- Write VAZX
1055 if (mmpi_myid == 0) then
1056 allocate(vazx_mpiglobal(nvadim_mpiglobal))
1057 else
1058 allocate(vazx_mpiglobal(1))
1059 end if
1060 call bmat_expandToMPIGlobal( vazx, & ! IN
1061 vazx_mpiglobal ) ! OUT
1062 if (mmpi_myid == 0) write(ireslun) vazx_mpiglobal(1:nvadim_mpiglobal)
1063 deallocate(vazx_mpiglobal)
1064
1065 !- Close the Hessian matrix file
1066 if (mmpi_myid == 0) ierr = fclos(ireslun)
1067 else
1068 call utl_abort('min_hessianIO: status not valid')
1069 endif
1070
1071 if (status == 0) then
1072 call utl_tmg_stop(93)
1073 elseif (status == 1) then
1074 call utl_tmg_stop(94)
1075 endif
1076
1077 end subroutine hessianIO
1078
1079
1080 subroutine grtest2(simul,na_dim,da_x0,na_range)
1081 !
1082 !:Purpose: To compare the variation of the functional against what the
1083 ! gradient gives for small changes in the control variable. This test
1084 ! should be accurate for values as small as DLALPHA = SQRT(machine
1085 ! precision). (see Courtier, 1987)
1086 !
1087 !:Arguments:
1088 ! :na_range: the test will be carried over values of ALPHA ranging between
1089 ! 10**(-NA_RANGE) < ALPHA < 0.1
1090 !
1091 implicit none
1092
1093 ! Arguments:
1094 external :: simul ! simulator: return cost function estimate and its gradient
1095 integer, intent(in) :: na_dim ! Size of the control vector
1096 real(8), intent(in) :: da_x0(na_dim) ! Control vector
1097 integer, intent(in) :: na_range
1098
1099 ! Locals:
1100 integer :: nl_indic, nl_j
1101 real(8) :: dl_wrk(na_dim),dl_gradj0(na_dim), dl_x(na_dim)
1102 real(8) :: dl_J0, dl_J, dl_test, dl_start,dl_end
1103 real(8) :: dl_alpha, dl_gnorm0
1104
1105 ! 1. Initialize dl_gradj0 at da_x0
1106 ! ------------------------------------
1107
1108 nl_indic = 2
1109 call simul(nl_indic,na_dim,da_x0,dl_j0,dl_gradj0)
1110 dl_gnorm0 = dot_product(dl_gradj0,dl_gradj0)
1111 call mmpi_allreduce_sumreal8scalar(dl_gnorm0,"GRID")
1112 dl_start = 1.d0
1113 dl_end = 10.0d0**(-na_range)
1114 write(*,FMT=9100) dl_start,dl_end, dl_j0, dl_gnorm0
1115 ! 2. Perform the test
1116 ! ----------------
1117
1118 if(dl_gnorm0 == 0.d0)then
1119 write(*,FMT=9101)
1120 return
1121 end if
1122 write(*,FMT=9200)
1123 do nl_j = 1, na_range
1124 dl_alpha = 10.0d0**(- nl_j)
1125 dl_x(:) = da_x0(:) - dl_alpha*dl_gradJ0(:)
1126 call simul(nl_indic,na_dim,dl_x,dl_j,dl_wrk)
1127 dl_test = (dl_j-dl_j0)/(-dl_alpha * dl_gnorm0)
1128 write(*,FMT=9201)nl_j, dl_alpha, dl_j, dl_test
1129 end do
1130
11319100 format(//,4X,&
1132 'GRTEST- The gradient is being tested for',&
1133 G23.16,' <= ALPHA <= ',G23.16,/,12X,&
1134 'Initial value of J(X):',1x,G23.16,4x,&
1135 'Norm of GRAD J(X)**2: ',G23.16)
11369101 format(/,4X,'-In GRTEST: gradient vanishes exactly',&
1137 '. Gradient test cannot be performed at this point')
11389200 format(/,4X,'J',8X,'ALPHA',11X,'J(X)',12X,'TEST')
1139
11409201 format(2X,'GRTEST: step',2X,I3,4X,G23.16,4X,G23.16,4X,&
1141 G23.16)
1142
1143 end subroutine grtest2
1144
1145end module minimization_mod
