doc/manual/source/bibliography.rst

@@ -6,6 +6,7 @@ Bibliography
 
 * Abramowitz, Milton and Irene A. Stegun (1964): “Handbook of Mathematical Functions”, Courier Dover Publications.
 * Adjemian, Stéphane, Matthieu Darracq Parriès and Stéphane Moyen (2008): “Towards a monetary policy evaluation framework”, *European Central Bank Working Paper*, 942.
+* Adjemian, Stéphane and Michel Juillard (2025): “Stochastic Extended Path”, *Dynare Working Papers*, 84, CEPREMAP.
 * Aguiar, Mark and Gopinath, Gita (2004): “Emerging Market Business Cycles: The Cycle is the Trend,” *NBER* Working Paper, 10734.
 * Amisano, Gianni and Tristani, Oreste (2010): “Euro area inflation persistence in an estimated nonlinear DSGE model”, *Journal of Economic Dynamics and Control*, 34(10), 1837–1858.
 * Andreasen, Martin M., Jesús Fernández-Villaverde, and Juan Rubio-Ramírez (2018): “The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications,” *Review of Economic Studies*, 85(1), 1-49.

doc/manual/source/installation-and-configuration.rst

@@ -112,6 +112,7 @@ On macOS
    Several versions of Dynare can coexist (by default in ``/Applications/Dynare``),
    as long as you correctly adjust your path settings (see :ref:`words-warning`).
 
+
 With MATLAB
 ^^^^^^^^^^^
 
@@ -125,6 +126,10 @@ The default installation directory is ``/Applications/Dynare/x.y-arch``.
 It is recommended to install the Xcode Command Line Tools (this is an Apple product)
 and GCC via Homebrew_ (see :ref:`prerequisites-macos`).
 
+To deinstall Dynare, simply delete the folder where you installed the program. The package installer does 
+not put any files anywhere else in the system.
+
+
 With Octave
 ^^^^^^^^^^^
 

doc/manual/source/running-dynare.rst

@@ -188,10 +188,13 @@ by the ``dynare`` command.
         Instructs Dynare to no create a logfile of this run in
         ``FILENAME.log.`` The default is to create the logfile.
 
-    .. option:: output=second|third
+    .. option:: output=first|second|third
 
        Instructs the preprocessor to output derivatives of the dynamic model at
-       least up to the given order.
+       least up to the given order. The `first` option is useful in
+       larger models when debugging steady state computation, because it allows
+       overriding the default computation and output of dynamic second order derivatives 
+       in case of the mod-file not containing commands for further computations.
 
     .. option:: language=matlab|julia
 

doc/manual/source/the-model-file.rst

@@ -959,6 +959,15 @@ The model is declared inside a ``model`` block:
         MODEL_EXPRESSION;
+    .. warning::
+
+        In Dynare, only equality signs can delineate the left and right-hand side of an 
+        equation. If Dynare encounters an expression like ``a>=b``, this will therefore not 
+        define an inequality constraint. Rather, it is interpreted as the homogenous equation
+        `(a>=b)=0;`, i.e., the Boolean `(a>=b)` must evaluate to 0. Inequality constraints 
+        in Dynare instead need to be set up either via OccBin or as mixed complementarity problems.
+
+
     |br| Inside the model block, Dynare allows the creation of
     *model-local variables*, which constitute a simple way to share a
     common expression between several equations. The syntax consists
@@ -2981,8 +2990,9 @@ Finding the steady state with Dynare nonlinear solver
            ``5``
                 Newton algorithm with a sparse Gaussian elimination (SPE)
-                solver at each iteration (requires ``bytecode`` option, see
-                :ref:`model-decl`).
+                solver at each iteration. This algorithm requires the
+                :opt:`bytecode` option. The :ref:`markowitz <steady_markowitz>`
+                option can be used to control the behaviour of the algorithm.
            ``6``
@@ -3775,13 +3785,20 @@ 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.
+               of all equations in all periods. The following options can be
+               used to control the behaviour of the algorithm:
+               :opt:`preconditioner <preconditioner = OPTION>`, :opt:`iter_tol
+               <iter_tol = DOUBLE>`, :opt:`iter_maxit <iter_maxit = INTEGER>`,
+               :opt:`gmres_restart <gmres_restart = INTEGER>`.
            ``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.
+               system of all equations in all periods. The following options
+               can be used to control the behaviour of the algorithm:
+               :opt:`preconditioner <preconditioner = OPTION>`, :opt:`iter_tol
+               <iter_tol = DOUBLE>`, :opt:`iter_maxit <iter_maxit = INTEGER>`.
            ``4``
@@ -3795,8 +3812,11 @@ speed-up on large models.
                Use the Laffargue-Boucekkine-Juillard (LBJ) algorithm proposed
                in *Juillard (1996)* on top of a sparse Gaussian elimination
                (SPE) solver. The latter takes advantage of the similarity of
