diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 929079add12b208366a19ca4451c5d887a9cb333..10747727c73411cecd61dad7400b136ca05abfca 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -328,6 +328,17 @@ 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.itertol = 1e-12;
+options_.simul.itermaxit = 500;
+options_.simul.gmres_restart = 100;
+options_.simul.iterstack.maxlu = 20000; % Similar to TROLL’s ITERSTACK
+options_.simul.iterstack.relu = 0.5; % Similar to TROLL’s ITERSTACK
+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/newton_solve.m b/matlab/optimization/newton_solve.m
index de1d2e48286ebd95973d632d1d7d3873e7f9a3bb..120bb475db3258258c0f69850e0c01de34f265a2 100644
--- a/matlab/optimization/newton_solve.m
+++ b/matlab/optimization/newton_solve.m
@@ -23,7 +23,7 @@ function [x, errorflag, errorcode] = newton_solve(func, x, jacobian_flag, gstep,
% 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,6 +103,7 @@ for it = 1:maxit
if solve_algo == 6
p = fjac\fvec;
else
+ % Should be the same options as in solve_one_boundary.m and bytecode/Interpreter.cc (static case)
ilu_setup.type = 'ilutp';
ilu_setup.droptol = 1e-10;
[L1, U1] = ilu(fjac, ilu_setup);
diff --git a/matlab/optimization/solve_one_boundary.m b/matlab/optimization/solve_one_boundary.m
index e9177db7d02acb3d508cc5eb9a7ef8d566b17c3d..0bd84cff39d5dc59bb0fc9e6622c8d5d4cfecdfd 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,7 @@ 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 && (stack_solve_algo==1 || stack_solve_algo==0 || stack_solve_algo==6 || (ismember(stack_solve_algo, [2, 3]) && strcmp(options_.simul.preconditioner, 'iterstack') && size(g1, 1) < options_.simul.iterstack.maxlu))) || (~is_dynamic && options_.solve_algo==6) % Fallback to LU if block too small for iterstack
if verbose && ~is_dynamic
disp('steady: Sparse LU ')
end
@@ -203,64 +199,45 @@ 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, 'iterstack')
+ [L, U, P, Q] = iterstack_preconditioner(g1, options_);
+ 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)
+ if (stack_solve_algo==2 && is_dynamic) || (options_.solve_algo==7 && ~is_dynamic)
+ zb = gmres(P*g1*Q, P*r, options_.simul.gmres_restart, options_.simul.itertol, options_.simul.itermaxit, L, U);
else
- ya = ya + lambda*dx;
- if is_dynamic
- y(y_index_eq, it_) = ya;
- else
- y(y_index_eq) = ya';
- end
+ zb = bicgstab(P*g1*Q, P*r, options_.simul.itertol, options_.simul.itermaxit, L, U);
end
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)
+ ilu_setup.type = 'ilutp';
+ ilu_setup.droptol = 1e-10;
+ [L, U] = ilu(g1, ilu_setup);
+ 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/iterstack_preconditioner.m b/matlab/perfect-foresight-models/iterstack_preconditioner.m
new file mode 100644
index 0000000000000000000000000000000000000000..dafa60622fd26919c0c231dc37b0e730ab1db2ac
--- /dev/null
+++ b/matlab/perfect-foresight-models/iterstack_preconditioner.m
@@ -0,0 +1,41 @@
+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/>.
+
+if size(A, 1) < options_.simul.iterstack.maxlu
+ error('iterstack preconditionner: too small problem')
+end
+ny = size(A, 1) / options_.periods;
+if ny > options_.simul.iterstack.maxlu
+ error('iterstack preconditionner: too large problem; either increase maxlu or switch to another preconditionner/algorithm')
+end
+lusize = floor(options_.simul.iterstack.maxlu / ny) * ny;
+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
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 974cda2d5ccdabf650ba0a9384e0b187b6383db6..df86221031d926f368d566c83d02e97f3fb9c1ad 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,37 @@ 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')
+ if options_.stack_solve_algo == 2
+ y = gmres(P*A*Q, P*b, options_.simul.gmres_restart, options_.simul.itertol, options_.simul.itermaxit, L, U);
+ elseif options_.stack_solve_algo == 3
+ y = bicgstab(P*A*Q, P*b, options_.simul.itertol, options_.simul.itermaxit, L, U);
+ else
+ error('sim1: invalid value for options_.stack_solve_algo')
+ end
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..16840c89f426b893d1640c8b65e5c7e85f87785f 100644
--- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
+++ b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
@@ -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..67a36b559ad0e67f55874348199c8d94b7016fa2 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,31 @@ 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') && 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)
+ if stack_solve_algo == 2
+ zb = gmres(P*g1a*Q, P*b, options_.simul.gmres_restart, options_.simul.itertol, options_.simul.itermaxit, 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 = bicgstab(P*g1a*Q, P*b, options_.simul.itertol, options_.simul.itermaxit, L, U);
end
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..cf7b609dafeb4198f5eb22992b78408cc99525d2 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -30,6 +30,50 @@
#include "Interpreter.hh"
+void
+iter_solver_opts_t::set_static_defaults()
+{
+ // Should be the same options as in newton_solve.m and solve_one_boundary.m (static case)
+ itertol = mxCreateDoubleMatrix(0, 0, mxREAL);
+ itermaxit = mxCreateDoubleMatrix(0, 0, mxREAL);
+ gmres_restart = mxCreateDoubleMatrix(0, 0, mxREAL);
+
+ std::array field_names {"droptol", "type"};
+ std::array dims {static_cast<mwSize>(1)};
+ ilu = mxCreateStructArray(dims.size(), dims.data(), field_names.size(), field_names.data());
+ mxSetFieldByNumber(ilu, 0, 0, mxCreateDoubleScalar(1e-10));
+ mxSetFieldByNumber(ilu, 0, 1, mxCreateString("ilutp"));
+}
+
+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, "itertol");
+ if (field < 0)
+ mexErrMsgTxt("itertol is not a field of options_.simul");
+ itertol = mxGetFieldByNumber(simul, 0, field);
+
+ field = mxGetFieldNumber(simul, "itermaxit");
+ if (field < 0)
+ mexErrMsgTxt("itermaxit is not a field of options_.simul");
+ itermaxit = 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 +82,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 +106,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 +2305,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 +2736,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.itertol,
+ iter_solver_opts.itermaxit,
L1,
U1,
x0_m};
@@ -2776,28 +2757,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 +2796,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 +2806,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.itertol, iter_solver_opts.itermaxit, 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 +2863,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 +3819,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 +3887,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 +3960,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 +4077,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 +4236,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 +4298,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..e89461a46b58ec4b1e2c220b89f0232b62bf5cec 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* itertol;
+ mxArray* itermaxit;
+ mxArray* gmres_restart;
+ mxArray* ilu;
+
+ void set_static_defaults();
+ 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..d5e531dcc50db76fbc4b0e20bcd5a613609f711d 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_defaults();
+ 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;