diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 183c916a3f3e7c33db25fd9180e83ac8a58cb904..cd36b3ecb6e5ea501ebfb5f441b8b3a3ce6b8cbb 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -3785,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``
@@ -3832,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
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 7147d8de24bf223ff92c2b0504e6e7a0dab25800..febf0223e76b6026637a954bebdaac99a823bd1c 100644
--- a/matlab/default_option_values.m
+++ b/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;
diff --git a/matlab/optimization/dynare_solve.m b/matlab/optimization/dynare_solve.m
index 3779aae7c091a15a2c0e63895bc4e89d8b377700..ebd305e7728b4f5591b8fd66c2fb911c9d81795d 100644
--- a/matlab/optimization/dynare_solve.m
+++ b/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
diff --git a/matlab/optimization/newton_solve.m b/matlab/optimization/newton_solve.m
index de1d2e48286ebd95973d632d1d7d3873e7f9a3bb..cdfc953a73e7543cc66952b8ff561498f962d60c 100644
--- a/matlab/optimization/newton_solve.m
+++ b/matlab/optimization/newton_solve.m
@@ -1,4 +1,4 @@
-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
diff --git a/matlab/optimization/solve_one_boundary.m b/matlab/optimization/solve_one_boundary.m
index e9177db7d02acb3d508cc5eb9a7ef8d566b17c3d..f73f29b6eaa76dac6b8dbceea2b6ce45225ee50b 100644
--- a/matlab/optimization/solve_one_boundary.m
+++ b/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'])
diff --git a/matlab/perfect-foresight-models/iter_solver_error_flag.m b/matlab/perfect-foresight-models/iter_solver_error_flag.m
new file mode 100644
index 0000000000000000000000000000000000000000..f7e61dbb791987e030238784f8fe05c6bbac381a
--- /dev/null
+++ b/matlab/perfect-foresight-models/iter_solver_error_flag.m
@@ -0,0 +1,30 @@
+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
diff --git a/matlab/perfect-foresight-models/iter_solver_params.m b/matlab/perfect-foresight-models/iter_solver_params.m
new file mode 100644
index 0000000000000000000000000000000000000000..8f1dc94c480a53b9d11142a4ef341218899fb245
--- /dev/null
+++ b/matlab/perfect-foresight-models/iter_solver_params.m
@@ -0,0 +1,33 @@
+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));
diff --git a/matlab/perfect-foresight-models/iterstack_preconditioner.m b/matlab/perfect-foresight-models/iterstack_preconditioner.m
new file mode 100644
index 0000000000000000000000000000000000000000..8c69e9de476dcf53b58581bb323208934b792f96
--- /dev/null
+++ b/matlab/perfect-foresight-models/iterstack_preconditioner.m
@@ -0,0 +1,66 @@
+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
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 974cda2d5ccdabf650ba0a9384e0b187b6383db6..bc3c70131278c381b1d9ac606e4f35cf2ab7fe5f 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/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);
diff --git a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
index 14db6dee693ee24108b535248d5c5b5147395361..cbfd37b82a3fef5395f3a021cb2e416270108c80 100644
--- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
+++ 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
diff --git a/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m b/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
index e3934e2579cefe5092d4b7f6f0bf59559cf56302..2c09bb1e4644b9726a2bbb7ce5cc616ec3b6cae6 100644
--- 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.
%
@@ -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)]);
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index c7d318d94c3fdb7ef911a75c3ae1feff9f1a68ed..1b361f0ffcb52e6d90d9e69a0a9dcf83823a8217 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/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);
}
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index 0259d5e8e0b3e95a8ce427a8ef615b8e95431948..393bbaeaba63aed3ff168b3f2b13ae932c988b2b 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -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,
- int precond);
+ void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, bool is_two_boundaries, mxArray* x0_m);
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,10 +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 is sparse and B is dense. The result is dense.
- static mxArray* mult_SAT_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);
@@ -240,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,
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index 86bbda4d57c59b6231b295f178c2f98e2e0950a1..1733afe30028235993862af7e60968145a4cf756 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -1,5 +1,5 @@
/*
- * 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;
diff --git a/preprocessor b/preprocessor
index 5c75643443ab51e49e18f7bc6546e85cd3f55536..a7d747aedac558d2d70067a284439e0130e40fa9 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 5c75643443ab51e49e18f7bc6546e85cd3f55536
+Subproject commit a7d747aedac558d2d70067a284439e0130e40fa9
diff --git a/tests/block_bytecode/ls2003.mod b/tests/block_bytecode/ls2003.mod
index 3fa05a9c265ff7d4a052b557e4147eeec57c7724..db48f2058eb16ad007caaf75bbdf07b2fe34baba 100644
--- a/tests/block_bytecode/ls2003.mod
+++ b/tests/block_bytecode/ls2003.mod
@@ -215,4 +215,18 @@ shocks;
end;
perfect_foresight_setup(periods=20);
-perfect_foresight_solver(no_homotopy, stack_solve_algo = @{stack_solve_algo});
+perfect_foresight_solver(no_homotopy, stack_solve_algo = @{stack_solve_algo}
+@#if length(preconditioner) > 0
+, preconditioner = @{preconditioner}
+@#endif
+@#if preconditioner == "iterstack"
+@# if !(block && stack_solve_algo == 3)
+// Ensure that problem is not too small for iterstack, and also that there is a “residual” in the block diagonal matrix (since 3 does not divide 20)
+, iterstack_nperiods = 3
+@# else
+// For some reason iterstack can’t solve this model with block decomposition + BiCGStab
+// Hence do the equivalent of a full LU
+, iterstack_nperiods = 20
+@# endif
+@#endif
+);
diff --git a/tests/block_bytecode/run_ls2003.m b/tests/block_bytecode/run_ls2003.m
index fd7018f2a629b9ce183284b5607c1e588f850331..a68c4e97828ac705aa1c0da9659f5f7e239f897e 100644
--- a/tests/block_bytecode/run_ls2003.m
+++ b/tests/block_bytecode/run_ls2003.m
@@ -1,6 +1,6 @@
-function run_ls2003(block, storage, solve_algo, stack_solve_algo)
+function run_ls2003(block, storage, solve_algo, stack_solve_algo, preconditioner)
-% Copyright © 2010-2013 Dynare Team
+% Copyright © 2010-2025 Dynare Team
%
% This file is part of Dynare.
%
@@ -17,17 +17,23 @@ function run_ls2003(block, storage, solve_algo, stack_solve_algo)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
- disp(['TEST: ls2003 (block=' num2str(block) ', bytecode=' ...
- num2str(storage==2) ', use_dll=' num2str(storage==1) ...
- ', solve_algo=' num2str(solve_algo) ...
- ', stack_solve_algo=' num2str(stack_solve_algo) ')...']);
+ if nargin < 5
+ preconditioner = '';
+ end
+
+ fprintf('TEST: ls2003 (block=%d, bytecode=%d, use_dll=%d, solve_algo=%d, stack_solve_algo=%d', block, storage==2, storage==1, solve_algo, stack_solve_algo)
+ if ~isempty(preconditioner)
+ fprintf(', preconditioner=%s', preconditioner)
+ end
+ fprintf(')...\n')
fid = fopen('ls2003_tmp.mod', 'w');
assert(fid > 0);
fprintf(fid, ['@#define block = %d\n@#define bytecode = %d\n' ...
'@#define use_dll = %d\n' ...
'@#define solve_algo = %d\n@#define stack_solve_algo = %d\n' ...
+ '@#define preconditioner = "%s"\n' ...
'@#include \"ls2003.mod\"\n'], block, storage==2, storage==1, ...
- solve_algo, stack_solve_algo);
+ solve_algo, stack_solve_algo, preconditioner);
fclose(fid);
dynare('ls2003_tmp.mod','console')
end
diff --git a/tests/run_block_bytecode_tests.m b/tests/run_block_bytecode_tests.m
index 0c95b794ab539939b7c58eefbb4d0d51b730c674..cf65a1868ba99e68d08e7d72879850feeaf5877f 100644
--- a/tests/run_block_bytecode_tests.m
+++ b/tests/run_block_bytecode_tests.m
@@ -1,4 +1,4 @@
-% Copyright © 2011-2024 Dynare Team
+% Copyright © 2011-2025 Dynare Team
%
% This file is part of Dynare.
%
@@ -111,32 +111,43 @@ for blockFlag = 0:1
end
end
for i = 1:length(stack_solve_algos)
- try
- old_path = path;
- clear oo_ % Ensure that oo_.endo_simul won’t be overwritten when loading wsMat
- save wsMat
- run_ls2003(blockFlag, storageFlag, default_solve_algo, stack_solve_algos(i))
- load wsMat
- path(old_path);
- % Test against the reference simulation path
- load('test.mat','y_ref');
- diff = oo_.endo_simul - y_ref;
- if max(max(abs(diff))) > options_.dynatol.x
- failedBlock{size(failedBlock,2)+1} = ['block_bytecode' filesep 'run_ls2003.m(' num2str(blockFlag) ',' num2str(storageFlag) ',' num2str(default_solve_algo) ',' num2str(stack_solve_algos(i)) ')'];
- if isoctave
- exception.message = 'ERROR: simulation path differs from the reference path';
- else
- exception = MException('Dynare:simerr', 'ERROR: simulation path difers from the reference path');
+ if ismember(stack_solve_algos(i), [2 3])
+ if storageFlag ~= 2
+ preconditioners = {'umfiter', 'iterstack', 'ilu'};
+ else % bytecode
+ preconditioners = {'ilu'};
+ end
+ else
+ preconditioners = {''};
+ end
+ for j = 1:length(preconditioners)
+ try
+ old_path = path;
+ clear oo_ % Ensure that oo_.endo_simul won’t be overwritten when loading wsMat
+ save wsMat
+ run_ls2003(blockFlag, storageFlag, default_solve_algo, stack_solve_algos(i), preconditioners{j})
+ load wsMat
+ path(old_path);
+ % Test against the reference simulation path
+ load('test.mat','y_ref');
+ diff = oo_.endo_simul - y_ref;
+ if max(max(abs(diff))) > options_.dynatol.x
+ failedBlock{size(failedBlock,2)+1} = ['block_bytecode' filesep 'run_ls2003.m(' num2str(blockFlag) ',' num2str(storageFlag) ',' num2str(default_solve_algo) ',' num2str(stack_solve_algos(i)) ',' preconditioners{j} ')'];
+ if isoctave
+ exception.message = 'ERROR: simulation path differs from the reference path';
+ else
+ exception = MException('Dynare:simerr', 'ERROR: simulation path difers from the reference path');
+ end
+ printTestError(['block_bytecode' filesep 'run_ls2003.m(' num2str(blockFlag) ',' num2str(storageFlag) ',' num2str(default_solve_algo) ',' num2str(stack_solve_algos(i)) ',' preconditioners{j} ')'], exception);
+ clear exception
end
- printTestError(['block_bytecode' filesep 'run_ls2003.m(' num2str(blockFlag) ',' num2str(storageFlag) ',' num2str(default_solve_algo) ',' num2str(stack_solve_algos(i)) ')'], exception);
+ catch exception
+ load wsMat
+ path(old_path);
+ failedBlock{size(failedBlock,2)+1} = ['block_bytecode' filesep 'run_ls2003.m(' num2str(blockFlag) ',' num2str(storageFlag) ',' num2str(default_solve_algo) ',' num2str(stack_solve_algos(i)) ',' preconditioners{j} ')'];
+ printTestError(['block_bytecode' filesep 'run_ls2003.m(' num2str(blockFlag) ',' num2str(storageFlag) ',' num2str(default_solve_algo) ',' num2str(stack_solve_algos(i)) ',' preconditioners{j} ')'], exception);
clear exception
end
- catch exception
- load wsMat
- path(old_path);
- failedBlock{size(failedBlock,2)+1} = ['block_bytecode' filesep 'run_ls2003.m(' num2str(blockFlag) ',' num2str(storageFlag) ',' num2str(default_solve_algo) ',' num2str(stack_solve_algos(i)) ')'];
- printTestError(['block_bytecode' filesep 'run_ls2003.m(' num2str(blockFlag) ',' num2str(storageFlag) ',' num2str(default_solve_algo) ',' num2str(stack_solve_algos(i)) ')'], exception);
- clear exception
end
end
end