-               the Jacobian across periods when searching for the pivots
-               (requires ``bytecode`` option, see :ref:`model-decl`).
+               the Jacobian across periods when searching for the pivots. This
+               algorithm requires the :opt:`bytecode` option. The following
+               options can be used to control the behaviour of the algorithm:
+               :ref:`markowitz <dynamic_markowitz>`,
+               :opt:`minimal_solving_periods <minimal_solving_periods = INTEGER>`.
            ``6``
@@ -3819,6 +3839,109 @@ speed-up on large models.
                trigger the computation of the solution with a trust
                region algorithm.
+    .. option:: preconditioner = OPTION
+
+       When :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` is equal to
+       ``2`` or ``3``, this option specifies which preconditioner will be used
+       in combination with the iterative sparse linear solver (either GMRES or
+       BiCGStab). Possible values for OPTION are:
+
+           ``umfiter``
+
+               At the first iteration of the nonlinear Newton solver, compute
+               the full LU decomposition with complete pivoting of the linear
+               system (and use it to solve that first iteration rather than
+               using the iterative solver). This LU decomposition is then used
+               as the preconditioner in further Newton iterations. Inspired
+               from TROLL’s option with the same name.
+
+           ``iterstack``
+
+               Compute the LU decomposition with complete pivoting for only a
+               few simulation periods within the stacked Jacobian (which is a
+               block tridiagonal matrix). Then repeat that LU decomposition
+               over the block diagonal to construct a preconditioner for the
+               linear system with all simulation periods. If the total number
+               of simulations periods is not a multiple of the number of
+               periods used for the small LU, then an additional LU is computed
+               for the remainder. The following options can be used to control
+               the construction of this preconditioner:
+               :opt:`iterstack_maxlu <iterstack_maxlu = INTEGER>`,
+               :opt:`iterstack_nperiods <iterstack_nperiods = INTEGER>`,
+               :opt:`iterstack_nlu <iterstack_nlu = INTEGER>`,
+               :opt:`iterstack_relu <iterstack_relu = DOUBLE>`.
+               Inspired from TROLL’s solver with the same name.
+
+           ``ilu``
+
+               Use an incomple LU decomposition as the preconditioner,
+               recomputed at every iteration of the nonlinear Newton solver.
+
+       |br| Default value is ``umfiter``.
+
+    .. option:: iter_tol = DOUBLE
+
+       When :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` is equal to
+       ``2`` or ``3``, this option controls the relative tolerance of the
+       iterative linear solver (either GMRES or BiCGStab). It corresponds to
+       the ``tol`` option of the ``gmres`` and ``bicgstab`` MATLAB/Octave
+       functions. Note that the perfect foresight solver uses an *absolute*
+       tolerance for determining convergence, so this option should be used
+       with care, and the default is meant to suit most situations.
+       Default: the value of the :opt:``tolf <tolf = DOUBLE>`` option, divided
+       by 10 times the infinite norm of the right-hand side of the linear system.
+
+    .. option:: iter_maxit = INTEGER
+
+       When :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` is equal to
+       ``2`` or ``3``, this option controls the maximum number of iterations of
+       the iterative linear solver (either GMRES or BiCGStab). It corresponds
+       to the ``maxit`` option of the ``gmres`` and ``bicgstab`` MATLAB/Octave
+       functions. It should not be confused with the :opt:`maxit <maxit = INTEGER>`
+       option, which controls the (outer) nonlinear Newton loop, while
+       ``iter_maxit`` controls the (inner) linear loop. Default: ``500``.
+
+    .. option:: gmres_restart = INTEGER
+
+       When :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` is equal to
+       ``2``, this option controls the number of iterations before restart of
+       the GMRES algorithm. It corresponds to the ``restart`` option of the
+       ``gmres`` MATLAB/Octave function. Default: ``100``.
+
+    .. option:: iterstack_maxlu = INTEGER
+
+       When :opt:`preconditioner <preconditioner = OPTION>` is equal to
+       ``iterstack``, controls the maximum size of the matrix for which the
+       small LU will computed. The actual size of the matrix will be determined
+       by the largest number of periods that, multiplied by the number of
+       equations, is less than the value of the option. Note that when combined
+       with block decomposition (see :opt:`block`), blocks for which the
+       whole stacked system is less than this option will be solved using a
+       regular LU decomposition instead of the iterative linear solver.
+       Default: ``20000``.
+
+    .. option:: iterstack_nperiods = INTEGER
+
+       When :opt:`preconditioner <preconditioner = OPTION>` is equal to
+       ``iterstack``, controls the number of periods used for the small LU.
+       If nonzero, this option overrides the :opt:`iterstack_maxlu <iterstack_maxlu = INTEGER>`
+       and :opt:`iterstack_nlu <iterstack_nlu = INTEGER>` options. Default: ``0``.
+
+    .. option:: iterstack_nlu = INTEGER
+
+       When :opt:`preconditioner <preconditioner = OPTION>` is equal to
+       ``iterstack``, specifies the number of times that the small LU should be
+       repeated in the large preconditioner. If nonzero, this option overrides
+       the :opt:`iterstack_maxlu <iterstack_maxlu = INTEGER>` option. Default:
+       ``0``.
+
+    .. option:: iterstack_relu = DOUBLE
+
+       When :opt:`preconditioner <preconditioner = OPTION>` is equal to
+       ``iterstack``, controls the relative position of the small LU within the
+       whole stacked system. Must be a number between ``0`` and ``1``. Default:
+       ``0.5``.
+
     .. option:: robust_lin_solve
        Triggers the use of a robust linear solver for the default
