Block decomposition: no longer use global variables for temporary terms in the dynamic M-file

Rather use a single vector as in non-block mode.

By the way, change the order of output arguments in static functions, to be
closer to the dynamic ones.
parent ea4d3f4d
......@@ -21,5 +21,5 @@ 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);']);
eval(['[r,~,~,g1] = ' M.fname '.static(b, y_all, exo, params, T);']);
g1 = full(g1);
......@@ -37,18 +37,18 @@ 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, check] = solve_one_boundary([M.fname '.block.static_' int2str(b)], ss, exo, ...
params, [], T, M.block_structure_stat.block(b).variable, n, 1, false, b, 0, options.simul.maxit, ...
options.solve_tolf, ...
options.slowc, 0, options.solve_algo, true, false, false, M, options);
[ss, T, check] = solve_one_boundary([M.fname '.block.static_' int2str(b)], ss, exo, ...
params, [], T, M.block_structure_stat.block(b).variable, n, 1, false, b, 0, options.simul.maxit, ...
options.solve_tolf, ...
options.slowc, 0, options.solve_algo, true, false, false, M, options);
if check
info = 1;
return
end
end
end
% The following updates the temporary terms vector
[r, g1, x, ~, T] = feval([M.fname '.static'], b, ss, exo, params, T);
% In particular, the following updates the temporary terms vector (for the dynare_solve case)
[r, x, T, g1] = feval([M.fname '.static'], b, ss, exo, params, T);
end
elseif options.bytecode
if options.solve_algo > 4
......
......@@ -52,12 +52,12 @@ else
% have zero residuals by construction
if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
M.block_structure_stat.block(b).Simulation_Type ~= 2
[r, ~, ~, ~, T] = feval(fh_static,b,ys,exo_ss,params,T);
[r, ~, T] = feval(fh_static,b,ys,exo_ss,params,T);
residuals(mfsb) = r;
else
%need to evaluate the recursive blocks to compute the
%temporary terms
[~, ~, ~, ~, T] = feval(fh_static,b,ys,exo_ss,params,T);
[~, ~, T] = feval(fh_static,b,ys,exo_ss,params,T);
end
end
else
......
function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, steady_state, it_)
function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, steady_state, T, it_)
% wrapper for solve_one_boundary m-file when it is used with a dynamic
% model
%
......@@ -11,6 +11,8 @@ function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, stea
% 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
% steady_state [vector] steady state of the model
% T [matrix] Temporary terms
% OUTPUTS
% r [vector] The residuals of the current block
%
......@@ -21,7 +23,7 @@ function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, stea
% none.
%
% Copyright (C) 2009-2017 Dynare Team
% Copyright (C) 2009-2020 Dynare Team
%
% This file is part of Dynare.
%
......@@ -40,4 +42,4 @@ function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, stea
%reshape the input arguments of the dynamic function
y(it_, :) = ya;
[r, y, g1, g2, g3]=feval(fname, y, x, params, steady_state, it_, 0);
[r, y, T, g1, g2, g3]=feval(fname, y, x, params, steady_state, T, it_, false);
function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
params, steady_state, periods, y_kmin, y_size,Periods)
params, steady_state, T, periods, y_kmin, y_size,Periods)
% wrapper for solve_one_boundary m-file when it is used with a dynamic
% model
%
......@@ -12,6 +12,8 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
% 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
% steady_state [vector] steady state of the model
% T [matrix] Temporary terms
% periods [int] The number of periods
% y_kmin [int] The maximum number of lag on en endogenous variables
% y_size [int] The number of endogenous variables
......@@ -26,7 +28,7 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
% none.
%
% Copyright (C) 2009-2017 Dynare Team
% Copyright (C) 2009-2020 Dynare Team
%
% This file is part of Dynare.
%
......@@ -45,5 +47,5 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
%reshape the input arguments of the dynamic function
y(y_kmin+1:y_kmin+periods, y_index) = reshape(ya',length(y_index),periods)';
[r, y, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, periods, 0, y_kmin, y_size, Periods);
[r, y, T, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, y_size, Periods);
ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*y_size, 1);
......@@ -71,7 +71,7 @@ if options_.block && ~options_.bytecode
z = zeros(M_.endo_nbr,1);
T = NaN(M_.block_structure_stat.tmp_nbr, 1);
for i = 1:length(M_.block_structure_stat.block)
[r, g, yy, var_indx, T] = feval([M_.fname '.static'],...
[r, yy, T, g, var_indx] = feval([M_.fname '.static'],...
i,...
oo_.steady_state,...
[oo_.exo_steady_state; ...
......
function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, T, ...
y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo, is_forward, is_dynamic, verbose, M, options, oo)
function [y, T, info] = solve_one_boundary(fname, y, x, params, steady_state, T, ...
y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, 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
%
......@@ -42,6 +42,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, T, ..
%
% OUTPUTS
% y [matrix] All endogenous variables of the model
% T [matrix] Temporary terms
% info [integer] >=0 no error
% <0 error
%
......@@ -93,9 +94,9 @@ for it_=start:incr:finish
g1=spalloc( Blck_size, Blck_size, nze);
while ~(cvg==1 || iter>maxit_)
if is_dynamic
[r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, it_, 0);
[r, y, T, g1, g2, g3] = feval(fname, y, x, params, steady_state, T, it_, false);
else
[r, y, g1] = feval(fname, y, x, params, T);
[r, y, T, g1] = feval(fname, y, x, params, T);
end
if ~isreal(r)
max_res=(-(max(max(abs(r))))^2)^0.5;
......@@ -204,7 +205,9 @@ for it_=start:incr:finish
p = -g1\r ;
[ya,f,r,check]=lnsrch1(y(it_,:)',f,g,p,stpmax, ...
'lnsrch1_wrapper_one_boundary',nn, ...
y_index_eq, options.solve_tolx, y_index_eq, fname, y, x, params, steady_state, it_);
y_index_eq, options.solve_tolx, y_index_eq, fname, y, x, params, steady_state, T, it_);
%% Recompute temporary terms, since they are not given as output of lnsrch1
[~, ~, T] = feval(fname, y, x, params, steady_state, T, it_, false);
dx = ya' - y(it_, :);
y(it_,:) = ya';
elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0)) || (~is_dynamic && options.solve_algo==6)
......@@ -261,10 +264,10 @@ for it_=start:incr:finish
y(y_index_eq) = phat;
end
if is_dynamic
[r, y, g1, g2, g3] = feval(fname, y, x, params, ...
steady_state, it_, 0);
[r, y, T, g1, g2, g3] = feval(fname, y, x, params, ...
steady_state, T, it_, false);
else
[r, y, g1] = feval(fname, y, x, params, T);
[r, y, T, g1] = feval(fname, y, x, params, T);
end
if max(abs(r))>=options.solve_tolf
[dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
......
function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_index, nze, periods, y_kmin_l, y_kmax_l, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo,options,M, oo)
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, lambda, 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
%
......@@ -9,6 +9,7 @@ function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_inde
% 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 [matrix] Temporary terms
% y_index [vector of int] The index of the endogenous variables of
% the block
% nze [integer] number of non-zero elements in the
......@@ -36,6 +37,7 @@ function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_inde
%
% OUTPUTS
% y [matrix] All endogenous variables of the model
% T [matrix] Temporary terms
% oo [structure] Results
%
% ALGORITHM
......@@ -82,7 +84,7 @@ Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods);
g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
reduced = 0;
while ~(cvg==1 || iter>maxit_)
[r, y, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, periods, 0, y_kmin, Blck_size,options.periods);
[r, y, T, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, Blck_size,options.periods);
preconditioner = 2;
g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
term1 = g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)';
......@@ -306,7 +308,7 @@ while ~(cvg==1 || 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, y, y_index,x, params, steady_state, periods, y_kmin, Blck_size,options.periods);
[yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, options.solve_tolx, fname, y, y_index,x, params, steady_state, T, periods, y_kmin, Blck_size,options.periods);
dx = ya - yn;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(yn',length(y_index),periods)';
end
......
Subproject commit 04b7d4386dd8a6dbd74c032e92d09e6ca0758fd6
Subproject commit faa6666abef430f777fd43d175de668611d67260
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment