! Copyright (c) 2004-2024 Lars Nerger ! ! This file is part of PDAF. ! ! PDAF is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License ! as published by the Free Software Foundation, either version ! 3 of the License, or (at your option) any later version. ! ! PDAF is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public ! License along with PDAF. If not, see <http://www.gnu.org/licenses/>. ! !$Id$ !BOP ! ! !ROUTINE: PDAF_hyb3dvar_optim_lbfgs --- Optimization with LBFGS for Hyb3dVar ! ! !INTERFACE: SUBROUTINE PDAF_hyb3dvar_optim_lbfgs(step, dim_p, dim_ens, dim_cv_par_p, dim_cv_ens_p, & dim_obs_p, ens_p, obs_p, dy_p, v_par_p, v_ens_p, U_prodRinvA, & U_cvt, U_cvt_adj, U_cvt_ens, U_cvt_adj_ens, U_obs_op_lin, U_obs_op_adj, & opt_parallel, beta_3dvar, screen) ! !DESCRIPTION: ! Optimization routine for ensemble 3D-Var using the LBFGS solver ! ! Variant for domain decomposed states. ! ! ! This is a core routine of PDAF and ! should not be changed by the user ! ! ! !REVISION HISTORY: ! 2021-03 - Lars Nerger - Initial code ! Later revisions - see svn log ! ! !USES: USE PDAF_timer, & ONLY: PDAF_timeit USE PDAF_memcounting, & ONLY: PDAF_memcount USE PDAF_mod_filtermpi, & ONLY: mype USE PDAF_mod_filter, & ONLY: m_lbfgs_var, factr_lbfgs_var, pgtol_lbfgs_var, debug IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: step ! Current time step INTEGER, INTENT(in) :: dim_p ! PE-local state dimension INTEGER, INTENT(in) :: dim_ens ! ensemble size INTEGER, INTENT(in) :: dim_cv_par_p ! Size of control vector (parameterized) INTEGER, INTENT(in) :: dim_cv_ens_p ! Size of control vector (ensemble) INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector REAL, INTENT(in) :: ens_p(dim_p, dim_ens) ! PE-local state ensemble REAL, INTENT(in) :: obs_p(dim_obs_p) ! Vector of observations REAL, INTENT(in) :: dy_p(dim_obs_p) ! Background innovation REAL, INTENT(inout) :: v_par_p(dim_cv_par_p) ! Control vector (parameterized part) REAL, INTENT(inout) :: v_ens_p(dim_cv_ens_p) ! Control vector (ensemble part) INTEGER, INTENT(in) :: screen ! Verbosity flag INTEGER, INTENT(in) :: opt_parallel ! Whether to use a decomposed control vector REAL, INTENT(in) :: beta_3dvar ! Hybrid weight ! ! External subroutines ! ! (PDAF-internal names, real names are defined in the call to PDAF) EXTERNAL :: U_prodRinvA, & ! Provide product R^-1 A U_cvt, & ! Apply control vector transform matrix to control vector (parameterized) U_cvt_adj, & ! Apply adjoint control vector transform matrix (parameterized) U_cvt_ens, & ! Apply control vector transform matrix to control vector (ensemble) U_cvt_adj_ens, & ! Apply adjoint control vector transform matrix (ensemble) U_obs_op_lin, & ! Linearized observation operator U_obs_op_adj ! Adjoint observation operator ! !CALLING SEQUENCE: ! Called by: PDAF_3dvar_analysis_cvt ! Calls: PDAF_timeit ! Calls: PDAF_memcount !EOP ! *** local variables *** INTEGER :: iter ! Counter INTEGER, SAVE :: allocflag = 0 ! Flag whether first time allocation is done INTEGER :: dim_cv_p ! Full size of control vector REAL :: J_tot ! Cost function REAL, ALLOCATABLE :: gradJ_p(:) ! PE-local part of gradient of J REAL, ALLOCATABLE :: v_p(:) ! PE-local full control vector ! Variables for LFBGS INTEGER :: m = 5 ! Number of corrections used in limited memory matrix; 3<=m<=20 recommended INTEGER :: iprint CHARACTER(len=60) :: task, csave LOGICAL :: lsave(4) INTEGER :: isave(44) REAL :: factr = 1.0e+7 ! Tolerance in termination test REAL :: pgtol = 1.0e-5 ! Limit for stopping iterations REAL :: dsave(29) INTEGER, ALLOCATABLE :: nbd(:), iwa(:) REAL, ALLOCATABLE :: lvec(:), uvec(:), wa(:) ! ********************** ! *** INITIALIZATION *** ! ********************** IF (debug>0) & WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_hyb3dvar_optim_LBFGS -- START' ! Initialize overall dimension of control vector dim_cv_p = dim_cv_par_p + dim_cv_ens_p ! Set verbosity of solver IF (screen>0 .AND. screen<2) THEN iprint = -1 ELSEIF (screen<3) THEN iprint = 0 IF (mype>0) iprint = -1 ELSE iprint = 99 END IF IF (debug>0) iprint = 99 ! Allocate arrays ALLOCATE(v_p(dim_cv_p)) ALLOCATE(nbd(dim_cv_p), lvec(dim_cv_p), uvec(dim_cv_p)) ALLOCATE (iwa(3*dim_cv_p)) ALLOCATE (wa(2*m*dim_cv_p + 5*dim_cv_p + 11*m*m + 8*m)) IF (allocflag == 0) CALL PDAF_memcount(3, 'r', 12*dim_cv_p + 2*m*dim_cv_p + 11*m*m + 8*m) ! Settings for LBGFS m = m_lbfgs_var factr = factr_lbfgs_var pgtol = pgtol_lbfgs_var nbd = 0 ! Values are unbounded task = 'START' iter = 1 IF (debug>0) THEN WRITE (*,*) '++ PDAF-debug PDAF_hyb3dvar_optim_LBFGS', debug, & 'Solver config: m ', m WRITE (*,*) '++ PDAF-debug PDAF_hyb3dvar_optim_LBFGS', debug, & 'Solver config: factr', factr WRITE (*,*) '++ PDAF-debug PDAF_hyb3dvar_optim_LBFGS', debug, & 'Solver config: pgtol', pgtol END IF ! Initialize control vector v_p = 0.0 ! *************************** ! *** Iterative solving *** ! *************************** ! Prepare arrays for iterations ALLOCATE(gradJ_p(dim_cv_p)) IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_cv_p) IF (mype==0 .AND. screen > 0) & WRITE (*, '(a, 5x, a)') 'PDAF', '--- OPTIMIZE' minloop: DO IF (.NOT.(task(1:2).EQ.'FG'.OR.task.EQ.'NEW_X'.OR. & task.EQ.'START') ) THEN IF (mype==0 .AND. screen > 0) & WRITE (*,'(a, 5x, a, a)') 'PDAF', '--- Exit optimization, status ', task EXIT minloop END IF ! LBFGS CALL PDAF_timeit(54, 'new') CALL setulb(dim_cv_p, m, v_p, lvec, uvec, nbd, & J_tot, gradJ_p, factr, pgtol, & wa, iwa, task, iprint,& csave, lsave, isave, dsave ) CALL PDAF_timeit(54, 'old') ! ******************************** ! *** Evaluate cost function *** ! ******************************** CALL PDAF_timeit(53, 'new') CALL PDAF_hyb3dvar_costf_cvt(step, iter, dim_p, dim_ens, & dim_cv_p, dim_cv_par_p, dim_cv_ens_p, dim_obs_p, ens_p, obs_p, & dy_p, v_par_p, v_ens_p, v_p, J_tot, gradJ_p, & U_prodRinvA, U_cvt, U_cvt_adj, U_cvt_ens, U_cvt_adj_ens, & U_obs_op_lin, U_obs_op_adj, opt_parallel, beta_3dvar) CALL PDAF_timeit(53, 'old') IF (mype==0 .AND. screen >2) & WRITE (*,'(a, 8x, a, i5, es12.4)') 'PDAF', '--- iter, J: ', iter, J_tot iter = iter + 1 END DO minloop ! ******************** ! *** Finishing up *** ! ******************** DEALLOCATE(gradJ_p) DEALLOCATE(nbd, lvec, uvec, iwa, wa) IF (allocflag == 0) allocflag = 1 IF (debug>0) & WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_hyb3dvar_optim_LBFGS -- END' END SUBROUTINE PDAF_hyb3dvar_optim_lbfgs