@@ -3919,11 +4042,13 @@ speed-up on large models.
        procedure, *i.e.* which must be kept at their value corresponding to
        100% of the shock during all homotopy iterations.
+    .. _dynamic_markowitz:
+
     .. option:: markowitz = DOUBLE
-       Value of the Markowitz criterion, used to select the
-       pivot. Only used when ``stack_solve_algo = 5``. Default:
-       ``0.5``.
+       Value of the Markowitz criterion, used to select the pivot (see
+       :ref:`markowitz <steady_markowitz>` for more details). Only used when
+       ``stack_solve_algo = 5``. Default: ``0.5``.
     .. option:: minimal_solving_periods = INTEGER
@@ -5478,7 +5603,8 @@ which is described below.
        If order is greater than ``0`` Dynare uses a gaussian
        quadrature to take into account the effects of future
-       uncertainty. If ``order`` :math:`=S` then the time series for
+       uncertainty; this is called *stochastic* extended path, see *Adjemian
+       and Juillard (2025)*. If ``order`` :math:`=S` then the time series for
        the endogenous variables are generated by assuming that the
        agents believe that there will no more shocks after period
        :math:`t+S`. This is an experimental feature and can be quite

matlab/default_option_values.m

@@ -58,6 +58,9 @@ options_.dr_display_tol=1e-6;
 options_.dp.maxit = 3000;
 options_.steady.maxit = 50;
 options_.steady.non_zero = false;
+options_.steady.ilu.type = 'ilutp';
+options_.steady.ilu.droptol = 1e-10;
+options_.steady.ilu.udiag = true;
 options_.simul.maxit = 50;
 options_.simul.robust_lin_solve = false;
 
@@ -328,6 +331,19 @@ options_.simul.endval_steady = false;
 options_.simul.first_simulation_period = dates();
 options_.simul.last_simulation_period = dates();
 
+% For stack_solve_algo={2,3}
+options_.simul.preconditioner = 'umfiter';
+options_.simul.iter_tol = []; % See perfect-foresight-models/iter_solver_params.m for details
+options_.simul.iter_maxit = 500;
+options_.simul.gmres_restart = 100;
+options_.simul.iterstack_maxlu = 20000;
+options_.simul.iterstack_nperiods = 0;
+options_.simul.iterstack_nlu = 0;
+options_.simul.iterstack_relu = 0.5;
+options_.simul.ilu.type = 'ilutp';
+options_.simul.ilu.droptol = 1e-12;
+options_.simul.ilu.udiag = true;
+
 options_.simul.homotopy_max_completion_share = 1;
 options_.simul.homotopy_min_step_size = 1e-3;
 options_.simul.homotopy_step_size_increase_success_count = 3;

matlab/getPowerDeriv.m

@@ -30,7 +30,7 @@ function dxp=getPowerDeriv(x,p,k)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if (abs(x) < 1e-12) && (p > 0) && (k > p) && (abs(p - round(p)) < 1e-12)
+if (abs(x) < 1e-12) && (p >= 0) && (k > p) && (abs(p - round(p)) < 1e-12)
     dxp = 0;
 else
     dxp = x^(p-k);
@@ -39,4 +39,23 @@ else
         p = p-1;
     end
 end
-end
+
+return % --*-- Unit tests --*--
+
+%@test:1
+x=getPowerDeriv(2,3,1);
+t(1)=all(abs(x-3*4)<1e-10);
+x=getPowerDeriv(0,2,2);
+t(2)=all(abs(x-2)<1e-10);
+x=getPowerDeriv(0,2,3); %special case evaluates to 0
+t(3)=all(abs(x-0)<1e-10);
+x=getPowerDeriv(1e-13,2,3-1e-13); %0 within tolerance
+t(4)=all(abs(x-0)<1e-10);
+x=getPowerDeriv(0,0,1);
+t(5)=all(abs(x-0)<1e-10);
+x=getPowerDeriv(0,0,0);
+t(6)=all(abs(x-1)<1e-10);
+x=getPowerDeriv(0,1/3,1); %derivative evaluating to Inf due to division by 0
+t(7)= isinf(x);
+T = all(t);
+%@eof:1
\ No newline at end of file

matlab/kalman/likelihood/missing_observations_kalman_filter.m

@@ -72,13 +72,14 @@ if isequal(H,0)
     H = zeros(pp,pp);
 end
 
+P=tril(P)+transpose(tril(P,-1)); % make sure P is symmetric
+
 % Get sample size.
 smpl = last-start+1;
 
 % Initialize some variables.
-dF   = 1;
 isqvec = false;
-if ndims(Q)>2
+if ~ismatrix(Q)
     Qvec = Q;
     Q=Q(:,:,1);
     isqvec = true;

matlab/kalman/likelihood/missing_observations_kalman_filter_d.m

