Compare revisions

9e8427e7 · 9e8427e7 · 9e8427e7 · 9e8427e7 · 9e8427e7 · 9e8427e7
--- a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
@@ -17,17 +17,19 @@ function [endogenousvariables, success] = solve_stacked_linear_problem(endogenou
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+periods = get_simulation_periods(options_);
 if M_.maximum_lag > 0
    y0 = endogenousvariables(:, M_.maximum_lag);
 else
    y0 = NaN(M_.endo_nbr, 1);
 end
 if M_.maximum_lead > 0
-    yT = endogenousvariables(:, M_.maximum_lag+options_.periods+1);
+    yT = endogenousvariables(:, M_.maximum_lag+periods+1);
 else
    yT = NaN(M_.endo_nbr, 1);
 end
-z = endogenousvariables(:,M_.maximum_lag+(1:options_.periods));
+z = endogenousvariables(:,M_.maximum_lag+(1:periods));
 % Evaluate the residuals and Jacobian of the dynamic model at the deterministic steady state.
 y3n = repmat(steadystate_y, 3, 1);
@@ -50,7 +52,7 @@ x = bsxfun(@minus, exogenousvariables, steadystate_x');
                                           options_, ...
                                           jacobian, y0-steadystate_y, yT-steadystate_y, ...
                                           x, M_.params, steadystate_y, ...
-                                           M_.maximum_lag, options_.periods, M_.endo_nbr);
+                                           M_.maximum_lag, periods, M_.endo_nbr);
 if all(imag(y)<.1*options_.dynatol.x)
    if ~isreal(y)
@@ -60,7 +62,7 @@ else
    check = 1;
 end
-endogenousvariables = [y0 bsxfun(@plus,reshape(y,M_.endo_nbr,options_.periods), steadystate_y) yT];
+endogenousvariables = [y0 bsxfun(@plus,reshape(y,M_.endo_nbr,periods), steadystate_y) yT];
 success = ~check;

--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
-function [endogenousvariables, success, maxerror] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M_, options_)
+function [endogenousvariables, success, maxerror, exogenousvariables] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, controlled_paths_by_period, M_, options_)
 % [endogenousvariables, success, maxerror] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M_, options_)
 % Solves the perfect foresight model using dynare_solve
 %
 % INPUTS
-% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial guess).
+% - endogenousvariables [double] N*(T+M_.maximum_lag+M_.maximum_lead) array, paths for the endogenous variables (initial guess).
-% - exogenousvariables  [double] T*M array, paths for the exogenous variables.
+% - exogenousvariables  [double] (T+M_.maximum_lag+M_.maximum_lead)*M array, paths for the exogenous variables.
 % - steadystate         [double] N*1 array, steady state for the endogenous variables.
+%   - controlled_paths_by_period [struct] data from perfect_foresight_controlled_paths block
 % - M_                   [struct] contains a description of the model.
 % - options_             [struct] contains various options.
 %
 % OUTPUTS
-% - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
+% - endogenousvariables [double] N*(T+M_.maximum_lag+M_.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model).
 % - success             [logical] Whether a solution was found
 % - maxerror            [double] 1-norm of the residual
+% - exogenousvariables  [double] (T+M_.maximum_lag+M_.maximum_lead)*M array, paths for the exogenous variables
+%                                (may be modified if perfect_foresight_controlled_paths present)
-% Copyright © 2015-2024 Dynare Team
+% Copyright © 2015-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -31,30 +34,33 @@ function [endogenousvariables, success, maxerror] = solve_stacked_problem(endoge
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+periods = get_simulation_periods(options_);
 if M_.maximum_lag > 0
    y0 = endogenousvariables(:, M_.maximum_lag);
 else
    y0 = NaN(M_.endo_nbr, 1);
 end
 if M_.maximum_lead > 0
-    yT = endogenousvariables(:, M_.maximum_lag+options_.periods+1);
+    yT = endogenousvariables(:, M_.maximum_lag+periods+1);
 else
    yT = NaN(M_.endo_nbr, 1);
 end
-z = endogenousvariables(:,M_.maximum_lag+(1:options_.periods));
+z = endogenousvariables(:,M_.maximum_lag+(1:periods));
 if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementarity problem
+    assert(isempty(controlled_paths_by_period))
    [lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), M_.params);
    if options_.linear_approximation
        lb = lb - steadystate_y;
        ub = ub - steadystate_y;
    end
    if options_.solve_algo == 10
-        options_.lmmcp.lb = repmat(lb,options_.periods,1);
+        options_.lmmcp.lb = repmat(lb,periods,1);
-        options_.lmmcp.ub = repmat(ub,options_.periods,1);
+        options_.lmmcp.ub = repmat(ub,periods,1);
    elseif options_.solve_algo == 11
-        options_.mcppath.lb = repmat(lb,options_.periods,1);
+        options_.mcppath.lb = repmat(lb,periods,1);
-        options_.mcppath.ub = repmat(ub,options_.periods,1);
+        options_.mcppath.ub = repmat(ub,periods,1);
    end
    dynamic_resid_function = str2func([M_.fname,'.sparse.dynamic_resid']);
    dynamic_g1_function = str2func([M_.fname,'.sparse.dynamic_g1']);
@@ -64,14 +70,33 @@ if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementari
                                                 dynamic_resid_function, dynamic_g1_function, y0, yT, ...
                                                 exogenousvariables, M_.params, steadystate, ...
                                                 M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, M_.dynamic_g1_sparse_colptr, ...
-                                                 M_.maximum_lag, options_.periods, M_.endo_nbr, ...
+                                                 M_.maximum_lag, periods, M_.endo_nbr, ...
                                                 M_.dynamic_mcp_equations_reordering);
    eq_to_ignore=find(isfinite(lb) | isfinite(ub));
-else
+elseif isempty(controlled_paths_by_period)
    [y, check, res, ~, errorcode] = dynare_solve(@perfect_foresight_problem, z(:), ...
                                               options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
-                                               options_, y0, yT, exogenousvariables, M_.params, steadystate, options_.periods, M_, options_);
+                                               options_, y0, yT, exogenousvariables, M_.params, steadystate, periods, M_, options_);
+else % controlled paths
+    for p = 1:periods
+        if isempty(controlled_paths_by_period(p).exogenize_id)
+            continue
+        end
+        z(controlled_paths_by_period(p).exogenize_id,p) = exogenousvariables(p+M_.maximum_lag,controlled_paths_by_period(p).endogenize_id);
+    end
+    [y, check, res, ~, errorcode] = dynare_solve(@controlled_paths_wrapper, z(:), ...
+                                                 options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
+                                                 options_, y0, yT, exogenousvariables, M_.params, steadystate, periods, controlled_paths_by_period, M_, options_);
+    for p = 1:periods
+        if isempty(controlled_paths_by_period(p).exogenize_id)
+            continue
+        end
+        exogenousvariables(p+M_.maximum_lag,controlled_paths_by_period(p).endogenize_id) = y(controlled_paths_by_period(p).exogenize_id+(p-1)*M_.endo_nbr);
+        y(controlled_paths_by_period(p).exogenize_id+(p-1)*M_.endo_nbr) = controlled_paths_by_period(p).values;
+    end
 end
 if all(imag(y)<.1*options_.dynatol.x)
@@ -82,9 +107,9 @@ else
    check = 1;
 end
-endogenousvariables(:, M_.maximum_lag+(1:options_.periods)) = reshape(y, M_.endo_nbr, options_.periods);
+endogenousvariables(:, M_.maximum_lag+(1:periods)) = reshape(y, M_.endo_nbr, periods);
 residuals=zeros(size(endogenousvariables));
-residuals(:, M_.maximum_lag+(1:options_.periods)) = reshape(res, M_.endo_nbr, options_.periods);
+residuals(:, M_.maximum_lag+(1:periods)) = reshape(res, M_.endo_nbr, periods);
 if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementarity problem
    residuals(eq_to_ignore,bsxfun(@le, endogenousvariables(eq_to_ignore,:), lb(eq_to_ignore)+eps) | bsxfun(@ge,endogenousvariables(eq_to_ignore,:),ub(eq_to_ignore)-eps))=0;
 end
