From aa263e6a2adc9414142551a01bbcc46abb6a3894 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 15 Mar 2024 18:43:54 +0100
Subject: [PATCH] Perfect foresight: allow stack_solve_algo={2,3} without block
nor bytecode
---
doc/manual/source/the-model-file.rst | 6 ++---
.../perfect_foresight_solver_core.m | 7 ++++--
.../private/check_input_arguments.m | 8 +++----
matlab/perfect-foresight-models/sim1.m | 23 +++++++++++++++----
tests/run_block_bytecode_tests.m | 4 ++--
5 files changed, 30 insertions(+), 18 deletions(-)
diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index a47791f929..aa9d5fd3aa 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -3706,15 +3706,13 @@ speed-up on large models.
Use a Newton algorithm with a Generalized Minimal Residual
(GMRES) solver at each iteration, applied on the stacked system
- of all equations in all periods (requires ``bytecode`` and/or
- ``block`` option, see :ref:`model-decl`)
+ of all equations in all periods.
``3``
Use a Newton algorithm with a Stabilized Bi-Conjugate Gradient
(BiCGStab) solver at each iteration, applied on the stacked
- system of all equations in all periods (requires ``bytecode``
- and/or ``block`` option, see :ref:`model-decl`).
+ system of all equations in all periods.
``4``
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 8d6ec41b2f..2c7202f6a1 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -17,7 +17,7 @@ function [y, success, maxerror, iter, per_block_status] = perfect_foresight_solv
% - iter [integer] Number of iterations of the underlying nonlinear solver (empty for non-iterative methods)
% - per_block_status [struct] In the case of block decomposition, provides per-block solver status information (empty if no block decomposition)
-% Copyright © 2015-2023 Dynare Team
+% Copyright © 2015-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -96,8 +96,11 @@ else
[y, success] = sim1_purely_static(y, exo_simul, steady_state, M_, options_);
else % General case
switch options_.stack_solve_algo
- case 0
+ case {0 2 3}
if options_.linear_approximation
+ if ismember(options_.stack_solve_algo, [2 3])
+ error('Invalid value of stack_solve_algo option!')
+ end
[y, success, maxerror] = sim1_linear(y, exo_simul, steady_state, exo_steady_state, M_, options_);
else
[y, success, maxerror, iter] = sim1(y, exo_simul, steady_state, M_, options_);
diff --git a/matlab/perfect-foresight-models/private/check_input_arguments.m b/matlab/perfect-foresight-models/private/check_input_arguments.m
index ccf4a34a56..b77de81f46 100644
--- a/matlab/perfect-foresight-models/private/check_input_arguments.m
+++ b/matlab/perfect-foresight-models/private/check_input_arguments.m
@@ -7,7 +7,7 @@ function check_input_arguments(options_, M_, oo_)
% M_ [structure] describing the model
% oo_ [structure] storing the results
-% Copyright © 2015-2023 Dynare Team
+% Copyright © 2015-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -28,10 +28,8 @@ if options_.stack_solve_algo < 0 || options_.stack_solve_algo > 7
error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: stack_solve_algo must be between 0 and 7')
end
-if ~options_.block && ~options_.bytecode && options_.stack_solve_algo ~= 0 ...
- && options_.stack_solve_algo ~= 1 && options_.stack_solve_algo ~= 6 ...
- && options_.stack_solve_algo ~= 7
- error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: you must use stack_solve_algo={0,1,6,7} when not using block nor bytecode option')
+if ~options_.block && ~options_.bytecode && ~ismember(options_.stack_solve_algo, [0:3 6 7])
+ error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: you must use stack_solve_algo={0,1,2,3,6,7} when not using block nor bytecode option')
end
if options_.block && ~options_.bytecode && options_.stack_solve_algo == 5
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index ad7eb61dc3..6fae2f151e 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -16,7 +16,7 @@ function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, e
% - err [double] ∞-norm of the residual
% - iter [integer] Number of iterations
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -128,7 +128,7 @@ for iter = 1:options_.simul.maxit
if options_.simul.robust_lin_solve
dy = -lin_solve_robust(A, res, verbose, options_);
else
- dy = -lin_solve(A, res, verbose);
+ dy = -lin_solve(A, res, verbose, options_);
end
if any(isnan(dy)) || any(isinf(dy))
if verbose
@@ -180,13 +180,26 @@ if verbose
skipline();
end
-function x = lin_solve(A, b, verbose)
+function x = lin_solve(A, b, verbose, options_)
if norm(b) < sqrt(eps) % then x = 0 is a solution
x = 0;
return
end
-x = A\b;
+if options_.stack_solve_algo == 0
+ x = A\b;
+else
+ ilu_setup.type = 'ilutp';
+ ilu_setup.droptol = 1e-10;
+ [L1, U1] = ilu(A, ilu_setup);
+ if options_.stack_solve_algo == 2
+ x = gmres(A, b, [], [], [], L1, U1);
+ elseif options_.stack_solve_algo == 3
+ x = bicgstab(A, b, [], [], L1, U1);
+ else
+ error('sim1: invalid value for options_.stack_solve_algo')
+ end
+end
x(~isfinite(x)) = 0;
relres = norm(b - A*x) / norm(b);
if relres > 1e-6 && verbose
@@ -343,4 +356,4 @@ if rank_jacob < size(jacob,1)
fprintf('Equation %5u, period %5u\n',equation(ii),period(ii))
end
end
-end
\ No newline at end of file
+end
diff --git a/tests/run_block_bytecode_tests.m b/tests/run_block_bytecode_tests.m
index 820db3f843..e9aac91947 100644
--- a/tests/run_block_bytecode_tests.m
+++ b/tests/run_block_bytecode_tests.m
@@ -1,4 +1,4 @@
-% Copyright © 2011-2023 Dynare Team
+% Copyright © 2011-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -42,7 +42,7 @@ for blockFlag = 0:1
default_stack_solve_algo = 0;
if ~blockFlag && storageFlag ~= 2
solve_algos = [1:4 9];
- stack_solve_algos = [0 1 6];
+ stack_solve_algos = [0:3 6];
elseif blockFlag && storageFlag ~= 2
solve_algos = [1:4 6:9];
stack_solve_algos = [0:4 6];
--
GitLab