From 801a9d76616418a62d10869e1b85fd5651f4d365 Mon Sep 17 00:00:00 2001
From: Normann Rion <normann@dynare.org>
Date: Wed, 29 May 2024 12:07:12 +0200
Subject: [PATCH] Provide `dr_cycle_reduction_maxiter` option See #1920
---
doc/manual/source/the-model-file.rst | 15 +++++++++++++++
matlab/+mom/default_option_mom_values.m | 1 +
matlab/default_option_values.m | 1 +
.../missing/mex/cycle_reduction/cycle_reduction.m | 15 +++++++--------
matlab/stochastic_solver/dyn_first_order_solver.m | 2 +-
mex/sources/cycle_reduction/mexFunction.f08 | 11 ++++++-----
tests/cyclereduction.m | 4 ++--
7 files changed, 33 insertions(+), 16 deletions(-)
diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index a85bc37ca8..e85e77d745 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -4837,6 +4837,11 @@ Computing the stochastic solution
The convergence criterion used in the cycle reduction
algorithm. Its default value is ``1e-7``.
+ .. option:: dr_cycle_reduction_maxiter = INTEGER
+
+ The maximum number of iterations used in the cycle
+ reduction algorithm. Its default value is ``100``.
+
.. option:: dr_logarithmic_reduction_tol = DOUBLE
The convergence criterion used in the logarithmic reduction
@@ -7998,6 +8003,11 @@ observed variables.
See :opt:`dr_cycle_reduction_tol <dr_cycle_reduction_tol = DOUBLE>`.
Default: ``1e-7``.
+ .. option:: dr_cycle_reduction_maxiter = INTEGER
+
+ See :opt:`dr_cycle_reduction_maxiter <dr_cycle_reduction_maxiter = INTEGER>`.
+ Default: ``100``.
+
.. option:: dr_logarithmic_reduction_tol = DOUBLE
See :opt:`dr_logarithmic_reduction_tol <dr_logarithmic_reduction_tol = DOUBLE>`.
@@ -9926,6 +9936,11 @@ Numerical algorithms options
See :opt:`dr_cycle_reduction_tol <dr_cycle_reduction_tol = DOUBLE>`.
Default: ``1e-7``.
+ .. option:: dr_cycle_reduction_maxiter = INTEGER
+
+ See :opt:`dr_cycle_reduction_maxiter <dr_cycle_reduction_maxiter = INTEGER>`.
+ Default: ``100``.
+
.. option:: dr_logarithmic_reduction_tol = DOUBLE
See :opt:`dr_logarithmic_reduction_tol <dr_logarithmic_reduction_tol = DOUBLE>`.
diff --git a/matlab/+mom/default_option_mom_values.m b/matlab/+mom/default_option_mom_values.m
index da9affb5e7..e5f056a04c 100644
--- a/matlab/+mom/default_option_mom_values.m
+++ b/matlab/+mom/default_option_mom_values.m
@@ -175,6 +175,7 @@ options_mom_ = set_default_option(options_mom_,'pruning',false);
options_mom_ = set_default_option(options_mom_,'aim_solver',false); % use AIM algorithm to compute perturbation approximation instead of mjdgges
options_mom_ = set_default_option(options_mom_,'k_order_solver',false); % use k_order_perturbation instead of mjdgges
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction',false); % use cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
+options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_maxiter',100); % maximum number of iterations used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_tol',1e-7); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction',false); % use logarithmic reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_maxiter',100); % maximum number of iterations used in the logarithmic reduction algorithm
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 13c4275fda..5af50faa27 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -689,6 +689,7 @@ options_.dr_logarithmic_reduction_tol = 1e-12;
% convergence criterion for iteratives methods to solve the decision rule
options_.dr_logarithmic_reduction_maxiter = 100;
+options_.dr_cycle_reduction_maxiter = 100;
% dates for historical time series
options_.initial_date = dates();
diff --git a/matlab/missing/mex/cycle_reduction/cycle_reduction.m b/matlab/missing/mex/cycle_reduction/cycle_reduction.m
index d3e126de2b..8f740aa38d 100644
--- a/matlab/missing/mex/cycle_reduction/cycle_reduction.m
+++ b/matlab/missing/mex/cycle_reduction/cycle_reduction.m
@@ -1,4 +1,4 @@
-function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch)
+function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, max_it, ch)
%@info:
%! @deftypefn {Function File} {[@var{X}, @var{info}] =} cycle_reduction (@var{A0},@var{A1},@var{A2},@var{cvg_tol},@var{ch})
@@ -60,14 +60,13 @@ function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-max_it = 300;
it = 0;
info = 0;
X = [];
crit = Inf;
A0_0 = A0;
Ahat1 = A1;
-if (nargin == 5 && ~isempty(ch) )
+if (nargin == 6 && ~isempty(ch) )
A1_0 = A1;
A2_0 = A2;
end
@@ -94,7 +93,7 @@ while cont
return
elseif isnan(crit)
info(1) = 402;
- info(2) = log(norm(A1,1))
+ info(2) = log(norm(A1,1));
return
end
it = it + 1;
@@ -102,12 +101,12 @@ end
X = -Ahat1\A0_0;
-if (nargin == 5 && ~isempty(ch) )
+if (nargin == 6 && ~isempty(ch) )
%check the solution
res = norm(A0_0 + A1_0 * X + A2_0 * X * X, 1);
if (res > cvg_tol)
- info(1) = 403
- info(2) = log(res)
+ info(1) = 403;
+ info(2) = log(res);
dprintf('The norm of the residual is %s whereas the tolerance criterion is %s', num2str(res), num2str(cvg_tol));
end
end
@@ -128,7 +127,7 @@ C = diag(15*ones(n,1)); C = C - diag(5*ones(n-1,1),-1); C = C - diag(5*ones(n-1,
% Solve the equation with the cycle reduction algorithm
try
- tic; X1 = cycle_reduction(C,B,A,1e-7); elapsedtime = toc;
+ tic; X1 = cycle_reduction(C,B,A,1e-7,100); elapsedtime = toc;
disp(['Elapsed time for cycle reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').'])
t(1) = 1;
catch
diff --git a/matlab/stochastic_solver/dyn_first_order_solver.m b/matlab/stochastic_solver/dyn_first_order_solver.m
index 7d39307bd9..3adfe0e70d 100644
--- a/matlab/stochastic_solver/dyn_first_order_solver.m
+++ b/matlab/stochastic_solver/dyn_first_order_solver.m
@@ -158,7 +158,7 @@ if task ~= 1 && (options_.dr_cycle_reduction || options_.dr_logarithmic_reductio
B1 = [aa(row_indx,index_0m) aa(row_indx,index_0p) ];
C1 = [zeros(ndynamic,npred) aa(row_indx,index_p)];
if options_.dr_cycle_reduction
- [ghx, info] = cycle_reduction(A1, B1, C1, options_.dr_cycle_reduction_tol);
+ [ghx, info] = cycle_reduction(A1, B1, C1, options_.dr_cycle_reduction_tol, options_.dr_cycle_reduction_maxiter);
else
[ghx, info] = logarithmic_reduction(C1, B1, A1, options_.dr_logarithmic_reduction_tol, options_.dr_logarithmic_reduction_maxiter);
end
diff --git a/mex/sources/cycle_reduction/mexFunction.f08 b/mex/sources/cycle_reduction/mexFunction.f08
index 5e507ba195..8399f32a81 100644
--- a/mex/sources/cycle_reduction/mexFunction.f08
+++ b/mex/sources/cycle_reduction/mexFunction.f08
@@ -28,16 +28,17 @@ module c_reduction
contains
! Cycle reduction algorithm
- subroutine cycle_reduction(A0, A1, A2, X, cvg_tol, check, info)
+ subroutine cycle_reduction(A0, A1, A2, X, cvg_tol, check, max_it, info)
real(real64), dimension(:,:), intent(in) :: A0, A1, A2
real(real64), dimension(:,:), intent(inout) :: X
real(real64), intent(in) :: cvg_tol
logical, intent(in) :: check
+ integer, intent(in) :: max_it
real(c_double), dimension(2), intent(inout) :: info
real(real64), dimension(:,:), allocatable :: A02, &
Q0, Q2, Ahat, invA1_A02, A1i, A0_tmp, A1_tmp
- integer :: it, n, dn, max_it
+ integer :: it, n, dn
integer(blint) :: info_inv
integer(blint), dimension(:), allocatable :: ipiv
real(real64) :: residual, crit
@@ -51,7 +52,6 @@ contains
Ahat = A1
A1i = A1
it = 0
- max_it = 300
loop: do
! Computing [A0;A2]*(A1\[A0 A2])
A1_tmp = A1i
@@ -118,7 +118,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
type(c_ptr), dimension(*), intent(out) :: plhs
integer(c_int), intent(in), value :: nlhs, nrhs
- integer :: i, n
+ integer :: i, n, max_it
character(kind=c_char, len=2) :: num2str
real(real64) :: cvg_tol
real(real64), dimension(:,:), pointer, contiguous :: A0, A1, A2, X
@@ -174,13 +174,14 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
A1(1:n,1:n) => mxGetPr(prhs(2))
A2(1:n,1:n) => mxGetPr(prhs(3))
cvg_tol = mxGetScalar(prhs(4))
+ max_it = int(mxGetScalar(prhs(5)))
info = [0._c_double,0._c_double]
plhs(1) = mxCreateDoubleMatrix(int(n, mwSize), int(n, mwSize), mxREAL)
X(1:n,1:n) => mxGetPr(plhs(1))
! 2. Calling the Cycle Reduction algorithm
- call cycle_reduction(A0, A1, A2, X, cvg_tol, check, info)
+ call cycle_reduction(A0, A1, A2, X, cvg_tol, check, max_it, info)
! 3. Editing the information output if necessary
if (nlhs == 2) then
diff --git a/tests/cyclereduction.m b/tests/cyclereduction.m
index 5c93a18fbc..a3f8ea3739 100644
--- a/tests/cyclereduction.m
+++ b/tests/cyclereduction.m
@@ -28,7 +28,7 @@ X2 = zeros(n,n);
% 1. Solve the equation with the Matlab cycle reduction algorithm
tElapsed1 = 0.;
try
- tic; [X1,info] = cycle_reduction_matlab(C,B,A,cvg_tol,[0.]); tElapsed1 = toc;
+ tic; [X1,info] = cycle_reduction_matlab(C,B,A,cvg_tol,100); tElapsed1 = toc;
disp(['Elapsed time for the Matlab cycle reduction algorithm is: ' num2str(tElapsed1) ' (n=' int2str(n) ').'])
R = norm(C+B*X1+A*X1*X1,1);
if (R > cvg_tol)
@@ -43,7 +43,7 @@ end
% 2. Solve the equation with the Fortran cycle reduction algorithm
tElapsed2 = 0.;
try
- tic; [X2,info] = cycle_reduction(C,B,A,cvg_tol,[0.]); tElapsed2 = toc;
+ tic; [X2,info] = cycle_reduction(C,B,A,cvg_tol,100); tElapsed2 = toc;
disp(['Elapsed time for the Fortran cycle reduction algorithm is: ' num2str(tElapsed2) ' (n=' int2str(n) ').'])
R = norm(C+B*X2+A*X2*X2,1);
if (R > cvg_tol)
--
GitLab