@@ -64,9 +64,8 @@ function [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(data_index,n
 smpl = last-start+1;
 
 % Initialize some variables.
-dF   = 1;
 isqvec = false;
-if ndims(Q)>2
+if ~ismatrix(Q)
     Qvec = Q;
     Q=Q(:,:,1);
     isqvec = true;
@@ -75,7 +74,6 @@ QQ   = R*Q*transpose(R);   % Variance of R times the vector of structural innova
 t    = start;              % Initialization of the time index.
 dlik = zeros(smpl,1);      % Initialization of the vector gathering the densities.
 dLIK = Inf;                % Default value of the log likelihood.
-oldK = Inf;
 
 if isequal(H,0)
     H = zeros(pp,pp);

matlab/optimization/dynare_solve.m

@@ -22,7 +22,7 @@ function [x, errorflag, fvec, fjac, errorcode] = dynare_solve(f, x, maxit, tolf,
 %        -10  -> System of equation ill-behaved at the initial guess (Inf, Nans or complex numbers).
 %        -11  -> Initial guess is a solution of the system of equations.
 
-% Copyright © 2001-2024 Dynare Team
+% Copyright © 2001-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -294,7 +294,7 @@ elseif options_.solve_algo==3
     end
     [fvec, fjac] = feval(f, x, varargin{:});
 elseif ismember(options_.solve_algo, [6 7 8])
-    [x, errorflag, errorcode] = newton_solve(f, x, jacobian_flag, options_.gstep, tolf, tolx, maxit, options_.solve_algo, varargin{:});
+    [x, errorflag, errorcode] = newton_solve(f, x, jacobian_flag, options_.gstep, tolf, tolx, maxit, options_.solve_algo, options_.steady.ilu, varargin{:});
     [fvec, fjac] = feval(f, x, varargin{:});
 elseif options_.solve_algo==10
     % LMMCP

matlab/optimization/newton_solve.m

-function [x, errorflag, errorcode] = newton_solve(func, x, jacobian_flag, gstep, tolf, tolx, maxit, solve_algo, varargin)
+function [x, errorflag, errorcode] = newton_solve(func, x, jacobian_flag, gstep, tolf, tolx, maxit, solve_algo, ilu_opts, varargin)
 
 % Solves systems of non linear equations of several variables using a Newton solver, with three
 % variants for the inner linear solver:
@@ -17,13 +17,14 @@ function [x, errorflag, errorcode] = newton_solve(func, x, jacobian_flag, gstep,
 %    tolx             tolerance for solution variation
 %    maxit            maximum number of iterations
 %    solve_algo       from options_
+%    ilu_opts         structure of options for ilu
 %    varargin:        list of extra arguments to the function
 %
 % OUTPUTS
 %    x:               results
 %    errorflag=true:  the model can not be solved
 
-% Copyright © 2024 Dynare Team
+% Copyright © 2024-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -103,9 +104,8 @@ for it = 1:maxit
     if solve_algo == 6
         p = fjac\fvec;
     else
-        ilu_setup.type = 'ilutp';
-        ilu_setup.droptol = 1e-10;
-        [L1, U1] = ilu(fjac, ilu_setup);
+        % Should be the same options as in solve_one_boundary.m and bytecode/Interpreter.cc (static case)
+        [L1, U1] = ilu(fjac, ilu_opts);
         if solve_algo == 7
             p = gmres(fjac, fvec, [], [], [], L1, U1);
         elseif solve_algo == 8

matlab/optimization/solve_one_boundary.m

@@ -40,7 +40,7 @@ function [y, T, success, max_res, iter] = solve_one_boundary(fh, y, x, params, s
 % ALGORITHM
 %   Newton with LU or GMRES or BicGstab for dynamic block
 
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -59,11 +59,8 @@ function [y, T, success, max_res, iter] = solve_one_boundary(fh, y, x, params, s
 
 Blck_size=size(y_index_eq,2);
 correcting_factor=0.01;
-ilu_setup.type='crout';
-ilu_setup.droptol=1e-10;
 max_resa=1e100;
 lambda = 1; % Length of Newton step
-reduced = 0;
 if is_forward
     incr = 1;
     start = y_kmin+1;
@@ -141,7 +138,6 @@ for it_=start:incr:finish
                         end
                     elseif lambda>1e-8
                         lambda=lambda/2;
-                        reduced = 1;
                         if verbose
                             disp(['reducing the path length: lambda=' num2str(lambda,'%f')])
                         end
@@ -192,7 +188,9 @@ for it_=start:incr:finish
                                     M_.block_structure.block(Block_Num).g1_sparse_rowval, ...
                                     M_.block_structure.block(Block_Num).g1_sparse_colval, ...
                                     M_.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
-            elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0 || stack_solve_algo==6)) || (~is_dynamic && options_.solve_algo==6)
+            elseif (is_dynamic && (ismember(stack_solve_algo, [0 1 6]) || ...
+                                   (ismember(stack_solve_algo, [2 3]) && strcmp(options_.simul.preconditioner, 'iterstack')))) ... % Iterstack does not make sense for one boundary problems
+                   || (~is_dynamic && options_.solve_algo==6)
                 if verbose && ~is_dynamic
                     disp('steady: Sparse LU ')
                 end
@@ -203,64 +201,43 @@ for it_=start:incr:finish
                 else
                     y(y_index_eq) = ya;
                 end
-            elseif (stack_solve_algo==2 && is_dynamic) || (options_.solve_algo==7 && ~is_dynamic)
-                flag1=1;
-                if verbose && ~is_dynamic
-                    disp('steady: GMRES ')
+            elseif is_dynamic && ismember(stack_solve_algo, [2 3])
+                if strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0
+                    [L, U, P, Q] = lu(g1);
+                    zc = L\(P*r);
+                    zb = U\zc;
+                elseif strcmp(options_.simul.preconditioner, 'ilu')
+                    [L, U, P] = ilu(g1, options_.simul.ilu);
+                    Q = speye(size(g1));
                 end
-                while flag1>0
-                    [L1, U1]=ilu(g1,ilu_setup);
-                    [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
-                    if  flag1>0 || reduced
-                        if verbose
-                            if flag1==1
-                                disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')])
-                            elseif(flag1==2)
-                                disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')])
-                            elseif(flag1==3)
-                                disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')])
-                            end
-                        end
-                        ilu_setup.droptol = ilu_setup.droptol/10;
-                        reduced = 0;
+                if ~(strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0)
+                    [iter_tol, iter_maxit, gmres_restart] = iter_solver_params(options_, g1, r);
+                    if stack_solve_algo==2
+                        [zb, flag] = gmres(P*g1*Q, P*r, gmres_restart, iter_tol, iter_maxit, L, U);
                     else
-                        ya = ya + lambda*dx;
-                        if is_dynamic
-                            y(y_index_eq, it_) = ya;
-                        else
-                            y(y_index_eq) = ya';
-                        end
+                        [zb, flag] = bicgstab(P*g1*Q, P*r, iter_tol, iter_maxit, L, U);
                     end
+                    iter_solver_error_flag(flag)
                 end
-            elseif (stack_solve_algo==3 && is_dynamic) || (options_.solve_algo==8 && ~is_dynamic)
-                flag1=1;
-                if verbose && ~is_dynamic
-                    disp('steady: BiCGStab')
-                end
-                while flag1>0
-                    [L1, U1]=ilu(g1,ilu_setup);
-                    [dx,flag1] = bicgstab(g1,-r,1e-6,Blck_size,L1,U1);
-                    if flag1>0 || reduced
-                        if verbose
-                            if(flag1==1)
-                                disp(['Error in simul: No convergence inside BiCGStab after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')])
-                            elseif(flag1==2)
-                                disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')])
-                            elseif(flag1==3)
-                                disp(['Error in simul: BiCGStab stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')])
-                            end
-                        end
-                        ilu_setup.droptol = ilu_setup.droptol/10;
-                        reduced = 0;
+                ya = ya - lambda*Q*zb;
+                y(y_index_eq, it_) = ya;
+            elseif ~is_dynamic && ismember(options_.solve_algo, [7 8])
+                if verbose
+                    if options_.solve_algo == 7
+                        disp('steady: GMRES ')
                     else
-                        ya = ya + lambda*dx;
-                        if is_dynamic
-                            y(y_index_eq, it_) = ya;
-                        else
-                            y(y_index_eq) = ya';
-                        end
+                        disp('steady: BiCGStab')
                     end
                 end
+                %% Should be the same options as in newton_solve.m and bytecode/Interpreter.cc (static case)
+                [L, U] = ilu(g1, options_.steady.ilu);
+                if options_.solve_algo == 7
+                    zb = gmres(g1, r, [], [], [], L, U);
+                else
+                    zb = bicgstab(g1, r, [], [], L, U);
+                end
+                ya = ya - lambda*zb;
+                y(y_index_eq) = ya;
             else
                 if is_dynamic
                     error(['options_.stack_solve_algo = ' num2str(stack_solve_algo) ' not implemented'])

matlab/perfect-foresight-models/iter_solver_error_flag.m

+function iter_solver_error_flag(flag)
+% Interprets error flag from gmres or bicgstab functions
+
+% Copyright © 2025 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+switch flag
+    case 1
+        error('Maximum number of iterations exceeded in GMRES/BiCGStab. You may want to increase the iter_maxit option or the gmres_restart option (if using stack_solve_algo=2 for the latter), or decrease iter_tol.')
+    case 2
+        error('The preconditioner matrix is ill conditioned in GMRES/BiCGStab. You should try another preconditioner.')
+    case 3
+        error('No progress between two successive iterations in GMRES/BiCGStab. You may want to increase the iter_tol option.')
+    case 4
+        error('One of the scalar quantities calculated by GMRES/BiCGStab became too small or too large. Try another preconditioner or algorithm.')
+end

matlab/perfect-foresight-models/iter_solver_params.m

+function [iter_tol, iter_maxit, gmres_restart] = iter_solver_params(options_, A, b)
+% Computes tolerance used for gmres or bicgstab functions in perfect foresight solver.
+% b is the RHS of the linear equation to solve. Also ensures that maxit and restart
+% are within acceptable range.
+%
+% The tolerance passed to gmres and bicgstab is a relative one (‖Ax−b‖/‖b‖).
+% However, we test the convergence of algorithms via options_.dynatol.f, which is an
+% absolute error (‖Ax−b‖). Hence the need to rescale by the norm of the RHS (‖b‖).
+
+% Copyright © 2025 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+iter_tol = options_.simul.iter_tol;
+if isempty(iter_tol)
+    iter_tol = options_.dynatol.f / max(abs(b)) / 10;
+end
+
+iter_maxit = min(options_.simul.iter_maxit, size(A, 1));
+gmres_restart = min(options_.simul.gmres_restart, size(A, 1));

matlab/perfect-foresight-models/iterstack_preconditioner.m

+function [L, U, P, Q] = iterstack_preconditioner(A, options_)
+% Copyright © 2025 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+ny = size(A, 1) / options_.periods;
+
+if options_.simul.iterstack_nperiods > 0
+    if options_.simul.iterstack_nperiods > options_.periods
+        error('iterstack preconditionner: iterstack_nperiods option is greater than the periods options')
+    end
+    lusize = ny * options_.simul.iterstack_nperiods;
+elseif options_.simul.iterstack_nlu ~= 0
+    if options_.simul.iterstack_nlu < 0 % To avoid misleading TROLL users
+        error('iterstack preconditionner: negative value of iterstack_nlu option is not supported')
+    end
+    lusize = floor(floor(size(A, 1) / options_.simul.iterstack_nlu) / ny) * ny;
+    if lusize == 0
+        error('iterstack preconditionner: iterstack_nlu option is too large')
+    end
+else
+    if size(A, 1) < options_.simul.iterstack_maxlu
+        error('iterstack preconditionner: too small problem; either decrease iterstack_maxlu option, or use iterstack_nperiods or iterstack_nlu options')
+    end
+    if ny > options_.simul.iterstack_maxlu
+        error('iterstack preconditionner: too large problem; either increase iterstack_maxlu option, or use iterstack_nperiods or iterstack_nlu options')
+    end
+    lusize = floor(options_.simul.iterstack_maxlu / ny) * ny;
+end
+
+if options_.simul.iterstack_relu < 0 || options_.simul.iterstack_relu > 1
+    error('iterstack preconditionner: iterstack_relu option must be between 0 and 1')
+end
+
+luidx = floor((floor(size(A, 1) / lusize) - 1) * options_.simul.iterstack_relu) * lusize + (1:lusize);
+[L1, U1, P1, Q1] = lu(A(luidx,luidx));
+eyek = speye(floor(size(A, 1) / lusize));
+L = kron(eyek, L1);
+U = kron(eyek, U1);
+P = kron(eyek, P1);
+Q = kron(eyek, Q1);
+r = rem(size(A, 1), lusize);
+if r > 0 % Compute additional smaller LU for remainder, if any
+    [L2, U2, P2, Q2] = lu(A(end-r+1:end,end-r+1:end));
+    L = blkdiag(L, L2);
+    U = blkdiag(U, U2);
+    P = blkdiag(P, P2);
+    Q = blkdiag(Q, Q2);
+end
+
+if options_.debug
+    fprintf('iterstack preconditionner: main LU size = %d, remainder LU size = %d\n', lusize, r)
+end

matlab/perfect-foresight-models/perfect_foresight_solver.m

@@ -215,7 +215,7 @@ elseif options_.simul.homotopy_marginal_linearization_fallback > 0 && completed_
             if is_numerical_exception(ME)
                 extra_success = false;
             else
-                rethrow ME
+                rethrow(ME)
             end
         end
     end
@@ -375,7 +375,7 @@ while step > options_.simul.homotopy_min_step_size
                 solver_iter = [];
                 per_block_status = [];
             else
-                rethrow ME
+                rethrow(ME)
             end
         end
     else
@@ -646,4 +646,4 @@ function r = is_numerical_exception(ME)
 r = (~isoctave && (strcmp(ME.identifier, 'MATLAB:erf:notFullReal') ...
                    || strcmp(ME.identifier, 'MATLAB:erfc:notFullReal'))) ...
     || (isoctave && (strcmp(ME.message, 'normcdf: X, MU, and SIGMA must not be complex.') ...
-                     || strcmp(ME.message, 'normpdf: X, MU, and SIGMA must not be complex.')))
+                     || strcmp(ME.message, 'normpdf: X, MU, and SIGMA must not be complex.')));

matlab/perfect-foresight-models/sim1.m

@@ -82,6 +82,7 @@ end
 h1 = clock;
 iter = 1;
 converged = false;
+umfiter_precond = [];
 
 while ~(converged || iter > options_.simul.maxit)
     h2 = clock;
@@ -148,7 +149,8 @@ while ~(converged || iter > options_.simul.maxit)
         if options_.simul.robust_lin_solve
             dy = -lin_solve_robust(A, res, verbose, options_);
         else
-            dy = -lin_solve(A, res, verbose, options_);
+            [mdy, umfiter_precond] = lin_solve(A, res, verbose, options_, umfiter_precond);
+            dy = -mdy;
         end
         if any(isnan(dy)) || any(isinf(dy))
             if verbose
@@ -217,7 +219,8 @@ if verbose
     skipline();
 end
 
-function x = lin_solve(A, b, verbose, options_)
+function [x, umfiter_precond] = lin_solve(A, b, verbose, options_, umfiter_precond)
+
 if norm(b) < sqrt(eps) % then x = 0 is a solution
     x = 0;
     return
@@ -226,16 +229,39 @@ end
 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);
