From e2500679597a258fa4a9bd71e0b46df4d1c89b12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 10 Jan 2023 18:01:27 +0100 Subject: [PATCH] Use sparse representation for block-decomposed perfect foresight and steady state computation Ref. #1859 --- matlab/block_mfs_steadystate.m | 13 ++-- matlab/dynare_solve_block_or_bytecode.m | 9 ++- matlab/lnsrch1_wrapper_one_boundary.m | 45 ------------- matlab/lnsrch1_wrapper_two_boundaries.m | 55 ---------------- .../solve_block_decomposed_problem.m | 24 +++++-- matlab/solve_one_boundary.m | 63 ++++++++++++------- matlab/solve_two_boundaries.m | 32 ++++++---- 7 files changed, 93 insertions(+), 148 deletions(-) delete mode 100644 matlab/lnsrch1_wrapper_one_boundary.m delete mode 100644 matlab/lnsrch1_wrapper_two_boundaries.m diff --git a/matlab/block_mfs_steadystate.m b/matlab/block_mfs_steadystate.m index 2c1ac532a..b58ef2404 100644 --- a/matlab/block_mfs_steadystate.m +++ b/matlab/block_mfs_steadystate.m @@ -1,8 +1,8 @@ -function [r, g1] = block_mfs_steadystate(y, b, y_all, exo, params, T, M) -% Wrapper around the *_static.m file, for use with dynare_solve, -% when block_mfs option is given to steady. +function [r, g1] = block_mfs_steadystate(y, fh_static, b, y_all, exo, params, T, M) +% Wrapper around the static files, for use with dynare_solve, +% when block option is given to steady. -% Copyright © 2009-2020 Dynare Team +% Copyright © 2009-2023 Dynare Team % % This file is part of Dynare. % @@ -21,5 +21,6 @@ function [r, g1] = block_mfs_steadystate(y, b, y_all, exo, params, T, M) y_all(M.block_structure_stat.block(b).variable) = y; -eval(['[r,~,~,g1] = ' M.fname '.static(b, y_all, exo, params, T);']); -g1 = full(g1); +[~,~,r,g1] = fh_static(y_all, exo, params, M.block_structure_stat.block(b).g1_sparse_rowval, ... + M.block_structure_stat.block(b).g1_sparse_colval, ... + M.block_structure_stat.block(b).g1_sparse_colptr, T); diff --git a/matlab/dynare_solve_block_or_bytecode.m b/matlab/dynare_solve_block_or_bytecode.m index 17db87311..b5b55561e 100644 --- a/matlab/dynare_solve_block_or_bytecode.m +++ b/matlab/dynare_solve_block_or_bytecode.m @@ -22,6 +22,7 @@ x = y; if options.block && ~options.bytecode T = NaN(M.block_structure_stat.tmp_nbr, 1); for b = 1:length(M.block_structure_stat.block) + fh_static = str2func(sprintf('%s.sparse.block.static_%d', M.fname, b)); ss = x; if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ... M.block_structure_stat.block(b).Simulation_Type ~= 2 @@ -29,7 +30,7 @@ if options.block && ~options.bytecode [y, errorflag] = dynare_solve('block_mfs_steadystate', ... ss(M.block_structure_stat.block(b).variable), ... options.simul.maxit, options.solve_tolf, options.solve_tolx, ... - options, b, ss, exo, params, T, M); + options, fh_static, b, ss, exo, params, T, M); if errorflag info = 1; return @@ -37,7 +38,7 @@ if options.block && ~options.bytecode ss(M.block_structure_stat.block(b).variable) = y; else n = length(M.block_structure_stat.block(b).variable); - [ss, T, ~, check] = solve_one_boundary([M.fname '.static' ], ss, exo, ... + [ss, T, ~, check] = solve_one_boundary(fh_static, ss, exo, ... params, [], T, M.block_structure_stat.block(b).variable, n, 1, false, b, 0, options.simul.maxit, ... options.solve_tolf, ... 0, options.solve_algo, true, false, false, M, options, []); @@ -49,7 +50,9 @@ if options.block && ~options.bytecode end % Compute endogenous if the block is of type evaluate forward/backward or if there are recursive variables in a solve block. % Also update the temporary terms vector (needed for the dynare_solve case) - [~, x, T, g1] = feval([M.fname '.static'], b, ss, exo, params, T); + [x, T] = fh_static(ss, exo, params, M.block_structure_stat.block(b).g1_sparse_rowval, ... + M.block_structure_stat.block(b).g1_sparse_colval, ... + M.block_structure_stat.block(b).g1_sparse_colptr, T); end elseif options.bytecode if options.solve_algo >= 5 && options.solve_algo <= 8 diff --git a/matlab/lnsrch1_wrapper_one_boundary.m b/matlab/lnsrch1_wrapper_one_boundary.m deleted file mode 100644 index 3684b21dd..000000000 --- a/matlab/lnsrch1_wrapper_one_boundary.m +++ /dev/null @@ -1,45 +0,0 @@ -function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, blk, y, x, params, steady_state, T, it_, M_) -% wrapper for solve_one_boundary m-file when it is used with a dynamic -% model -% -% INPUTS -% ya [vector] The endogenous of the current block -% y_index [vector of int] The index of the endogenous variables of the block -% fname [string] name of the static/dynamic file -% blk [int] block number -% y [vector] All endogenous variables of the model -% x [matrix] All the exogenous variables of the model -% params [vector] All the parameters of the model -% steady_state [vector] steady state of the model -% T [vector] Temporary terms -% M_ Model description structure -% -% OUTPUTS -% r [vector] The residuals of the current block -% -% ALGORITHM -% none. -% -% SPECIAL REQUIREMENTS -% none. -% - -% Copyright © 2009-2022 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 1 && it_ < size(y, 2) + y3n = reshape(y(:, it_+(-1:1)), 3*M_.endo_nbr, 1); + elseif it_ > 1 % Purely backward model (in last period) + y3n = [ reshape(y(:, it_+(-1:0)), 2*M_.endo_nbr, 1); NaN(M_.endo_nbr, 1) ]; + elseif it_ < size(y, 2) % Purely forward model (in first period) + y3n = [ NaN(M_.endo_nbr, 1); reshape(y(:, it_+(0:1)), 2*M_.endo_nbr, 1) ]; + else % Static model + y3n = [ NaN(M_.endo_nbr, 1); y(:, it_); NaN(M_.endo_nbr, 1) ] + end + [y3n, T(:, it_)] = fh_dynamic(y3n, oo_.exo_simul(it_, :), M_.params, oo_.steady_state, ... + M_.block_structure.block(blk).g1_sparse_rowval, ... + M_.block_structure.block(blk).g1_sparse_colval, ... + M_.block_structure.block(blk).g1_sparse_colptr, T(:, it_)); + y(:, it_) = y3n(M_.endo_nbr+(1:M_.endo_nbr)); end elseif M_.block_structure.block(blk).Simulation_Type == 3 || ... % solveForwardSimple M_.block_structure.block(blk).Simulation_Type == 4 || ... % solveBackwardSimple M_.block_structure.block(blk).Simulation_Type == 6 || ... % solveForwardComplete M_.block_structure.block(blk).Simulation_Type == 7 % solveBackwardComplete is_forward = M_.block_structure.block(blk).Simulation_Type == 3 || M_.block_structure.block(blk).Simulation_Type == 6; - [y, T, oo_] = solve_one_boundary(funcname, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, is_forward, true, false, M_, options_, oo_); + [y, T, oo_] = solve_one_boundary(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, is_forward, true, false, M_, options_, oo_); elseif M_.block_structure.block(blk).Simulation_Type == 5 || ... % solveTwoBoundariesSimple M_.block_structure.block(blk).Simulation_Type == 8 % solveTwoBoundariesComplete - [y, T, oo_] = solve_two_boundaries(funcname, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).maximum_lag, M_.block_structure.block(blk).maximum_lead, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, options_, M_, oo_); + [y, T, oo_] = solve_two_boundaries(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).maximum_lag, M_.block_structure.block(blk).maximum_lead, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, options_, M_, oo_); end tmp = y(M_.block_structure.block(blk).variable, :); diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m index 9c1a9ac1b..932d50ff8 100644 --- a/matlab/solve_one_boundary.m +++ b/matlab/solve_one_boundary.m @@ -1,11 +1,10 @@ -function [y, T, oo_, info] = solve_one_boundary(fname, y, x, params, steady_state, T, ... +function [y, T, oo_, info] = solve_one_boundary(fh, y, x, params, steady_state, T, ... y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo, is_forward, is_dynamic, verbose, M, options, oo_) -% Computes the deterministic simulation of a block of equation containing -% lead or lag variables +% Computes the deterministic simulation or the steady state for a block of equations containing +% only lags or only leads (but not both). % % INPUTS -% fname [string] name of the file containing the block -% to simulate +% fh [handle] function handle to the static/dynamic file for the block % y [matrix] All the endogenous variables of the model % x [matrix] All the exogenous variables of the model % params [vector] All the parameters of the model @@ -26,11 +25,10 @@ function [y, T, oo_, info] = solve_one_boundary(fname, y, x, params, steady_stat % stack_solve_algo [integer] linear solver method used in the Newton algorithm % is_forward [logical] Whether the block has to be solved forward % If false, the block is solved backward -% is_dynamic [logical] If true, the block belongs to the dynamic file -% file and the oo_.deterministic_simulation field has to be uptated -% If false, the block belongs to the static -% file and the oo_.detereministic_simulation -% field remains unchanged +% is_dynamic [logical] If true, this is a deterministic simulation +% and the oo_.deterministic_simulation field is updated. +% If false, this is a steady state computation +% (oo_.detereministic_simulation remains unchanged). % verbose [logical] Whether iterations are to be printed % % OUTPUTS @@ -41,12 +39,8 @@ function [y, T, oo_, info] = solve_one_boundary(fname, y, x, params, steady_stat % % ALGORITHM % Newton with LU or GMRES or BicGstab for dynamic block -% -% SPECIAL REQUIREMENTS -% none. -% -% Copyright © 1996-2022 Dynare Team +% Copyright © 1996-2023 Dynare Team % % This file is part of Dynare. % @@ -63,7 +57,6 @@ function [y, T, oo_, info] = solve_one_boundary(fname, y, x, params, steady_stat % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . - Blck_size=size(y_index_eq,2); correcting_factor=0.01; ilu_setup.type='crout'; @@ -87,10 +80,15 @@ for it_=start:incr:finish g1=spalloc( Blck_size, Blck_size, nze); while ~(cvg || iter>maxit_) if is_dynamic - [r, yy, T(:, it_), g1] = feval(fname, Block_Num, dynvars_from_endo_simul(y, it_, M), x, params, steady_state, T(:, it_), it_, false); - y(:, it_) = yy(M.lead_lag_incidence(M.maximum_endo_lag+1,:)); + [yy, T(:, it_), r, g1] = fh(dynendo(y, it_, M), x(it_, :), params, steady_state, ... + 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_)); + y(:, it_) = yy(M.endo_nbr+(1:M.endo_nbr)); else - [r, y, T, g1] = feval(fname, Block_Num, y, x, params, T); + [y, T, r, g1] = fh(y, x, params, M.block_structure_stat.block(Block_Num).g1_sparse_rowval, ... + M.block_structure_stat.block(Block_Num).g1_sparse_colval, ... + M.block_structure_stat.block(Block_Num).g1_sparse_colptr, T); end if ~isreal(r) max_res=(-(max(max(abs(r))))^2)^0.5; @@ -201,12 +199,15 @@ for it_=start:incr:finish f = 0.5*r'*r; p = -g1\r ; [ya,f,r,check]=lnsrch1(ya,f,g,p,stpmax, ... - 'lnsrch1_wrapper_one_boundary',nn, ... - nn, options.solve_tolx, y_index_eq, fname, Block_Num, y, x, params, steady_state, T(:, it_), it_, M); + @lnsrch1_wrapper_one_boundary,nn, ... + nn, options.solve_tolx, y_index_eq, fh, Block_Num, y, x, params, steady_state, T(:, it_), it_, M); dx = ya - y(y_index_eq, it_); y(y_index_eq, it_) = ya; %% Recompute temporary terms, since they are not given as output of lnsrch1 - [~, ~, T(:, it_)] = feval(fname, Block_Num, dynvars_from_endo_simul(y, it_, M), x, params, steady_state, T(:, it_), it_, false); + [~, T(:, it_)] = fh(dynendo(y, it_, M), x(it_, :), params, steady_state, ... + 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) if verbose && ~is_dynamic disp('steady: Sparse LU ') @@ -319,3 +320,21 @@ if is_dynamic else info = 0; end + +function y3n = dynendo(y, it_, M) + if it_ > 1 && it_ < size(y, 2) + y3n = reshape(y(:, it_+(-1:1)), 3*M.endo_nbr, 1); + elseif it_ > 1 % Purely backward model (in last period) + y3n = [ reshape(y(:, it_+(-1:0)), 2*M.endo_nbr, 1); NaN(M_.endo_nbr, 1) ]; + elseif it_ < size(y, 2) % Purely forward model (in first period) + y3n = [ NaN(M_.endo_nbr, 1); reshape(y(:, it_+(0:1)), 2*M.endo_nbr, 1) ]; + else % Static model + y3n = [ NaN(M_.endo_nbr, 1); y(:, it_); NaN(M_.endo_nbr, 1) ] + end + +function r = lnsrch1_wrapper_one_boundary(ya, y_index, fh, Block_Num, y, x, params, steady_state, T, it_, M) + y(y_index, it_) = ya; + [~, ~, r] = fh(dynendo(y, it_, M), x(it_, :), params, steady_state, ... + 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); diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index 46b793ca7..946e34e6f 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -1,10 +1,9 @@ -function [y, T, oo]= solve_two_boundaries(fname, y, x, params, steady_state, T, y_index, nze, periods, y_kmin_l, y_kmax_l, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo,options,M, oo) +function [y, T, oo]= solve_two_boundaries(fh, y, x, params, steady_state, T, y_index, nze, periods, y_kmin_l, y_kmax_l, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo,options,M, oo) % Computes the deterministic simulation of a block of equation containing % both lead and lag variables using relaxation methods % % INPUTS -% fname [string] name of the file containing the block -% to simulate +% fh [handle] function handle to the dynamic file for the block % y [matrix] All the endogenous variables of the model % x [matrix] All the exogenous variables of the model % params [vector] All the parameters of the model @@ -40,12 +39,8 @@ function [y, T, oo]= solve_two_boundaries(fname, y, x, params, steady_state, T, % % ALGORITHM % Newton with LU or GMRES or BicGstab -% -% SPECIAL REQUIREMENTS -% none. -% -% Copyright © 1996-2022 Dynare Team +% Copyright © 1996-2023 Dynare Team % % This file is part of Dynare. % @@ -82,8 +77,11 @@ while ~(cvg || iter>maxit_) r = NaN(Blck_size, periods); g1a = spalloc(Blck_size*periods, Blck_size*periods, nze*periods); for it_ = y_kmin+(1:periods) - [r(:, it_-y_kmin), yy, T(:, it_), g1]=feval(fname, Block_Num, dynvars_from_endo_simul(y, it_, M), x, params, steady_state, T(:, it_), it_, false); - y(:, it_) = yy(M.lead_lag_incidence(M.maximum_endo_lag+1,:)); + [yy, T(:, it_), r(:, it_-y_kmin), g1]=fh(dynendo(y, it_, M), x(it_, :), params, steady_state, ... + 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_)); + y(:, it_) = yy(M.endo_nbr+(1:M.endo_nbr)); if periods == 1 g1a = g1(:, Blck_size+(1:Blck_size)); elseif it_ == y_kmin+1 @@ -321,7 +319,7 @@ while ~(cvg || iter>maxit_) g = (ra'*g1a)'; f = 0.5*ra'*ra; p = -g1a\ra; - [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, options.solve_tolx, fname, Block_Num, y, y_index,x, params, steady_state, T, periods, Blck_size, M); + [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,@lnsrch1_wrapper_two_boundaries,nn,nn, options.solve_tolx, fh, Block_Num, y, y_index,x, params, steady_state, T, periods, Blck_size, M); dx = ya - yn; y(y_index, y_kmin+(1:periods))=reshape(yn',length(y_index),periods); end @@ -352,3 +350,15 @@ oo.deterministic_simulation.iterations = iter; oo.deterministic_simulation.block(Block_Num).status = true;% Convergency obtained. oo.deterministic_simulation.block(Block_Num).error = max_res; oo.deterministic_simulation.block(Block_Num).iterations = iter; + +function y3n = dynendo(y, it_, M) + y3n = reshape(y(:, it_+(-1:1)), 3*M.endo_nbr, 1); + +function ra = lnsrch1_wrapper_two_boundaries(ya, fh, Block_Num, y, y_index, x, ... + params, steady_state, T, periods, ... + y_size, M) + y(y_index, M.maximum_lag+(1:periods)) = reshape(ya',length(y_index),periods); + ra = NaN(periods*y_size, 1); + for it_ = M.maximum_lag+(1:periods) + [~, ~, ra((it_-M.maximum_lag-1)*y_size+(1:y_size)), g1] = fh(dynendo(y, it_, M), x(it_, :), params, steady_state, 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_)); + end