@@ -94,3 +119,24 @@ success = ~check;
 if ~success && options_.debug
    dprintf('solve_stacked_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
 end
+end
+function [res, A] = controlled_paths_wrapper(z, y0, yT, exogenousvariables, params, steadystate, periods, controlled_paths_by_period, M_, options_)
+y = z;
+for p = 1:periods
+    if isempty(controlled_paths_by_period(p).exogenize_id)
+        continue
+    end
+    exogenousvariables(p+M_.maximum_lag,controlled_paths_by_period(p).endogenize_id) = z(controlled_paths_by_period(p).exogenize_id+(p-1)*M_.endo_nbr);
+    y(controlled_paths_by_period(p).exogenize_id+(p-1)*M_.endo_nbr) = controlled_paths_by_period(p).values;
+end
+[res, A] = perfect_foresight_problem(y, y0, yT, exogenousvariables, params, steadystate, periods, M_, options_);
+A = controlled_paths_substitute_stacked_jacobian(A, y, y0, yT, exogenousvariables, steadystate, controlled_paths_by_period, M_);
+end
--- a/matlab/perfect-foresight-models/solve_two_boundaries_lbj.m
+++ b/matlab/perfect-foresight-models/solve_two_boundaries_lbj.m
@@ -25,7 +25,7 @@ function [y, T, success, err, iter] = solve_two_boundaries_lbj(fh, y, x, steady_
 %   simulation of dynamic models with forward variables through the use
 %   of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602.
-% Copyright © 2023 Dynare Team
+% Copyright © 2023-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -46,7 +46,7 @@ sparse_rowval = M_.block_structure.block(blk).g1_sparse_rowval;
 sparse_colval = M_.block_structure.block(blk).g1_sparse_colval;
 sparse_colptr = M_.block_structure.block(blk).g1_sparse_colptr;
-periods = options_.periods;
+periods = get_simulation_periods(options_);
 % NB: notations are deliberately similar to those of sim1_lbj.m

--- a/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
+++ b/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
@@ -25,7 +25,7 @@ function [y, T, success, max_res, iter] = solve_two_boundaries_stacked(fh, y, x,
 % ALGORITHM
 %   Newton with LU or GMRES or BiCGStab
-% Copyright © 1996-2024 Dynare Team
+% Copyright © 1996-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -44,7 +44,7 @@ function [y, T, success, max_res, iter] = solve_two_boundaries_stacked(fh, y, x,
 Blck_size = M_.block_structure.block(Block_Num).mfs;
 y_index = M_.block_structure.block(Block_Num).variable(end-Blck_size+1:end);
-periods = options_.periods;
+periods = get_simulation_periods(options_);
 y_kmin = M_.maximum_lag;
 stack_solve_algo = options_.stack_solve_algo;
@@ -57,11 +57,8 @@ verbose = options_.verbosity;
 cvg=false;
 iter=0;
 correcting_factor=0.01;
-ilu_setup.droptol=1e-10;
-ilu_setup.type = 'ilutp';
 max_resa=1e100;
 lambda = 1; % Length of Newton step (unused for stack_solve_algo=4)
-reduced = 0;
 while ~(cvg || iter > options_.simul.maxit)
    r = NaN(Blck_size, periods);
    g1a = spalloc(Blck_size*periods, Blck_size*periods, M_.block_structure.block(Block_Num).NNZDerivatives*periods);
@@ -81,7 +78,6 @@ while ~(cvg || iter > options_.simul.maxit)
            g1a((it_-y_kmin-1)*Blck_size+(1:Blck_size), (it_-y_kmin-2)*Blck_size+(1:3*Blck_size)) = g1;
        end
    end
-    preconditioner = 2;
    ya = reshape(y(y_index, y_kmin+(1:periods)), 1, periods*Blck_size)';
    ra = reshape(r, periods*Blck_size, 1);
    b=-ra+g1a*ya;
@@ -98,7 +94,7 @@ while ~(cvg || iter > options_.simul.maxit)
    end
    if ~cvg
        if iter>0
-            if ~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1)
+            if ~isreal(max_res) || isnan(max_res) %|| (max_resa<max_res && iter>1)
                if verbose && ~isreal(max_res)
                    disp(['Variable ' M_.endo_names{max_indx} ' (' int2str(max_indx) ') returns an undefined value']);
                end
@@ -128,7 +124,6 @@ while ~(cvg || iter > options_.simul.maxit)
                    end
                elseif lambda>1e-8 && stack_solve_algo ~= 4
                    lambda=lambda/2;
-                    reduced = 1;
                    if verbose
                        disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
                    end
@@ -155,102 +150,34 @@ while ~(cvg || iter > options_.simul.maxit)
        g1aa=g1a;
        ba=b;
        max_resa=max_res;
-        if stack_solve_algo==0
+        if stack_solve_algo==0 || (ismember(stack_solve_algo, [2 3]) && strcmp(options_.simul.preconditioner, 'iterstack') ...
+                                   && options_.simul.iterstack_nperiods == 0 && options_.simul.iterstack_nlu == 0 && size(g1a, 1) < options_.simul.iterstack_maxlu) % Fallback to LU if block too small for iterstack
            dx = g1a\b- ya;
            ya = ya + lambda*dx;
            y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
-        elseif stack_solve_algo==2
+        elseif ismember(stack_solve_algo, [2, 3])
-            flag1=1;
+            if strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0
-            while flag1>0
+                [L, U, P, Q] = lu(g1a);
-                if preconditioner==2
+                zc = L\(P*b);
-                    [L1, U1]=ilu(g1a,ilu_setup);
+                zb = U\zc;
-                elseif preconditioner==3
+            elseif strcmp(options_.simul.preconditioner, 'iterstack')
-                    Size = Blck_size;
+                [L, U, P, Q] = iterstack_preconditioner(g1a, options_);
-                    gss1 =  g1a(Size + 1: 2*Size,Size + 1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
+            elseif strcmp(options_.simul.preconditioner, 'ilu')
-                    [L1, U1]=lu(gss1);
+                [L, U, P] = ilu(g1a, options_.simul.ilu);
-                    L(1:Size,1:Size) = L1;
+                Q = speye(size(g1a));
-                    U(1:Size,1:Size) = U1;
-                    gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
-                    [L2, U2]=lu(gss2);
-                    L(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), L2);
-                    U(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), U2);
-                    gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size);
-                    [L3, U3]=lu(gss2);
-                    L((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = L3;
-                    U((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = U3;
-                    L1 = L;
-                    U1 = U;
-                elseif preconditioner==4
-                    Size = Blck_size;
-                    gss1 =  g1a(1: 3*Size, 1: 3*Size);
-                    [L, U] = lu(gss1);
-                    L1 = kron(eye(ceil(periods/3)),L);
-                    U1 = kron(eye(ceil(periods/3)),U);
-                    L1 = L1(1:periods * Size, 1:periods * Size);
-                    U1 = U1(1:periods * Size, 1:periods * Size);
-                end
-                [za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);
-                if (flag1>0 || reduced)
-                    if verbose
-                        if flag1==1
-                            disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
-                        elseif flag1==2
-                            disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
-                        elseif flag1==3
-                            disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
-                        end
-                    end
-                    ilu_setup.droptol = ilu_setup.droptol/10;
-                    reduced = 0;
-                else
-                    dx = za - ya;
-                    ya = ya + lambda*dx;
-                    y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
-                end
            end
-        elseif stack_solve_algo==3
+            if ~(strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0)
-            flag1=1;
+                [iter_tol, iter_maxit, gmres_restart] = iter_solver_params(options_, g1a, b);
-            while flag1>0
+                if stack_solve_algo == 2
-                if preconditioner==2
+                    [zb, flag] = gmres(P*g1a*Q, P*b, gmres_restart, iter_tol, iter_maxit, L, U);
-                    [L1, U1]=ilu(g1a,ilu_setup);
-                    [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
-                elseif preconditioner==3
-                    Size = Blck_size;
-                    gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
-                    [L1, U1]=lu(gss0);
-                    P1 = eye(size(gss0));
-                    Q1 = eye(size(gss0));
-                    L = kron(eye(periods),L1);
-                    U = kron(eye(periods),U1);
-                    P = kron(eye(periods),P1);
-                    Q = kron(eye(periods),Q1);
-                    [za,flag1] = bicgstab1(g1a,b,1e-7,Blck_size*periods,L,U, P, Q);
                else
-                    Size = Blck_size;
+                    [zb, flag] = bicgstab(P*g1a*Q, P*b, iter_tol, iter_maxit, L, U);
-                    gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
-                    [L1, U1]=lu(gss0);
-                    L1 = kron(eye(periods),L1);
-                    U1 = kron(eye(periods),U1);
-                    [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
-                end
-                if flag1>0 || reduced
-                    if verbose
-                        if flag1==1
-                            disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
-                        elseif flag1==2
-                            disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
-                        elseif flag1==3
-                            disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
-                        end
-                    end
-                    ilu_setup.droptol = ilu_setup.droptol/10;
-                    reduced = 0;
-                else
-                    dx = za - ya;
-                    ya = ya + lambda*dx;
-                    y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
                end
+                iter_solver_error_flag(flag)
            end
+            dx = Q*zb - ya;
+            ya = ya + lambda*dx;
+            y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
        elseif stack_solve_algo==4
            stpmx = 100 ;
            stpmax = stpmx*max([sqrt(ya'*ya);size(y_index,2)]);

--- a/matlab/set_dynare_threads.m
+++ b/matlab/set_dynare_threads.m
 function set_dynare_threads(mexname,n)
+% function set_dynare_threads(mexname,n)
 % This function sets the number of threads used by some MEX files using
 % OpenMP features or any other parallel library.
 %
@@ -8,8 +9,9 @@ function set_dynare_threads(mexname,n)
 %
 % OUTPUTS
 %  none.
+% Documented standalone function, not to be removed
-% Copyright © 2009-2019 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %

--- a/matlab/stochastic_solver/stochastic_solvers.m
+++ b/matlab/stochastic_solver/stochastic_solvers.m
@@ -258,11 +258,12 @@ end
 %exogenous deterministic variables
 if M_.exo_det_nbr > 0
    gx = dr.gx;
-    f1 = g1(:,2*M_.endo_nbr + order_var(nstatic+npred+(1:nsfwrd)));
+    % Algorithm at https://git.dynare.org/Dynare/dynare/-/wikis/varexo_det:-Implementation
-    f0 = g1(:,M_.endo_nbr + order_var(icurr_dr));
+    f1 = g1(:,2*M_.endo_nbr + order_var(nstatic+npred+(1:nsfwrd))); %A matrix, N ny N_forward
-    fudet = g1(:,3*M_.endo_nbr+M_.exo_nbr+1:end);
+    f0 = g1(:,M_.endo_nbr + order_var(icurr_dr)); %B matrix, N by N
-    M1 = inv(f0+[zeros(M_.endo_nbr,nstatic) f1*gx zeros(M_.endo_nbr,nsfwrd-nboth)]);
+    fudet = g1(:,3*M_.endo_nbr+M_.exo_nbr+1:end); %K matrix, N by N_exo_det
-    M2 = M1*f1;
+    M1 = inv(f0+[zeros(M_.endo_nbr,nstatic) f1*gx zeros(M_.endo_nbr,nsfwrd-nboth)]); %(AF+B)^(-1)
+    M2 = M1*f1; %A*(B+AF)^(-1)
    dr.exo_det_length = 0;
    for i = 1:length(M_.det_shocks)
@@ -272,9 +273,9 @@ if M_.exo_det_nbr > 0
    end
    dr.ghud = cell(dr.exo_det_length,1);
-    dr.ghud{1} = -M1*fudet;
+    dr.ghud{1} = -M1*fudet; %-K*(AF+B)^(-1)
    for i = 2:dr.exo_det_length
-        dr.ghud{i} = -M2*dr.ghud{i-1}(end-nsfwrd+1:end,:);
+        dr.ghud{i} = -M2*dr.ghud{i-1}(end-nsfwrd+1:end,:); %-AH_{i-1}*(AF+B)^(-1)
    end
    if local_order > 1

--- a/matlab/subst_auxvar.m
+++ b/matlab/subst_auxvar.m
@@ -59,8 +59,8 @@ if ~isempty(aux_index)
            % Expectation operator
            str = sprintf('%s', M_.aux_vars(aux_index).orig_expr);
            return
-        case 6
+        case {5,6}
-            % Ramsey's multipliers
+            % differentiate_forward_vars and Ramsey multipliers
            if ~isempty(aux_lead_lag)
                str = sprintf('%s(%d)', M_.endo_names{M_.aux_vars(aux_index).endo_index}, aux_lead_lag);
            else

--- a/meson.build
+++ b/meson.build
@@ -173,12 +173,12 @@ if get_option('build_for') == 'matlab'
  umfpack_dep = declare_dependency(link_args : '-lmwumfpack', dependencies : blas_dep)
  ut_dep = declare_dependency(link_args : '-lut')
-  # Workaround for Meson bug https://github.com/mesonbuild/meson/issues/12757
+  if meson.version().version_compare('<1.7.0') and get_option('prefer_static') and host_machine.system() == 'linux'
-  # Use the C compiler as a fallback for detecting SLICOT under Linux with
+    # Workaround for Meson bug https://github.com/mesonbuild/meson/issues/12757 (fixed in 1.7.0)
-  # prefer_static=true (but still try the Fortran compiler to honour the -B
+    # Use the C compiler as a fallback for detecting SLICOT under Linux with
-  # option in fortran_args, as documented). Needed for building the MATLAB
+    # prefer_static=true (but still try the Fortran compiler to honour the -B
-  # Online package.
+    # option in fortran_args, as documented). Needed for building the MATLAB
-  if get_option('prefer_static') and host_machine.system() == 'linux'
+    # Online package.
    slicot_dep_tmp = fortran_compiler.find_library('slicot64_pic', required : false)
    if not slicot_dep_tmp.found()
      slicot_dep_tmp = c_compiler.find_library('slicot64_pic')
@@ -199,7 +199,11 @@ else # Octave build
  octave_incflags = run_command(mkoctfile_exe, '-p', 'INCFLAGS', check : true).stdout().split()
  octlibdir = run_command(mkoctfile_exe, '-p', 'OCTLIBDIR', check : true).stdout().strip()
-  octave_libs = run_command(mkoctfile_exe, '-p', 'OCTAVE_LIBS', check : true).stdout().split()
+  if octave_version.version_compare('<10')
+    octave_libs = run_command(mkoctfile_exe, '-p', 'OCTAVE_LIBS', check : true).stdout().split()
+  else
+    octave_libs = run_command(mkoctfile_exe, '-p', 'LIBOCTMEX', check : true).stdout().split()
+  endif
  octave_link_args = []
@@ -378,10 +382,6 @@ shared_module('logarithmic_reduction', [ 'mex/sources/logarithmic_reduction/mexF
 shared_module('disclyap_fast', [ 'mex/sources/disclyap_fast/disclyap_fast.f08' ] + mex_blas_fortran_iface,
              kwargs : mex_kwargs, dependencies : [ blas_dep, lapack_dep ])
-# TODO: Same remark as A_times_B_kronecker_C
-shared_module('riccati_update', [ 'mex/sources/riccati_update/mexFunction.f08' ] + mex_blas_fortran_iface,
-              kwargs : mex_kwargs, dependencies : [ blas_dep, lapack_dep ])
 qmc_sequence_src = [ 'mex/sources/sobol/qmc_sequence.cc',
                     'mex/sources/sobol/sobol.f08' ]
 # Hack for statically linking libgfortran
@@ -611,7 +611,7 @@ endif
 dynare_manual_html = custom_target('dynare-manual.html',
                                   output : 'dynare-manual.html', input : sphinx_src,
-                                   command : [ sphinx_build_exe, '-b', 'html', sphinx_defines,
+                                   command : [ sphinx_build_exe, '-n', '-b', 'html', sphinx_defines,
                                               '-d', '@PRIVATE_DIR@',
                                               meson.current_source_dir() / 'doc/manual/source',
                                               '@OUTPUT@' ],
@@ -726,6 +726,10 @@ alias_target('doc',
 ### Tests implemented as .mod or .m (both integration and unit tests)
 mod_and_m_tests = [
+  { 'test' : [ 'hank/ks.mod' ],
+    'extra': [ 'hank/ks_ar1_ss.mat',
+               'hank/ks_iid_ss.mat',
+               'hank/expect_error.m' ] },
  { 'test' : [ 'decision_rules/first_order/walsh.mod' ],
    'extra' : [ 'decision_rules/first_order/+matlab/+namespace/y_k.m'] },
  { 'test' : [ 'occbin/model_irrcap_twoconstraints/dynrbc.mod',
@@ -823,6 +827,7 @@ mod_and_m_tests = [
  { 'test' : [ 'estimation/slice/fs2000_slice.mod' ],
    'extra' : [ 'estimation/fsdat_simul.m' ] },
  { 'test' : [ 'analytic_derivatives/fs2000_analytic_derivation.mod' ] },
+  { 'test' : [ 'analytic_derivatives/fs2000_analytic_derivation_newrat.mod' ] },
  { 'test' : [ 'analytic_derivatives/BrockMirman_PertParamsDerivs.mod' ],
    'extra' : [ 'analytic_derivatives/nBrockMirmanSYM.mat' ] },
  { 'test' : [ 'analytic_derivatives/burnside_3_order_PertParamsDerivs.mod' ] },
@@ -865,7 +870,6 @@ mod_and_m_tests = [
  { 'test' : [ 'fs2000/fs2000.mod' ],
    'extra' : [ 'fs2000/fsdat_simul.m' ] },
  { 'test' : [ 'ls2003/ls2003_hessian_zero.mod' ] },
-  { 'test' : [ 'ep/rbc.mod' ] },
  { 'test' : [ 'exogenous-observed-variables/preprocessor.mod' ] },
  { 'test' : [ 'estimation/fs2000_with_weibull_prior.mod' ],
    'extra' : [ 'estimation/fsdat_simul.m' ] },
@@ -1069,6 +1073,7 @@ mod_and_m_tests = [
               'discretionary_policy/dennis_1_estim.mod',
               'discretionary_policy/dennis_1_estim_MoM.mod' ] },
  { 'test' : [ 'discretionary_policy/Gali_discretion.mod' ] },
+  { 'test' : [ 'discretionary_policy/NK_discretion_forward.mod' ] },
  { 'test' : [ 'discretionary_policy/Gali_2015_chapter_3.mod',
               'discretionary_policy/Gali_2015_chapter_3_nonlinear.mod' ] },
  { 'test' : [ 'histval_initval_file/ramst_data_generate.mod',
@@ -1366,6 +1371,7 @@ mod_and_m_tests = [
  { 'test' : [ 'stochastic-backward-models/solow_ces.mod' ] },
  { 'test' : [ 'stochastic-backward-models/solow_cd_with_steadystate.mod' ] },
  { 'test' : [ 'deterministic_simulations/ramst.mod' ] },
+  { 'test' : [ 'deterministic_simulations/ramst_dseries.mod' ] },
  { 'test' : [ 'deterministic_simulations/ramst_a.mod' ] },
  { 'test' : [ 'deterministic_simulations/ramst_mshocks.mod' ] },
  { 'test' : [ 'deterministic_simulations/ramst_mshocks_vec.mod' ] },
@@ -1379,6 +1385,7 @@ mod_and_m_tests = [
  { 'test' : [ 'deterministic_simulations/Irreversible_investment.mod',
               'deterministic_simulations/Irreversible_investment_leadlag.mod' ] },
  { 'test' : [ 'deterministic_simulations/linear_state_space_arma.mod' ] },
+  { 'test' : [ 'deterministic_simulations/linear_state_space_arma_condpaths.mod' ] },
  { 'test' : [ 'deterministic_simulations/purely_forward/ar1.mod' ] },
  { 'test' : [ 'deterministic_simulations/purely_forward/nk.mod' ] },
  { 'test' : [ 'deterministic_simulations/purely_backward/ar1.mod' ] },
@@ -1393,8 +1400,15 @@ mod_and_m_tests = [
  { 'test' : [ 'deterministic_simulations/rbc_det6.mod' ] },
  { 'test' : [ 'deterministic_simulations/homotopy.mod' ] },
  { 'test' : [ 'deterministic_simulations/homotopy_histval.mod' ] },
-  { 'test' : [ 'deterministic_simulations/homotopy_endval_steady_linearization.mod' ] },
+  { 'test' : [
-  { 'test' : [ 'deterministic_simulations/homotopy_marginal_linearization.mod' ] },
+    'deterministic_simulations/homotopy_endval_steady_linearization.mod' ],
+    'extra' : [ 'deterministic_simulations/homotopy_linearization.inc' ] },
+  { 'test' : [ 'deterministic_simulations/homotopy_marginal_linearization.mod' ],
+    'extra' : [ 'deterministic_simulations/homotopy_linearization.inc' ] },
+  { 'test' : [ 'deterministic_simulations/homotopy_linearization_condpaths.mod' ],
+    'extra' : [ 'deterministic_simulations/homotopy_linearization.inc' ] },
+  { 'test' : [ 'deterministic_simulations/homotopy_marglinearization_condpaths.mod' ],
+    'extra' : [ 'deterministic_simulations/homotopy_linearization.inc' ] },
  { 'test' : [ 'deterministic_simulations/homotopy_exclude_varexo.mod' ] },
  { 'test' : [ 'deterministic_simulations/rbc_det_exo_lag_2a.mod',
               'deterministic_simulations/rbc_det_exo_lag_2b.mod',
@@ -1420,10 +1434,24 @@ mod_and_m_tests = [
    'extra' : [ 'deterministic_simulations/pfwee.csv' ] },
  { 'test' : [ 'deterministic_simulations/pfwee_constant_sim_length.mod' ],
    'extra' : [ 'deterministic_simulations/pfwee.csv' ] },
-  { 'test' : [ 'deterministic_simulations/pfwee_learnt_in.mod' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_learnt_in.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_learnt_in.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_learnt_in_bare_1st_period.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_learnt_in.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_learnt_in_dseries.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_learnt_in.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_learnt_in_dseries_bare_1st_period.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_learnt_in.inc' ] },
  { 'test' : [ 'deterministic_simulations/pfwee_multiple_shocks.mod' ] },
-  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_linearization.mod' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_linearization.mod' ],
-  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_marginal_linearization.mod' ] },
+    'extra' : [ 'deterministic_simulations/pfwee_homotopy.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_marginal_linearization.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_homotopy.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_linearization_condpaths.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_homotopy.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_marglin_condpaths.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_homotopy.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_marginal_linearization.mod' ] },
  { 'test' : [ 'deterministic_simulations/lbj/rbc.mod' ] },
  { 'test' : [ 'lmmcp/rbc.mod' ] },
  { 'test' : [ 'lmmcp/sw_lmmcp.mod',
@@ -1599,6 +1627,13 @@ mod_and_m_tests = [
               'deprecated/ramsey_ex.mod' ] },
  { 'test' : [ 'deterministic_simulations/ramst.mod',
               'deprecated/ramst.mod' ] },
+  { 'test' : [ 'deterministic_simulations/ramst_pf_controlled_paths.mod',
+               'deterministic_simulations/ramst_pf_controlled_paths_lbj.mod',
+               'deterministic_simulations/ramst_pf_controlled_paths_ssa7.mod',
+               'deterministic_simulations/ramst_pf_controlled_paths_check.mod' ],
+    'extra' : [ 'deterministic_simulations/ramst_pf_controlled_paths_common.inc'] },
+  { 'test' : [ 'deterministic_simulations/ramst_pfwee_controlled_paths.mod' ],
+    'extra' : [ 'deterministic_simulations/ramst_pf_controlled_paths_common.inc'] },
  # ECB modfiles
  { 'test' : [ 'var-expectations/1/example1.mod' ] },
@@ -1764,6 +1799,7 @@ mod_and_m_tests = [
                'particle/benchmark.m',
                'particle/mysample.m'] },
  { 'test' : [ 'particle/first_spec.mod',
+               'particle/first_spec_hssmc.mod',
               'particle/local_state_space_iteration_k_test.mod',
               'particle/local_state_space_iteration_3_test.mod' ],
    'extra' : [ 'particle/first_spec_common.inc' ] },
@@ -1863,7 +1899,6 @@ mod_and_m_tests = [
                'solver-test-functions/wood.m',] },
  { 'test' : [ 'cyclereduction.m' ] },
  { 'test' : [ 'logarithmicreduction.m' ] },
-  { 'test' : [ 'riccatiupdate.m' ] },
  { 'test' : [ 'kalman/likelihood/test_kalman_mex.m' ] },
  { 'test' : [ 'contribs.m' ],
    'extra' : [ 'sandbox.mod',

--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -30,6 +30,54 @@
 #include "Interpreter.hh"
+void
+iter_solver_opts_t::set_static_values(const mxArray* options_)
+{
+  // Should be the same options as in newton_solve.m and solve_one_boundary.m (static case)
+  iter_tol = mxCreateDoubleMatrix(0, 0, mxREAL);
+  iter_maxit = mxCreateDoubleMatrix(0, 0, mxREAL);
+  gmres_restart = mxCreateDoubleMatrix(0, 0, mxREAL);
+  int field {mxGetFieldNumber(options_, "steady")};
+  if (field < 0)
+    mexErrMsgTxt("steady is not a field of options_");
+  mxArray* steady {mxGetFieldByNumber(options_, 0, field)};
+  field = mxGetFieldNumber(steady, "ilu");
+  if (field < 0)
+    mexErrMsgTxt("ilu is not a field of options_.steady");
+  ilu = mxGetFieldByNumber(steady, 0, field);
+}
+void
+iter_solver_opts_t::set_dynamic_values(const mxArray* options_)
+{
+  int field {mxGetFieldNumber(options_, "simul")};
+  if (field < 0)
+    mexErrMsgTxt("simul is not a field of options_");
+  mxArray* simul {mxGetFieldByNumber(options_, 0, field)};
+  field = mxGetFieldNumber(simul, "iter_tol");
+  if (field < 0)
+    mexErrMsgTxt("iter_tol is not a field of options_.simul");
+  iter_tol = mxGetFieldByNumber(simul, 0, field);
+  field = mxGetFieldNumber(simul, "iter_maxit");
+  if (field < 0)
+    mexErrMsgTxt("iter_maxit is not a field of options_.simul");
+  iter_maxit = mxGetFieldByNumber(simul, 0, field);
+  field = mxGetFieldNumber(simul, "gmres_restart");
+  if (field < 0)
+    mexErrMsgTxt("gmres_restart is not a field of options_.simul");
+  gmres_restart = mxGetFieldByNumber(simul, 0, field);
+  field = mxGetFieldNumber(simul, "ilu");
+  if (field < 0)
+    mexErrMsgTxt("ilu is not a field of options_.simul");
+  ilu = mxGetFieldByNumber(simul, 0, field);
+}
 Interpreter::Interpreter(Evaluate& evaluator_arg, double* params_arg, double* y_arg, double* ya_arg,
                         double* x_arg, double* steady_y_arg, double* direction_arg, int y_size_arg,
                         int nb_row_x_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
@@ -38,11 +86,13 @@ Interpreter::Interpreter(Evaluate& evaluator_arg, double* params_arg, double* y_
                         int solve_algo_arg, bool print_arg,
                         const mxArray* GlobalTemporaryTerms_arg, bool steady_state_arg,
                         bool block_decomposed_arg, int col_x_arg, int col_y_arg,
-                         const BasicSymbolTable& symbol_table_arg, int verbosity_arg) :
+                         const BasicSymbolTable& symbol_table_arg, int verbosity_arg,
+                         iter_solver_opts_t iter_solver_opts_arg) :
    symbol_table {symbol_table_arg},
    steady_state {steady_state_arg},
    block_decomposed {block_decomposed_arg},
    evaluator {evaluator_arg},
+    iter_solver_opts {iter_solver_opts_arg},
    minimal_solving_periods {minimal_solving_periods_arg},
    y_size {y_size_arg},
    y_kmin {y_kmin_arg},
@@ -60,7 +110,6 @@ Interpreter::Interpreter(Evaluate& evaluator_arg, double* params_arg, double* y_
  start_compare = 0;
  restart = 0;
  IM_i.clear();
-  lu_inc_tol = 1e-10;
  Symbolic = nullptr;
  Numeric = nullptr;
  params = params_arg;
@@ -2198,7 +2247,7 @@ Interpreter::complete(int beg_t)
    }
  mxFree(save_code);
  mxFree(diff);
-  return (beg_t);
+  return beg_t;
 }
 void
@@ -2260,263 +2309,6 @@ Interpreter::simple_bksub()
    }
 }
-mxArray*
-Interpreter::subtract_A_B(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  size_t m_A = mxGetM(A_m);
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateDoubleMatrix(m_A, n_B, mxREAL);
-  double* C_d = mxGetPr(C_m);
-  for (int j = 0; j < static_cast<int>(n_A); j++)
-    for (unsigned int i = 0; i < m_A; i++)
-      {
-        size_t index = j * m_A + i;
-        C_d[index] = A_d[index] - B_d[index];
-      }
-  return C_m;
-}
-mxArray*
-Interpreter::Sparse_subtract_SA_SB(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  size_t m_A = mxGetM(A_m);
-  mwIndex* A_i = mxGetIr(A_m);
-  mwIndex* A_j = mxGetJc(A_m);
-  size_t total_nze_A = A_j[n_A];
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  mwIndex* B_i = mxGetIr(B_m);
-  mwIndex* B_j = mxGetJc(B_m);
-  size_t total_nze_B = B_j[n_B];
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateSparse(m_A, n_B, m_A * n_B, mxREAL);
-  mwIndex* C_i = mxGetIr(C_m);
-  mwIndex* C_j = mxGetJc(C_m);
-  double* C_d = mxGetPr(C_m);
-  unsigned int nze_B = 0, nze_C = 0, nze_A = 0;
-  unsigned int A_col = 0, B_col = 0, C_col = 0;
-  C_j[C_col] = 0;
-  while (nze_A < total_nze_A || nze_B < total_nze_B)
-    {
-      while (nze_A >= static_cast<unsigned int>(A_j[A_col + 1]) && (nze_A < total_nze_A))
-        A_col++;
-      size_t A_row = A_i[nze_A];
-      while (nze_B >= static_cast<unsigned int>(B_j[B_col + 1]) && (nze_B < total_nze_B))
-        B_col++;
-      size_t B_row = B_i[nze_B];
-      if (A_col == B_col)
-        {
-          if (A_row == B_row && (nze_B < total_nze_B && nze_A < total_nze_A))
-            {
-              C_d[nze_C] = A_d[nze_A++] - B_d[nze_B++];
-              C_i[nze_C] = A_row;
-              while (C_col < A_col)
-                C_j[++C_col] = nze_C;
-              C_j[A_col + 1] = nze_C++;
-              C_col = A_col;
-            }
-          else if ((A_row < B_row && nze_A < total_nze_A) || nze_B == total_nze_B)
-            {
-              C_d[nze_C] = A_d[nze_A++];
-              C_i[nze_C] = A_row;
-              while (C_col < A_col)
-                C_j[++C_col] = nze_C;
-              C_j[A_col + 1] = nze_C++;
-              C_col = A_col;
-            }
-          else
-            {
-              C_d[nze_C] = -B_d[nze_B++];
-              C_i[nze_C] = B_row;
-              while (C_col < B_col)
-                C_j[++C_col] = nze_C;
-              C_j[B_col + 1] = nze_C++;
-              C_col = B_col;
-            }
-        }
-      else if ((A_col < B_col && nze_A < total_nze_A) || nze_B == total_nze_B)
-        {
-          C_d[nze_C] = A_d[nze_A++];
-          C_i[nze_C] = A_row;
-          while (C_col < A_col)
-            C_j[++C_col] = nze_C;
-          C_j[A_col + 1] = nze_C++;
-          C_col = A_col;
-        }
-      else
-        {
-          C_d[nze_C] = -B_d[nze_B++];
-          C_i[nze_C] = B_row;
-          while (C_col < B_col)
-            C_j[++C_col] = nze_C;
-          C_j[B_col + 1] = nze_C++;
-          C_col = B_col;
-        }
-    }
-  while (C_col < n_B)
-    C_j[++C_col] = nze_C;
-  mxSetNzmax(C_m, nze_C);
-  return C_m;
-}
-mxArray*
-Interpreter::mult_SAT_B(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  mwIndex* A_i = mxGetIr(A_m);
-  mwIndex* A_j = mxGetJc(A_m);
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateDoubleMatrix(n_A, n_B, mxREAL);
-  double* C_d = mxGetPr(C_m);
-  for (int j = 0; j < static_cast<int>(n_B); j++)
-    for (unsigned int i = 0; i < n_A; i++)
-      {
-        double sum = 0;
-        size_t nze_A = A_j[i];
-        while (nze_A < static_cast<unsigned int>(A_j[i + 1]))
-          {
-            size_t i_A = A_i[nze_A];
-            sum += A_d[nze_A++] * B_d[i_A];
-          }
-        C_d[j * n_A + i] = sum;
-      }
-  return C_m;
-}
-mxArray*
-Interpreter::Sparse_mult_SAT_B(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  mwIndex* A_i = mxGetIr(A_m);
-  mwIndex* A_j = mxGetJc(A_m);
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  size_t m_B = mxGetM(B_m);
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateSparse(n_A, n_B, n_A * n_B, mxREAL);
-  mwIndex* C_i = mxGetIr(C_m);
-  mwIndex* C_j = mxGetJc(C_m);
-  double* C_d = mxGetPr(C_m);
-  unsigned int nze_C = 0;
-  // unsigned int nze_A = 0;
-  unsigned int C_col = 0;
-  C_j[C_col] = 0;
-  // #pragma omp parallel for
-  for (unsigned int j = 0; j < n_B; j++)
-    for (unsigned int i = 0; i < n_A; i++)
-      {
-        double sum = 0;
-        size_t nze_A = A_j[i];
-        while (nze_A < static_cast<unsigned int>(A_j[i + 1]))
-          {
-            size_t i_A = A_i[nze_A];
-            sum += A_d[nze_A++] * B_d[i_A];
-          }
-        if (fabs(sum) > 1e-10)
-          {
-            C_d[nze_C] = sum;
-            C_i[nze_C] = i;
-            while (C_col < j)
-              C_j[++C_col] = nze_C;
-            nze_C++;
-          }
-      }
-  while (C_col < m_B)
-    C_j[++C_col] = nze_C;
-  mxSetNzmax(C_m, nze_C);
-  return C_m;
-}
-mxArray*
-Interpreter::Sparse_mult_SAT_SB(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  mwIndex* A_i = mxGetIr(A_m);
-  mwIndex* A_j = mxGetJc(A_m);
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  mwIndex* B_i = mxGetIr(B_m);
-  mwIndex* B_j = mxGetJc(B_m);
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateSparse(n_A, n_B, n_A * n_B, mxREAL);
-  mwIndex* C_i = mxGetIr(C_m);
-  mwIndex* C_j = mxGetJc(C_m);
-  double* C_d = mxGetPr(C_m);
-  size_t nze_B = 0, nze_C = 0, nze_A = 0;
-  unsigned int C_col = 0;
-  C_j[C_col] = 0;
-  for (unsigned int j = 0; j < n_B; j++)
-    for (unsigned int i = 0; i < n_A; i++)
-      {
-        double sum = 0;
-        nze_B = B_j[j];
-        nze_A = A_j[i];
-        while (nze_A < static_cast<unsigned int>(A_j[i + 1])
-               && nze_B < static_cast<unsigned int>(B_j[j + 1]))
-          {
-            size_t i_A = A_i[nze_A];
-            size_t i_B = B_i[nze_B];
-            if (i_A == i_B)
-              sum += A_d[nze_A++] * B_d[nze_B++];
-            else if (i_A < i_B)
-              nze_A++;
-            else
-              nze_B++;
-          }
-        if (fabs(sum) > 1e-10)
-          {
-            C_d[nze_C] = sum;
-            C_i[nze_C] = i;
-            while (C_col < j)
-              C_j[++C_col] = nze_C;
-            nze_C++;
-          }
-      }
-  while (C_col < n_B)
-    C_j[++C_col] = nze_C;
-  mxSetNzmax(C_m, nze_C);
-  return C_m;
-}
-mxArray*
-Interpreter::Sparse_transpose(const mxArray* A_m)
-{
-  size_t n_A = mxGetN(A_m);
-  size_t m_A = mxGetM(A_m);
-  mwIndex* A_i = mxGetIr(A_m);
-  mwIndex* A_j = mxGetJc(A_m);
-  size_t total_nze_A = A_j[n_A];
-  double* A_d = mxGetPr(A_m);
-  mxArray* C_m = mxCreateSparse(n_A, m_A, total_nze_A, mxREAL);
-  mwIndex* C_i = mxGetIr(C_m);
-  mwIndex* C_j = mxGetJc(C_m);
-  double* C_d = mxGetPr(C_m);
-  unsigned int nze_C = 0, nze_A = 0;
-  ranges::fill_n(C_j, m_A + 1, 0);
-  map<pair<mwIndex, unsigned int>, double> B2;
-  for (unsigned int i = 0; i < n_A; i++)
-    while (nze_A < static_cast<unsigned int>(A_j[i + 1]))
-      {
-        C_j[A_i[nze_A] + 1]++;
-        B2[{A_i[nze_A], i}] = A_d[nze_A];
-        nze_A++;
-      }
-  for (unsigned int i = 0; i < m_A; i++)
-    C_j[i + 1] += C_j[i];
-  for (auto& [key, val] : B2)
-    {
-      C_d[nze_C] = val;
-      C_i[nze_C++] = key.second;
-    }
-  return C_m;
-}
 void
 Interpreter::compute_block_time(int my_Per_u_, bool evaluate, bool no_derivatives)
 {
@@ -2948,24 +2740,19 @@ void
 Interpreter::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, bool is_two_boundaries, mxArray* x0_m)
 {
  size_t n = mxGetM(A_m);
-  std::array field_names {"droptol", "type"};
-  std::array dims {static_cast<mwSize>(1)};
-  mxArray* Setup
-      = mxCreateStructArray(dims.size(), dims.data(), field_names.size(), field_names.data());
-  mxSetFieldByNumber(Setup, 0, 0, mxCreateDoubleScalar(lu_inc_tol));
-  mxSetFieldByNumber(Setup, 0, 1, mxCreateString("ilutp"));
  std::array<mxArray*, 2> lhs0;
-  std::array rhs0 {A_m, Setup};
+  std::array rhs0 {A_m, iter_solver_opts.ilu};
  if (mexCallMATLAB(lhs0.size(), lhs0.data(), rhs0.size(), rhs0.data(), "ilu"))
    throw FatalException("In GMRES, the incomplete LU decomposition (ilu) has failed");
  mxArray* L1 = lhs0[0];
  mxArray* U1 = lhs0[1];
-  /*[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);*/
  std::array rhs {A_m,
                  b_m,
-                  mxCreateDoubleScalar(size),
+                  iter_solver_opts.gmres_restart,
-                  mxCreateDoubleScalar(1e-6),
+                  iter_solver_opts.iter_tol,
-                  mxCreateDoubleScalar(static_cast<double>(n)),
+                  iter_solver_opts.iter_maxit,
                  L1,
                  U1,
                  x0_m};
@@ -2974,28 +2761,24 @@ Interpreter::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, bool is_two_boundari
  mxArray* z = lhs[0];
  mxArray* flag = lhs[1];
  double flag1 {mxGetScalar(flag)};
-  mxDestroyArray(rhs0[1]);
-  mxDestroyArray(rhs[2]);
-  mxDestroyArray(rhs[3]);
-  mxDestroyArray(rhs[4]);
  mxDestroyArray(rhs[5]);
  mxDestroyArray(rhs[6]);
  if (flag1 > 0)
    {
      if (flag1 == 1)
-        mexWarnMsgTxt(
+        mexErrMsgTxt(
            ("Error in bytecode: No convergence inside GMRES, in block " + to_string(block_num + 1))
                .c_str());
      else if (flag1 == 2)
-        mexWarnMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
+        mexErrMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
-                       + to_string(block_num + 1))
+                      + to_string(block_num + 1))
-                          .c_str());
+                         .c_str());
      else if (flag1 == 3)
-        mexWarnMsgTxt(("Error in bytecode: GMRES stagnated (Two consecutive iterates were the "
+        mexErrMsgTxt(("Error in bytecode: GMRES stagnated (Two consecutive iterates were the "
-                       "same.), in block "
+                      "same.), in block "
-                       + to_string(block_num + 1))
+                      + to_string(block_num + 1))
-                          .c_str());
+                         .c_str());
-      lu_inc_tol /= 10;
    }
  else
    {
@@ -3017,6 +2800,7 @@ Interpreter::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, bool is_two_boundari
            y[eq + it_ * y_size] += slowc * yy;
          }
    }
  mxDestroyArray(A_m);
  mxDestroyArray(b_m);
  mxDestroyArray(x0_m);
@@ -3026,143 +2810,42 @@ Interpreter::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, bool is_two_boundari
 void
 Interpreter::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, bool is_two_boundaries,
-                                   mxArray* x0_m, int preconditioner)
+                                   mxArray* x0_m)
 {
-  /* precond = 0  => Jacobi
-     precond = 1  => Incomplet LU decomposition*/
  size_t n = mxGetM(A_m);
-  mxArray *L1 = nullptr, *U1 = nullptr, *Diag = nullptr;
-  if (preconditioner == 0)
+  std::array<mxArray*, 2> lhs0;
-    {
+  std::array rhs0 {A_m, iter_solver_opts.ilu};
-      std::array<mxArray*, 1> lhs0;
+  if (mexCallMATLAB(lhs0.size(), lhs0.data(), rhs0.size(), rhs0.data(), "ilu"))
-      std::array rhs0 {A_m, mxCreateDoubleScalar(0)};
+    throw FatalException {"In BiCGStab, the incomplete LU decomposition (ilu) has failed"};
-      mexCallMATLAB(lhs0.size(), lhs0.data(), rhs0.size(), rhs0.data(), "spdiags");
+  mxArray* L1 = lhs0[0];
-      mxArray* tmp = lhs0[0];
+  mxArray* U1 = lhs0[1];
-      double* tmp_val = mxGetPr(tmp);
-      Diag = mxCreateSparse(n, n, n, mxREAL);
-      mwIndex* Diag_i = mxGetIr(Diag);
-      mwIndex* Diag_j = mxGetJc(Diag);
-      double* Diag_val = mxGetPr(Diag);
-      for (size_t i = 0; i < n; i++)
-        {
-          Diag_val[i] = tmp_val[i];
-          Diag_j[i] = i;
-          Diag_i[i] = i;
-        }
-      Diag_j[n] = n;
-    }
-  else if (preconditioner == 1)
-    {
-      /*[L1, U1] = ilu(g1a=;*/
-      std::array field_names {"type", "droptol"};
-      const int type {0}, droptol {1};
-      std::array dims {static_cast<mwSize>(1)};
-      mxArray* Setup
-          = mxCreateStructArray(dims.size(), dims.data(), field_names.size(), field_names.data());
-      mxSetFieldByNumber(Setup, 0, type, mxCreateString("ilutp"));
-      mxSetFieldByNumber(Setup, 0, droptol, mxCreateDoubleScalar(lu_inc_tol));
-      std::array<mxArray*, 2> lhs0;
-      std::array rhs0 {A_m, Setup};
-      if (mexCallMATLAB(lhs0.size(), lhs0.data(), rhs0.size(), rhs0.data(), "ilu"))
-        throw FatalException {"In BiCGStab, the incomplete LU decomposition (ilu) has failed"};
-      L1 = lhs0[0];
-      U1 = lhs0[1];
-      mxDestroyArray(Setup);
-    }
-  double flags = 2;
-  mxArray* z = nullptr;
-  if (steady_state) /*Octave BicStab algorihtm involves a 0 division in case of a preconditionner
-                       equal to the LU decomposition of A matrix*/
-    {
-      mxArray* res = mult_SAT_B(Sparse_transpose(A_m), x0_m);
-      double* resid = mxGetPr(res);
-      double* b = mxGetPr(b_m);
-      for (int i = 0; i < static_cast<int>(n); i++)
-        resid[i] = b[i] - resid[i];
-      std::array<mxArray*, 1> lhs;
-      std::array rhs {L1, res};
-      mexCallMATLAB(lhs.size(), lhs.data(), rhs.size(), rhs.data(), "mldivide");
-      std::array rhs2 {U1, lhs[0]};
-      mexCallMATLAB(lhs.size(), lhs.data(), rhs2.size(), rhs2.data(), "mldivide");
-      z = lhs[0];
-      double* phat = mxGetPr(z);
-      double* x0 = mxGetPr(x0_m);
-      for (int i = 0; i < static_cast<int>(n); i++)
-        phat[i] = x0[i] + phat[i];
-      /*Check the solution*/
-      res = mult_SAT_B(Sparse_transpose(A_m), z);
-      resid = mxGetPr(res);
-      double cum_abs = 0;
-      for (int i = 0; i < static_cast<int>(n); i++)
-        {
-          resid[i] = b[i] - resid[i];
-          cum_abs += fabs(resid[i]);
-        }
-      if (cum_abs > 1e-7)
-        flags = 2;
-      else
-        flags = 0;
-      mxDestroyArray(res);
-    }
-  if (flags == 2)
+  std::array rhs {A_m, b_m, iter_solver_opts.iter_tol, iter_solver_opts.iter_maxit, L1, U1, x0_m};
-    {
+  std::array<mxArray*, 2> lhs;
-      if (preconditioner == 0)
+  mexCallMATLAB(lhs.size(), lhs.data(), rhs.size(), rhs.data(), "bicgstab");
-        {
+  mxArray* z = lhs[0];
-          /*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
+  mxArray* flag = lhs[1];
-          std::array rhs {A_m, b_m, mxCreateDoubleScalar(1e-6),
+  double flags {mxGetScalar(flag)};
-                          mxCreateDoubleScalar(static_cast<double>(n)), Diag};
+  mxDestroyArray(flag);
-          std::array<mxArray*, 2> lhs;
+  mxDestroyArray(rhs[4]);
-          mexCallMATLAB(lhs.size(), lhs.data(), rhs.size(), rhs.data(), "bicgstab");
+  mxDestroyArray(rhs[5]);
-          z = lhs[0];
-          mxArray* flag = lhs[1];
-          flags = mxGetScalar(flag);
-          mxDestroyArray(flag);
-          mxDestroyArray(rhs[2]);
-          mxDestroyArray(rhs[3]);
-          mxDestroyArray(rhs[4]);
-        }
-      else if (preconditioner == 1)
-        {
-          /*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
-          std::array rhs {A_m,
-                          b_m,
-                          mxCreateDoubleScalar(1e-6),
-                          mxCreateDoubleScalar(static_cast<double>(n)),
-                          L1,
-                          U1,
-                          x0_m};
-          std::array<mxArray*, 2> lhs;
-          mexCallMATLAB(lhs.size(), lhs.data(), rhs.size(), rhs.data(), "bicgstab");
-          z = lhs[0];
-          mxArray* flag = lhs[1];
-          flags = mxGetScalar(flag);
-          mxDestroyArray(flag);
-          mxDestroyArray(rhs[2]);
-          mxDestroyArray(rhs[3]);
-          mxDestroyArray(rhs[4]);
-          mxDestroyArray(rhs[5]);
-        }
-    }
  if (flags > 0)
    {
      if (flags == 1)
-        mexWarnMsgTxt(("Error in bytecode: No convergence inside BiCGStab, in block "
+        mexErrMsgTxt(("Error in bytecode: No convergence inside BiCGStab, in block "
-                       + to_string(block_num + 1))
+                      + to_string(block_num + 1))
-                          .c_str());
+                         .c_str());
      else if (flags == 2)
-        mexWarnMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
+        mexErrMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
-                       + to_string(block_num + 1))
+                      + to_string(block_num + 1))
-                          .c_str());
+                         .c_str());
      else if (flags == 3)
-        mexWarnMsgTxt(("Error in bytecode: BiCGStab stagnated (Two consecutive iterates were the "
+        mexErrMsgTxt(("Error in bytecode: BiCGStab stagnated (Two consecutive iterates were the "
-                       "same.), in block "
+                      "same.), in block "
-                       + to_string(block_num + 1))
+                      + to_string(block_num + 1))
-                          .c_str());
+                         .c_str());
-      lu_inc_tol /= 10;
    }
  else
    {
@@ -3184,6 +2867,7 @@ Interpreter::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, bool is_two_bound
            y[eq + it_ * y_size] += slowc * yy;
          }
    }
  mxDestroyArray(A_m);
  mxDestroyArray(b_m);
  mxDestroyArray(x0_m);
@@ -4139,7 +3823,6 @@ Interpreter::Simulate_One_Boundary()
  mxArray *b_m = nullptr, *A_m = nullptr, *x0_m = nullptr;
  SuiteSparse_long *Ap = nullptr, *Ai = nullptr;
  double *Ax = nullptr, *b = nullptr;
-  int preconditioner = 1;
  try_at_iteration = 0;
  Clear_u();
@@ -4208,14 +3891,10 @@ Interpreter::Simulate_One_Boundary()
              mexPrintf("MODEL STEADY STATE: (method=Sparse LU)\n");
              break;
            case 7:
-              mexPrintf(preconditioner_print_out("MODEL STEADY STATE: (method=GMRES)\n",
+              mexPrintf("MODEL STEADY STATE: (method=GMRES with 'ilu' preconditioner)\n");
-                                                 preconditioner, true)
-                            .c_str());
              break;
            case 8:
-              mexPrintf(preconditioner_print_out("MODEL STEADY STATE: (method=BiCGStab)\n",
+              mexPrintf("MODEL STEADY STATE: (method=BiCGStab with 'ilu' preconditioner)\n");
-                                                 preconditioner, true)
-                            .c_str());
              break;
            }
        }
@@ -4285,7 +3964,7 @@ Interpreter::Simulate_One_Boundary()
      else if ((solve_algo == 7 && steady_state) || (stack_solve_algo == 2 && !steady_state))
        Solve_Matlab_GMRES(A_m, b_m, false, x0_m);
      else if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 3 && !steady_state))
-        Solve_Matlab_BiCGStab(A_m, b_m, false, x0_m, preconditioner);
+        Solve_Matlab_BiCGStab(A_m, b_m, false, x0_m);
      else if ((solve_algo == 6 && steady_state)
               || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4
                    || stack_solve_algo == 6)
@@ -4402,43 +4081,12 @@ Interpreter::Simulate_Newton_One_Boundary(bool forward)
  mxFree(r);
 }
-string
-Interpreter::preconditioner_print_out(string s, int preconditioner, bool ss)
-{
-  int n = s.length();
-  string tmp = ", preconditioner=";
-  switch (preconditioner)
-    {
-    case 0:
-      if (ss)
-        tmp.append("Jacobi on static jacobian");
-      else
-        tmp.append("Jacobi on dynamic jacobian");
-      break;
-    case 1:
-      if (ss)
-        tmp.append("incomplete lutp on static jacobian");
-      else
-        tmp.append("incomplete lu0 on dynamic jacobian");
-      break;
-    case 2:
-      tmp.append("incomplete lutp on dynamic jacobian");
-      break;
-    case 3:
-      tmp.append("lu on static jacobian");
-      break;
-    }
-  s.insert(n - 2, tmp);
-  return s;
-}
 void
 Interpreter::Simulate_Newton_Two_Boundaries(
    bool cvg, const vector_table_conditional_local_type& vector_table_conditional_local)
 {
  double top = 0.5;
  double bottom = 0.1;
-  int preconditioner = 2;
  if (start_compare == 0)
    start_compare = y_kmin;
  u_count_alloc_save = u_count_alloc;
@@ -4592,15 +4240,11 @@ Interpreter::Simulate_Newton_Two_Boundaries(
              break;
            case 2:
              mexPrintf(
-                  preconditioner_print_out("MODEL SIMULATION: (method=GMRES on stacked system)\n",
+                  "MODEL SIMULATION: (method=GMRES on stacked system with 'ilu' preconditioner)\n");
-                                           preconditioner, false)
-                      .c_str());
              break;
            case 3:
-              mexPrintf(preconditioner_print_out(
+              mexPrintf("MODEL SIMULATION: (method=BiCGStab on stacked system with 'ilu' "
-                            "MODEL SIMULATION: (method=BiCGStab on stacked system)\n",
+                        "preconditioner)\n");
-                            preconditioner, false)
-                            .c_str());
              break;
            case 4:
              mexPrintf("MODEL SIMULATION: (method=Sparse LU solver with optimal path length on "
@@ -4658,7 +4302,7 @@ Interpreter::Simulate_Newton_Two_Boundaries(
      else if (stack_solve_algo == 2)
        Solve_Matlab_GMRES(A_m, b_m, true, x0_m);
      else if (stack_solve_algo == 3)
-        Solve_Matlab_BiCGStab(A_m, b_m, true, x0_m, 1);
+        Solve_Matlab_BiCGStab(A_m, b_m, true, x0_m);
      else if (stack_solve_algo == 5)
        Solve_ByteCode_Symbolic_Sparse_GaussianElimination(symbolic);
    }

--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -69,6 +69,19 @@ constexpr double mem_increasing_factor = 1.1;
 constexpr int NO_ERROR_ON_EXIT {0}, ERROR_ON_EXIT {1};
+// Stores various options for iterative solvers (GMRES, BiCGStab), from options_.simul
+struct iter_solver_opts_t
+{
+  const char* preconditioner;
+  mxArray* iter_tol;
+  mxArray* iter_maxit;
+  mxArray* gmres_restart;
+  mxArray* ilu;
+  void set_static_values(const mxArray* options_);
+  void set_dynamic_values(const mxArray* options_);
+};
 class Interpreter
 {
 private:
@@ -110,7 +123,7 @@ private:
  int iter;
  int start_compare;
  int restart;
-  double lu_inc_tol;
+  iter_solver_opts_t iter_solver_opts;
  SuiteSparse_long *Ap_save, *Ai_save;
  double *Ax_save, *b_save;
@@ -199,13 +212,11 @@ private:
                                     double* b);
  void Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, bool is_two_boundaries, mxArray* x0_m);
-  void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, bool is_two_boundaries, mxArray* x0_m,
+  void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, bool is_two_boundaries, mxArray* x0_m);
-                             int precond);
  void Check_and_Correct_Previous_Iteration();
  bool Simulate_One_Boundary();
  bool solve_linear(bool do_check_and_correct);
  void solve_non_linear();
-  string preconditioner_print_out(string s, int preconditioner, bool ss);
  bool compare(int* save_op, int* save_opa, int* save_opaa, int beg_t, long nop4);
  void Insert(int r, int c, int u_index, int lag_index);
  void Delete(int r, int c);
@@ -222,18 +233,6 @@ private:
  int complete(int beg_t);
  void bksub(int tbreak, int last_period);
  void simple_bksub();
-  // Computes Aᵀ where A is are sparse. The result is sparse.
-  static mxArray* Sparse_transpose(const mxArray* A_m);
-  // Computes Aᵀ·B where A and B are sparse. The result is sparse.
-  static mxArray* Sparse_mult_SAT_SB(const mxArray* A_m, const mxArray* B_m);
-  // Computes Aᵀ·B where A is sparse and B is dense. The result is sparse.
-  static mxArray* Sparse_mult_SAT_B(const mxArray* A_m, const mxArray* B_m);
-  // Computes Aᵀ·B where A is sparse and B is dense. The result is dense.
-  static mxArray* mult_SAT_B(const mxArray* A_m, const mxArray* B_m);
-  // Computes A−B where A and B are sparse. The result is sparse.
-  static mxArray* Sparse_subtract_SA_SB(const mxArray* A_m, const mxArray* B_m);
-  // Computes A−B where A and B are dense. The result is dense.
-  static mxArray* subtract_A_B(const mxArray* A_m, const mxArray* B_m);
  void compute_block_time(int my_Per_u_, bool evaluate, bool no_derivatives);
  bool compute_complete(bool no_derivatives);
@@ -248,7 +247,8 @@ public:
              int stack_solve_algo_arg, int solve_algo_arg, bool print_arg,
              const mxArray* GlobalTemporaryTerms_arg, bool steady_state_arg,
              bool block_decomposed_arg, int col_x_arg, int col_y_arg,
-              const BasicSymbolTable& symbol_table_arg, int verbosity_arg);
+              const BasicSymbolTable& symbol_table_arg, int verbosity_arg,
+              iter_solver_opts_t iter_solver_opts_arg);
  pair<bool, vector<int>>
  extended_path(const string& file_name, bool evaluate, int block, int nb_periods,
                const vector<s_plan>& sextended_path,

--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -19,6 +19,7 @@
 #include <algorithm>
 #include <cmath>
+#include <cstring>
 #include <type_traits>
 #include "ErrorHandling.hh"
@@ -480,6 +481,34 @@ mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
      }
  }()};
+  iter_solver_opts_t iter_solver_opts;
+  if (steady_state)
+    iter_solver_opts.set_static_values(options_);
+  else
+    {
+      if (stack_solve_algo == 2 || stack_solve_algo == 3)
+        {
+          // Check that a preconditioner other than 'ilu' has not been requested
+          int field {mxGetFieldNumber(options_, "simul")};
+          if (field < 0)
+            mexErrMsgTxt("simul is not a field of options_");
+          mxArray* simul {mxGetFieldByNumber(options_, 0, field)};
+          field = mxGetFieldNumber(simul, "preconditioner");
+          if (field < 0)
+            mexErrMsgTxt("preconditioner is not a field of options_.simul");
+          mxArray* preconditioner {mxGetFieldByNumber(simul, 0, field)};
+          if (!mxIsChar(preconditioner))
+            mexErrMsgTxt("options_.simul.preconditioner should be a character array");
+          char* preconditioner_str {mxArrayToString(preconditioner)};
+          if (std::strcmp(preconditioner_str, "ilu"))
+            mexErrMsgTxt(
+                "The 'bytecode' option is not compatible with a preconditioner other than 'ilu'");
+          mxFree(preconditioner_str);
+        }
+      iter_solver_opts.set_dynamic_values(options_);
+    }
  field = mxGetFieldNumber(M_, "fname");
  if (field < 0)
    mexErrMsgTxt("fname is not a field of M_");
@@ -549,7 +578,8 @@ mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
                          static_cast<int>(col_x),
                          static_cast<int>(col_y),
                          symbol_table,
-                          verbosity};
+                          verbosity,
+                          iter_solver_opts};
  bool r;
  vector<int> blocks;

--- a/mex/sources/libkorder/dynamic_abstract_class.hh
+++ b/mex/sources/libkorder/dynamic_abstract_class.hh
@@ -47,7 +47,9 @@ public:
      dynamic_g1_sparse_rowval_mx {dynamic_g1_sparse_rowval_mx_arg},
      dynamic_g1_sparse_colval_mx {dynamic_g1_sparse_colval_mx_arg},
      dynamic_g1_sparse_colptr_mx {dynamic_g1_sparse_colptr_mx_arg},
-      dynamic_gN_sparse_indices {move(dynamic_gN_sparse_indices_arg)} {};
+      dynamic_gN_sparse_indices {move(dynamic_gN_sparse_indices_arg)}
+  {
+  }
  virtual ~DynamicModelAC() = default;
  virtual void eval(const Vector& y, const Vector& x, const Vector& params, const Vector& ySteady,
                    Vector& residual, const std::map<int, int>& dynToDynpp,

--- a/mex/sources/libkorder/dynamic_dll.cc
+++ b/mex/sources/libkorder/dynamic_dll.cc
@@ -40,16 +40,11 @@ DynamicModelDLL::DynamicModelDLL(const std::string& modName, int order_arg,
  for (int o {0}; o <= order; o++)
    {
      std::string func_name {"dynamic_"s + (o == 0 ? "resid" : "g"s + std::to_string(o))};
-      std::string mex_filename
+      std::string mex_filename {
-      {
 #if !defined(__CYGWIN32__) && !defined(_WIN32)
-        "./"s +
+          "./"s +
 #endif
-            "+"s + modName + "/+sparse/" + func_name + MEXEXT
+          "+"s + modName + "/+sparse/" + func_name + MEXEXT};
-        // clang-format off
-        // As of Clang 16, the RemoveSemicolon option incorrectly removes the following semicolon
-      };
-      // clang-format on
      mex_handle_t handle {load_mex(mex_filename)};
      if (handle)

--- a/mex/sources/libkorder/kord/first_order.cc
+++ b/mex/sources/libkorder/kord/first_order.cc
@@ -35,8 +35,8 @@ std::mutex FirstOrder::mut;
 lapack_int
 FirstOrder::order_eigs(const double* alphar, const double* alphai, const double* beta)
 {
-  return (*alphar * *alphar + *alphai * *alphai
+  return *alphar * *alphar + *alphai * *alphai
-          < *beta * *beta * qz_criterium_global * qz_criterium_global);
+         < *beta * *beta * qz_criterium_global * qz_criterium_global;
 }
 /* Here we solve the linear approximation. The result are the matrices

--- a/mex/sources/libkorder/sylv/GeneralMatrix.hh
+++ b/mex/sources/libkorder/sylv/GeneralMatrix.hh
@@ -44,7 +44,9 @@ private:
  T& orig;
 public:
-  TransposedMatrix(T& orig_arg) : orig {orig_arg} {};
+  TransposedMatrix(T& orig_arg) : orig {orig_arg}
+  {
+  }
 };
 // Syntactic sugar for representing a transposed matrix

--- a/mex/sources/libkorder/sylv/QuasiTriangular.cc
+++ b/mex/sources/libkorder/sylv/QuasiTriangular.cc
@@ -306,7 +306,7 @@ Diagonal::print() const
 bool
 Diagonal::isZero(double p)
 {
-  return (std::abs(p) < EPS);
+  return std::abs(p) < EPS;
 }
 QuasiTriangular::const_col_iter

--- a/mex/sources/libkorder/sylv/QuasiTriangular.hh
+++ b/mex/sources/libkorder/sylv/QuasiTriangular.hh
@@ -267,7 +267,9 @@ class _column_iter : public _matrix_iter<_TRef, _TPtr>
 public:
  _column_iter(_TPtr base, int ds, bool r, int rw) :
-      _matrix_iter<_TRef, _TPtr>(base, ds, r), row(rw) {};
+      _matrix_iter<_TRef, _TPtr>(base, ds, r), row(rw)
+  {
+  }
  _Self&
  operator++() override
  {
@@ -298,8 +300,9 @@ class _row_iter : public _matrix_iter<_TRef, _TPtr>
  int col;
 public:
-  _row_iter(_TPtr base, int ds, bool r, int cl) :
+  _row_iter(_TPtr base, int ds, bool r, int cl) : _matrix_iter<_TRef, _TPtr>(base, ds, r), col(cl)
-      _matrix_iter<_TRef, _TPtr>(base, ds, r), col(cl) {};
+  {
+  }
  _Self&
  operator++() override
  {

--- a/mex/sources/libkorder/sylv/SchurDecompEig.hh
+++ b/mex/sources/libkorder/sylv/SchurDecompEig.hh
@@ -33,7 +33,9 @@ public:
  SchurDecompEig(const SqSylvMatrix& m) : SchurDecomp(m)
  {
  }
-  SchurDecompEig(const QuasiTriangular& tr) : SchurDecomp(tr) {};
+  SchurDecompEig(const QuasiTriangular& tr) : SchurDecomp(tr)
+  {
+  }
  SchurDecompEig(QuasiTriangular& tr) : SchurDecomp(tr)
  {
  }

--- a/mex/sources/libkorder/sylv/SylvesterSolver.hh
+++ b/mex/sources/libkorder/sylv/SylvesterSolver.hh
@@ -42,7 +42,7 @@ private:
  static bool
  zeroPad(const SchurDecompZero& kdecomp)
  {
-    return ((kdecomp.getZeroCols() * 3 < kdecomp.getDim() * 2) || (kdecomp.getZeroCols() < 10));
+    return (kdecomp.getZeroCols() * 3 < kdecomp.getDim() * 2) || (kdecomp.getZeroCols() < 10);
  }
 public:

--- a/mex/sources/libkorder/sylv/tests/tests.cc
+++ b/mex/sources/libkorder/sylv/tests/tests.cc
@@ -142,7 +142,7 @@ TestRunnable::quasi_solve(bool trans, const std::string& mname, const std::strin
  xx.add(1.0, x);
  double norm = xx.getNorm();
  std::cout << "\terror norm = " << norm << std::endl;
-  return (norm < eps_norm);
+  return norm < eps_norm;
 }
 bool
@@ -174,7 +174,7 @@ TestRunnable::mult_kron(bool trans, const std::string& mname, const std::string&
  c.add(-1.0, v);
  double norm = c.getNorm();
  std::cout << "\terror norm = " << norm << std::endl;
-  return (norm < eps_norm);
+  return norm < eps_norm;
 }
 bool
@@ -208,7 +208,7 @@ TestRunnable::level_kron(bool trans, const std::string& mname, const std::string
  x.add(-1, c);
  double norm = x.getNorm();
  std::cout << "\terror norm = " << norm << std::endl;
-  return (norm < eps_norm);
+  return norm < eps_norm;
 }
 bool
@@ -241,7 +241,7 @@ TestRunnable::kron_power(const std::string& m1name, const std::string& m2name,
  x.add(-1, c);
  double norm = x.getNorm();
  std::cout << "\terror norm = " << norm << std::endl;
-  return (norm < eps_norm);
+  return norm < eps_norm;
 }
 bool
@@ -282,7 +282,7 @@ TestRunnable::lin_eval(const std::string& m1name, const std::string& m2name,
  double norm1 = x1.getNorm();
  double norm2 = x2.getNorm();
  std::cout << "\terror norm1 = " << norm1 << "\n\terror norm2 = " << norm2 << '\n';
-  return (norm1 * norm1 + norm2 * norm2 < eps_norm * eps_norm);
+  return norm1 * norm1 + norm2 * norm2 < eps_norm * eps_norm;
 }
 bool
@@ -323,7 +323,7 @@ TestRunnable::qua_eval(const std::string& m1name, const std::string& m2name,
  double norm1 = x1.getNorm();
  double norm2 = x2.getNorm();
  std::cout << "\terror norm1 = " << norm1 << "\n\terror norm2 = " << norm2 << std::endl;
-  return (norm1 * norm1 + norm2 * norm2 < 100 * eps_norm * eps_norm); // relax norm
+  return norm1 * norm1 + norm2 * norm2 < 100 * eps_norm * eps_norm; // relax norm
 }
 bool
@@ -362,7 +362,7 @@ TestRunnable::tri_sylv(const std::string& m1name, const std::string& m2name,
  double max = dcheck.getMax();
  double xmax = v.getMax();
  std::cout << "\trel. error max = " << max / xmax << std::endl;
-  return (norm < xnorm * eps_norm);
+  return norm < xnorm * eps_norm;
 }
 bool
@@ -388,8 +388,8 @@ TestRunnable::gen_sylv(const std::string& aname, const std::string& bname, const
  gs.check(mmd.getData());
  const SylvParams& pars = gs.getParams();
  pars.print("\t");
-  return (*(pars.mat_err1) < eps_norm && *(pars.mat_errI) < eps_norm && *(pars.mat_errF) < eps_norm
+  return *(pars.mat_err1) < eps_norm && *(pars.mat_errI) < eps_norm && *(pars.mat_errF) < eps_norm
-          && *(pars.vec_err1) < eps_norm && *(pars.vec_errI) < eps_norm);
+         && *(pars.vec_err1) < eps_norm && *(pars.vec_errI) < eps_norm;
 }
 bool
@@ -424,7 +424,7 @@ TestRunnable::eig_bubble(const std::string& aname, int from, int to)
            << "\tabs. error∞ = " << normInf << std::endl
            << "\trel. error1 = " << norm1 / onorm1 << std::endl
            << "\trel. error∞ = " << normInf / onormInf << std::endl;
-  return (norm1 < eps_norm * onorm1 && normInf < eps_norm * onormInf);
+  return norm1 < eps_norm * onorm1 && normInf < eps_norm * onormInf;
 }
 bool
@@ -463,7 +463,7 @@ TestRunnable::block_diag(const std::string& aname, double log10norm)
  std::cout << "\terror Q·Q⁻¹:" << std::endl
            << "\tabs. error1 = " << nor1 << std::endl
            << "\tabs. error∞ = " << norInf << std::endl;
-  return (norm1 < eps_norm * pow(10, log10norm) * onorm1);
+  return norm1 < eps_norm * pow(10, log10norm) * onorm1;
 }
 bool
@@ -503,7 +503,7 @@ TestRunnable::iter_sylv(const std::string& m1name, const std::string& m2name,
  double max = dcheck.getMax();
  double xmax = v.getMax();
  std::cout << "\trel. error max = " << max / xmax << std::endl;
-  return (cnorm < xnorm * eps_norm);
+  return cnorm < xnorm * eps_norm;
 }
 /**********************************************************/
No results found