+    if strcmp(options_.simul.preconditioner, 'umfiter')
+        if isempty(umfiter_precond)
+            [L, U, P, Q] = lu(A);
+        else
+            L = umfiter_precond.L;
+            U = umfiter_precond.U;
+            P = umfiter_precond.P;
+            Q = umfiter_precond.Q;
+        end
+    elseif strcmp(options_.simul.preconditioner, 'iterstack')
+        [L, U, P, Q] = iterstack_preconditioner(A, options_);
+    elseif strcmp(options_.simul.preconditioner, 'ilu')
+        [L, U, P] = ilu(A, options_.simul.ilu);
+        Q = speye(size(A));
+    end
+
+    if strcmp(options_.simul.preconditioner, 'umfiter') && isempty(umfiter_precond)
+        z = L\(P*b);
+        y = U\z;
+        umfiter_precond = struct('L', L, 'U', U, 'P', P, 'Q', Q);
     else
-        error('sim1: invalid value for options_.stack_solve_algo')
+        [iter_tol, iter_maxit, gmres_restart] = iter_solver_params(options_, A, b);
+        if options_.stack_solve_algo == 2
+            [y, flag] = gmres(P*A*Q, P*b, gmres_restart, iter_tol, iter_maxit, L, U);
+        elseif options_.stack_solve_algo == 3
+            [y, flag] = bicgstab(P*A*Q, P*b, iter_tol, iter_maxit, L, U);
+        else
+            error('sim1: invalid value for options_.stack_solve_algo')
+        end
+        iter_solver_error_flag(flag)
     end
