Add option for trust region algorithm (mex).

trust_region_initial_step_bound_factor is determinining thye initial
step bound. Default value is 100 (previous hard coded value was 1).

matlab/default_option_values.m

@@ -53,6 +53,7 @@ options_.qz_zero_threshold = 1e-6;
 options_.lyapunov_complex_threshold = 1e-15;
 options_.solve_tolf = eps^(1/3);
 options_.solve_tolx = eps^(2/3);
+options_.trust_region_initial_step_bound_factor = 100;
 options_.dr_display_tol=1e-6;
 options_.minimal_workspace = false;
 options_.dp.maxit = 3000;

matlab/dynare_solve.m

@@ -370,7 +370,7 @@ elseif ismember(options.solve_algo, [13, 14])
         auxstruct.isloggedlhs = isloggedlhs;
         auxstruct.isauxdiffloggedrhs = isauxdiffloggedrhs;
     end
-    [x, errorflag] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.debug, auxstruct, arguments{:});
+    [x, errorflag] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, auxstruct, arguments{:});
     [fvec, fjac] = feval(f, x, arguments{:});
 else
     error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')

mex/sources/block_trust_region/mexFunction.f08

@@ -31,7 +31,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   real(real64), dimension(:), allocatable, target :: x
   type(dm_block), dimension(:), allocatable, target :: blocks
   integer :: info, i
-  real(real64) :: tolf, tolx
+  real(real64) :: tolf, tolx, factor
   integer :: maxiter
   real(real64), dimension(:), allocatable :: fvec
   real(real64), dimension(:,:), allocatable :: fjac
@@ -67,22 +67,28 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
      call mexErrMsgTxt("Fifth argument (maxiter) should be a numeric scalar")
   end if
 
-  if (.not. (mxIsLogicalScalar(prhs(6)))) then
-     call mexErrMsgTxt("Sixth argument (debug) should be a logical scalar")
+  if (.not. (mxIsScalar(prhs(6)) .and. mxIsNumeric(prhs(6)))) then
+     call mexErrMsgTxt("Sixth argument (factor) should be a numeric scalar")
   end if
 
-  if (.not. (mxIsStruct(prhs(7)) .and. &
-       (mxGetNumberOfFields(prhs(7)) == 0 .or. mxGetNumberOfFields(prhs(7)) == 4))) then
-     call mexErrMsgTxt("Seventh argument should be a struct with either 0 or 4 fields")
+  if (.not. (mxIsLogicalScalar(prhs(7)))) then
+     call mexErrMsgTxt("Seventh argument (debug) should be a logical scalar")
   end if
-  specializedunivariateblocks = (mxGetNumberOfFields(prhs(7)) == 4)
+
+  if (.not. (mxIsStruct(prhs(8)) .and. &
+       (mxGetNumberOfFields(prhs(8)) == 0 .or. mxGetNumberOfFields(prhs(8)) == 4))) then
+     call mexErrMsgTxt("Eighth argument should be a struct with either 0 or 4 fields")
+  end if
+
+  specializedunivariateblocks = (mxGetNumberOfFields(prhs(8)) == 4)
 
   func => prhs(1)
   tolf = mxGetScalar(prhs(3))
   tolx = mxGetScalar(prhs(4))
   maxiter = int(mxGetScalar(prhs(5)))
-  debug = mxGetScalar(prhs(6)) == 1._c_double
-  extra_args => prhs(8:nrhs) ! Extra arguments to func are in argument 8 and subsequent ones
+  factor = mxGetScalar(prhs(6))
+  debug = mxGetScalar(prhs(7)) == 1._c_double
+  extra_args => prhs(9:nrhs) ! Extra arguments to func are in argument 9 and subsequent ones
   associate (x_mat => mxGetPr(prhs(2)))
     allocate(x(size(x_mat)))
     x = x_mat
@@ -90,25 +96,25 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   if (specializedunivariateblocks) then
      block
        type(c_ptr) :: tmp
