diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index c3d921b924541fe96f2e4ec42b6f672148b397a1..23813d1468a689f2a0d4d44d09f7098d31d0a480 100644
--- a/matlab/default_option_values.m
+++ b/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;
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 8f44513b3d7954bba306e6316252d5b46f3f617c..85d15d507f1329cbefbd5e12679835ac807a10c6 100644
--- a/matlab/dynare_solve.m
+++ b/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]')
diff --git a/mex/sources/block_trust_region/mexFunction.f08 b/mex/sources/block_trust_region/mexFunction.f08
index 507f477146df8b451348de2bc7b2e328cbb96aa3..f36238a372fe613160e6ab58dd970f0d06a06580 100644
--- a/mex/sources/block_trust_region/mexFunction.f08
+++ b/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
diff --git a/mex/sources/block_trust_region/trust_region.f08 b/mex/sources/block_trust_region/trust_region.f08
index 8ad879a170d26c88a39be0dd263c69414c70199a..dfb450cbc7a50de7273fa4d2ae9ef520d7db3eb0 100644
--- a/mex/sources/block_trust_region/trust_region.f08
+++ b/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