+
+    x = Q*y;
 end
 x(~isfinite(x)) = 0;
 relres = norm(b - A*x) / norm(b);

matlab/perfect-foresight-models/solve_block_decomposed_problem.m

@@ -14,7 +14,7 @@ function [y, success, maxerror, per_block_status] = solve_block_decomposed_probl
 %   maxerror         [double]    ∞-norm of the residual
 %   per_block_status [struct]    vector structure with per-block information about convergence
 
-% Copyright © 2020-2024 Dynare Team
+% Copyright © 2020-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -40,9 +40,9 @@ switch options_.stack_solve_algo
     case {1,6}
         mthd='LBJ with LU solver';
     case 2
-        mthd='GMRES on stacked system';
+        mthd=sprintf('GMRES on stacked system with ''%s'' preconditioner', options_.simul.preconditioner);
     case 3
-        mthd='BiCGStab on stacked system';
+        mthd=sprintf('BiCGStab on stacked system with ''%s'' preconditioner', options_.simul.preconditioner);
     case 4
         mthd='Sparse LU solver with optimal path length on stacked system';
     case 7

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.
 %
@@ -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
-            flag1=1;
-            while flag1>0
-                if preconditioner==2
-                    [L1, U1]=ilu(g1a,ilu_setup);
-                elseif preconditioner==3
-                    Size = Blck_size;
-                    gss1 =  g1a(Size + 1: 2*Size,Size + 1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
-                    [L1, U1]=lu(gss1);
-                    L(1:Size,1:Size) = L1;
-                    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
+        elseif ismember(stack_solve_algo, [2, 3])
+            if strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0
+                [L, U, P, Q] = lu(g1a);
+                zc = L\(P*b);
+                zb = U\zc;
+            elseif strcmp(options_.simul.preconditioner, 'iterstack')
+                [L, U, P, Q] = iterstack_preconditioner(g1a, options_);
+            elseif strcmp(options_.simul.preconditioner, 'ilu')
+                [L, U, P] = ilu(g1a, options_.simul.ilu);
+                Q = speye(size(g1a));
             end
-        elseif stack_solve_algo==3
-            flag1=1;
-            while flag1>0
-                if preconditioner==2
-                    [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);
+            if ~(strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0)
+                [iter_tol, iter_maxit, gmres_restart] = iter_solver_params(options_, g1a, b);
+                if stack_solve_algo == 2
+                    [zb, flag] = gmres(P*g1a*Q, P*b, gmres_restart, iter_tol, iter_maxit, L, U);
                 else
-                    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);
-                    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);
+                    [zb, flag] = bicgstab(P*g1a*Q, P*b, iter_tol, iter_maxit, L, U);
                 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)]);