-       tmp = mxGetField(prhs(7), 1_mwIndex, "isloggedlhs")
+       tmp = mxGetField(prhs(8), 1_mwIndex, "isloggedlhs")
        if (.not. (c_associated(tmp) .and. mxIsLogical(tmp) .and. mxGetNumberOfElements(tmp) == size(x))) then
           call mexErrMsgTxt("Seventh argument must have a 'isloggedlhs' field of type logical, of same size as second argument")
        end if
        isloggedlhs => mxGetLogicals(tmp)
 
-       tmp = mxGetField(prhs(7), 1_mwIndex, "isauxdiffloggedrhs")
+       tmp = mxGetField(prhs(8), 1_mwIndex, "isauxdiffloggedrhs")
        if (.not. (c_associated(tmp) .and. mxIsLogical(tmp) .and. mxGetNumberOfElements(tmp) == size(x))) then
           call mexErrMsgTxt("Seventh argument must have a 'isauxdiffloggedrhs' field of type &
                &logical, of same size as second argument")
        end if
        isauxdiffloggedrhs => mxGetLogicals(tmp)
 
-       lhs = mxGetField(prhs(7), 1_mwIndex, "lhs")
+       lhs = mxGetField(prhs(8), 1_mwIndex, "lhs")
        if (.not. (c_associated(lhs) .and. mxIsCell(lhs) .and. mxGetNumberOfElements(lhs) == size(x))) then
           call mexErrMsgTxt("Seventh argument must have a 'lhs' field of type cell, of same size as second argument")
        end if
 
-       endo_names = mxGetField(prhs(7), 1_mwIndex, "endo_names")
+       endo_names = mxGetField(prhs(8), 1_mwIndex, "endo_names")
        if (.not. (c_associated(endo_names) .and. mxIsCell(endo_names) .and. mxGetNumberOfElements(endo_names) == size(x))) then
           call mexErrMsgTxt("Seventh argument must have a 'endo_names' field of type cell, of same size as second argument")
        end if
@@ -203,7 +209,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
           call mexErrMsgTxt("Non-square block")
        end if
        x_block = x(x_indices)
-       call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter)
+       call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter, factor)
        x(x_indices) = x_block
      end block
 
@@ -222,7 +228,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   if (maxval(abs(fvec)) > tolf) then
      if (debug) &
           call mexPrintf_trim_newline("DYNARE_SOLVE (solve_algo=13|14): residuals still too large, solving for the whole model")
-     call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter)
+     call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter, factor)
   else
      info = 1
   end if

mex/sources/block_trust_region/trust_region.f08

@@ -28,7 +28,7 @@ module trust_region
   public :: trust_region_solve
 
 contains
-  subroutine trust_region_solve(x, f, info, tolx, tolf, maxiter)
+  subroutine trust_region_solve(x, f, info, tolx, tolf, maxiter, factor)
 
     real(real64), dimension(:), intent(inout) :: x ! On entry: guess value; on exit: solution
     interface
@@ -47,9 +47,11 @@ contains
     integer, intent(out) :: info
     real(real64), intent(in), optional :: tolx, tolf ! Tolerances in x and f
     integer, intent(in), optional :: maxiter ! Maximum number of iterations
+    real(real64), intent(in), optional :: factor ! Used in determining the initial step bound.
 
     real(real64) :: tolx_actual, tolf_actual
     integer :: maxiter_actual
+    real(real64) :: factor_actual
     integer, parameter :: maxslowiter = 15 ! Maximum number of consecutive iterations with slow progress
     real(real64), parameter :: macheps = epsilon(x)
     real(real64) :: delta ! Radius of the trust region
@@ -79,6 +81,11 @@ contains
     else
        maxiter_actual = 50
     end if
+    if (present(factor)) then
+       factor_actual = factor
+    else
+       factor_actual = 100.0_real64
+    end if
 
     niter = 1
     ncsucc = 0
@@ -123,6 +130,7 @@ contains
             else
                delta = 1
             end if
+            delta = delta*factor
          end if
 
          ! Get trust-region model (dogleg) minimizer