Perfect foresight: allow stack_solve_algo={2,3} without block nor bytecode

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``

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_);

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

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

tests/run_block_bytecode_tests.m

-% 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];