meson.build

@@ -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@' ],

mex/sources/bytecode/Interpreter.cc

@@ -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;
@@ -2260,65 +2309,6 @@ Interpreter::simple_bksub()
     }
 }
 
-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_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)
 {
@@ -2750,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),
-                  mxCreateDoubleScalar(1e-6),
-                  mxCreateDoubleScalar(static_cast<double>(n)),
+                  iter_solver_opts.gmres_restart,
+                  iter_solver_opts.iter_tol,
+                  iter_solver_opts.iter_maxit,
                   L1,
                   U1,
                   x0_m};
@@ -2776,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 "
-                       + to_string(block_num + 1))
-                          .c_str());
+        mexErrMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
+                      + to_string(block_num + 1))
+                         .c_str());
       else if (flag1 == 3)
-        mexWarnMsgTxt(("Error in bytecode: GMRES stagnated (Two consecutive iterates were the "
-                       "same.), in block "
-                       + to_string(block_num + 1))
-                          .c_str());
-      lu_inc_tol /= 10;
+        mexErrMsgTxt(("Error in bytecode: GMRES stagnated (Two consecutive iterates were the "
+                      "same.), in block "
+                      + to_string(block_num + 1))
+                         .c_str());
     }
   else
     {
@@ -2819,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);
@@ -2828,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*, 1> lhs0;
-      std::array rhs0 {A_m, mxCreateDoubleScalar(0)};
-      mexCallMATLAB(lhs0.size(), lhs0.data(), rhs0.size(), rhs0.data(), "spdiags");
-      mxArray* tmp = lhs0[0];
-      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);
-    }
+  std::array<mxArray*, 2> lhs0;
+  std::array rhs0 {A_m, iter_solver_opts.ilu};
+  if (mexCallMATLAB(lhs0.size(), lhs0.data(), rhs0.size(), rhs0.data(), "ilu"))
+    throw FatalException {"In BiCGStab, the incomplete LU decomposition (ilu) has failed"};
+  mxArray* L1 = lhs0[0];
+  mxArray* U1 = lhs0[1];
 
-  if (flags == 2)
-    {
-      if (preconditioner == 0)
-        {
-          /*[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)), Diag};
-          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]);
-        }
-      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]);
-        }
-    }
+  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;
+  mexCallMATLAB(lhs.size(), lhs.data(), rhs.size(), rhs.data(), "bicgstab");
+  mxArray* z = lhs[0];
+  mxArray* flag = lhs[1];
+  double flags {mxGetScalar(flag)};
+  mxDestroyArray(flag);
+  mxDestroyArray(rhs[4]);
+  mxDestroyArray(rhs[5]);
 
   if (flags > 0)
     {
       if (flags == 1)
-        mexWarnMsgTxt(("Error in bytecode: No convergence inside BiCGStab, in block "
-                       + to_string(block_num + 1))
-                          .c_str());
+        mexErrMsgTxt(("Error in bytecode: No convergence inside BiCGStab, in block "
+                      + to_string(block_num + 1))
+                         .c_str());
       else if (flags == 2)
-        mexWarnMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
-                       + to_string(block_num + 1))
-                          .c_str());
+        mexErrMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
+                      + to_string(block_num + 1))
+                         .c_str());
       else if (flags == 3)
-        mexWarnMsgTxt(("Error in bytecode: BiCGStab stagnated (Two consecutive iterates were the "
-                       "same.), in block "
-                       + to_string(block_num + 1))
-                          .c_str());
-      lu_inc_tol /= 10;
+        mexErrMsgTxt(("Error in bytecode: BiCGStab stagnated (Two consecutive iterates were the "
+                      "same.), in block "
+                      + to_string(block_num + 1))
+                         .c_str());
     }
   else
     {
@@ -2986,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);
@@ -3941,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();
@@ -4010,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",
-                                                 preconditioner, true)
-                            .c_str());
+              mexPrintf("MODEL STEADY STATE: (method=GMRES with 'ilu' preconditioner)\n");
               break;
             case 8:
-              mexPrintf(preconditioner_print_out("MODEL STEADY STATE: (method=BiCGStab)\n",
-                                                 preconditioner, true)
-                            .c_str());
+              mexPrintf("MODEL STEADY STATE: (method=BiCGStab with 'ilu' preconditioner)\n");
               break;
             }
         }
@@ -4087,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)
@@ -4204,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;
@@ -4394,15 +4240,11 @@ Interpreter::Simulate_Newton_Two_Boundaries(
               break;
             case 2:
               mexPrintf(
-                  preconditioner_print_out("MODEL SIMULATION: (method=GMRES on stacked system)\n",
-                                           preconditioner, false)
-                      .c_str());
+                  "MODEL SIMULATION: (method=GMRES on stacked system with 'ilu' preconditioner)\n");
               break;
             case 3:
-              mexPrintf(preconditioner_print_out(
-                            "MODEL SIMULATION: (method=BiCGStab on stacked system)\n",
-                            preconditioner, false)
-                            .c_str());
+              mexPrintf("MODEL SIMULATION: (method=BiCGStab on stacked system with 'ilu' "
+                        "preconditioner)\n");
               break;
             case 4:
               mexPrintf("MODEL SIMULATION: (method=Sparse LU solver with optimal path length on "
@@ -4460,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);
     }