stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Block decomposition: use same orientation convention as oo_.endo_simul for matrix of simulated paths 

Improves consistency, but also efficiency (less transpose needed).
time-shift
Sébastien Villemot 2020-06-17 17:28:59 +02:00
parent 20431ed312
commit 0305b427a6
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
4 changed files with 30 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
matlab/lnsrch1_wrapper_two_boundaries.m
View File

@ -48,8 +48,8 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
% along with Dynare. If not, see <http://www.gnu.org/licen
%reshape the input arguments of the dynamic function
y(M_.maximum_lag+(1:periods), y_index) = reshape(ya',length(y_index),periods)';
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]=feval(fname, dynvars_from_endo_simul(y', it_, M_), x, params, steady_state, T(:, it_), it_, false);
[ra((it_-M_.maximum_lag-1)*y_size+(1:y_size)), ~, g1]=feval(fname, dynvars_from_endo_simul(y, it_, M_), x, params, steady_state, T(:, it_), it_, false);
end

11
matlab/perfect-foresight-models/solve_block_decomposed_problem.m
View File

@ -39,7 +39,7 @@ if options_.verbosity
skipline()
end
y=oo_.endo_simul';
y=oo_.endo_simul;
T=NaN(M_.block_structure.dyn_tmp_nbr, options_.periods+M_.maximum_lag+M_.maximum_lead);
oo_.deterministic_simulation.status = 0;
@ -63,9 +63,9 @@ for blk = 1:length(M_.block_structure.block)
range = M_.maximum_lag+options_.periods:-1:M_.maximum_lag+1;
end
for it_ = range
y2 = dynvars_from_endo_simul(y', it_, M_);
y2 = dynvars_from_endo_simul(y, it_, M_);
[y2, T(:, it_)] = feval(funcname, y2, oo_.exo_simul, M_.params, oo_.steady_state, T(:, it_), it_, false);
y(it_, find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, :))) = y2(nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+1, :)));
y(find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, :)), it_) = y2(nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+1, :)));
end
elseif M_.block_structure.block(blk).Simulation_Type == 3 || ... % solveForwardSimple
M_.block_structure.block(blk).Simulation_Type == 4 || ... % solveBackwardSimple
@ -77,7 +77,7 @@ for blk = 1:length(M_.block_structure.block)
[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_.solve_tolf, options_.slowc, cutoff, options_.stack_solve_algo, options_, M_, oo_);
end
tmp = y(:,M_.block_structure.block(blk).variable);
tmp = y(M_.block_structure.block(blk).variable, :);
if any(isnan(tmp) | isinf(tmp))
disp(['Inf or Nan value during the resolution of block ' num2str(blk)]);
oo_.deterministic_simulation.status = false;
@ -90,5 +90,4 @@ for blk = 1:length(M_.block_structure.block)
end
end
oo_.endo_simul = y';
oo_.endo_simul = y;

28
matlab/solve_one_boundary.m
View File

@ -92,7 +92,7 @@ for it_=start:incr:finish
g1=spalloc( Blck_size, Blck_size, nze);
while ~(cvg==1 || iter>maxit_)
if is_dynamic
[r, T(:, it_), g1] = feval(fname, dynvars_from_endo_simul(y', it_, M), x, params, steady_state, T(:, it_), it_, false);
[r, T(:, it_), g1] = feval(fname, dynvars_from_endo_simul(y, it_, M), x, params, steady_state, T(:, it_), it_, false);
else
[r, T, g1] = feval(fname, y, x, params, T);
end
@ -104,7 +104,7 @@ for it_=start:incr:finish
if verbose
disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)])
if is_dynamic
disp([M.endo_names{y_index_eq} num2str([y(it_,y_index_eq)' r g1])])
disp([M.endo_names{y_index_eq} num2str([y(y_index_eq, it_) r g1])])
else
disp([M.endo_names{y_index_eq} num2str([y(y_index_eq) r g1])])
end
@ -135,7 +135,7 @@ for it_=start:incr:finish
disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')])
end
dx = - r/(g1+correcting_factor*speye(Blck_size));
y(it_,y_index_eq)=ya_save+lambda*dx;
y(y_index_eq, it_)=ya_save+lambda*dx;
continue
else
if verbose
@ -152,7 +152,7 @@ for it_=start:incr:finish
disp(['reducing the path length: lambda=' num2str(lambda,'%f')])
end
if is_dynamic
y(it_,y_index_eq)=ya_save-lambda*dx;
y(y_index_eq, it_)=ya_save-lambda*dx;
else
y(y_index_eq)=ya_save-lambda*dx;
end
@ -183,7 +183,7 @@ for it_=start:incr:finish
end
end
if is_dynamic
ya = y(it_,y_index_eq)';
ya = y(y_index_eq, it_);
else
ya = y(y_index_eq);
end
@ -203,11 +203,11 @@ for it_=start:incr:finish
p = -g1\r ;
[ya,f,r,check]=lnsrch1(ya,f,g,p,stpmax, ...
'lnsrch1_wrapper_one_boundary',nn, ...
y_index_eq, options.solve_tolx, M.lead_lag_incidence(M.maximum_endo_lag+1, :), fname, dynvars_from_endo_simul(y', it_, M), x, params, steady_state, T(:, it_), it_);
y_index_eq, options.solve_tolx, M.lead_lag_incidence(M.maximum_endo_lag+1, :), fname, dynvars_from_endo_simul(y, it_, M), x, params, steady_state, T(:, it_), it_);
%% Recompute temporary terms, since they are not given as output of lnsrch1
[~, T(:, it_)] = feval(fname, dynvars_from_endo_simul(y', it_, M), x, params, steady_state, T(:, it_), it_, false);
dx = ya' - y(it_, index_eq);
y(it_, index_eq) = ya';
[~, T(:, it_)] = feval(fname, dynvars_from_endo_simul(y, it_, M), x, params, steady_state, T(:, it_), it_, false);
dx = ya' - y(index_eq, it_);
y(index_eq, it_) = ya;
elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0)) || (~is_dynamic && options.solve_algo==6)
if verbose && ~is_dynamic
disp('steady: Sparse LU ')
@ -215,7 +215,7 @@ for it_=start:incr:finish
dx = g1\r;
ya = ya - lambda*dx;
if is_dynamic
y(it_,y_index_eq) = ya';
y(y_index_eq, it_) = ya;
else
y(y_index_eq) = ya;
end
@ -242,7 +242,7 @@ for it_=start:incr:finish
else
ya = ya + lambda*dx;
if is_dynamic
y(it_,y_index_eq) = ya';
y(y_index_eq, it_) = ya;
else
y(y_index_eq) = ya';
end
@ -257,12 +257,12 @@ for it_=start:incr:finish
[L1, U1]=ilu(g1,ilu_setup);
phat = ya - U1 \ (L1 \ r);
if is_dynamic
y(it_,y_index_eq) = phat;
y(y_index_eq, it_) = phat;
else
y(y_index_eq) = phat;
end
if is_dynamic
[r, T(:, it_), g1] = feval(fname, dynvars_from_endo_simul(y', it_, M), x, params, ...
[r, T(:, it_), g1] = feval(fname, dynvars_from_endo_simul(y, it_, M), x, params, ...
steady_state, T(:, it_), it_, false);
else
[r, T, g1] = feval(fname, y, x, params, T);
@ -288,7 +288,7 @@ for it_=start:incr:finish
else
ya = ya + lambda*dx;
if is_dynamic
y(it_,y_index_eq) = ya';
y(y_index_eq, it_) = ya;
else
y(y_index_eq) = ya';
end

18
matlab/solve_two_boundaries.m
View File

@ -83,7 +83,7 @@ while ~(cvg==1 || 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), T(:, it_), g1]=feval(fname, dynvars_from_endo_simul(y', it_, M), x, params, steady_state, T(:, it_), it_, false);
[r(:, it_-y_kmin), T(:, it_), g1]=feval(fname, dynvars_from_endo_simul(y, it_, M), x, params, steady_state, T(:, it_), it_, false);
if periods == 1
g1a = g1(:, Blck_size+(1:Blck_size));
elseif it_ == y_kmin+1
@ -95,7 +95,7 @@ while ~(cvg==1 || iter>maxit_)
end
end
preconditioner = 2;
ya = reshape(y(y_kmin+(1:periods),y_index)', 1, periods*Blck_size)';
ya = reshape(y(y_index, y_kmin+(1:periods)), 1, periods*Blck_size)';
ra = reshape(r, periods*Blck_size, 1);
b=-ra+g1a*ya;
[max_res, max_indx]=max(max(abs(r')));
@ -131,7 +131,7 @@ while ~(cvg==1 || iter>maxit_)
disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]);
end
dx = (g1aa+correcting_factor*speye(periods*Blck_size))\ba- ya_save;
y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';
y(y_index, y_kmin+(1:periods))=reshape((ya_save+lambda*dx)',length(y_index),periods);
continue
else
disp('The singularity of the jacobian matrix could not be corrected');
@ -144,7 +144,7 @@ while ~(cvg==1 || iter>maxit_)
if verbose
disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
end
y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';
y(y_index, y_kmin+(1:periods))=reshape((ya_save+lambda*dx)',length(y_index),periods);
continue
else
if verbose
@ -175,7 +175,7 @@ while ~(cvg==1 || iter>maxit_)
if stack_solve_algo==0
dx = g1a\b- ya;
ya = ya + lambda*dx;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
elseif stack_solve_algo==1
for t=1:periods
first_elem = (t-1)*Blck_size+1;
@ -213,7 +213,7 @@ while ~(cvg==1 || iter>maxit_)
y_Elem = Blck_size * (t-1)+1:Blck_size * (t);
dx(y_Elem) = za - ya(y_Elem);
ya(y_Elem) = ya(y_Elem) + lambda*dx(y_Elem);
y(y_kmin + t, y_index) = ya(y_Elem);
y(y_index, y_kmin + t) = ya(y_Elem);
end
elseif stack_solve_algo==2
flag1=1;
@ -261,7 +261,7 @@ while ~(cvg==1 || iter>maxit_)
else
dx = za - ya;
ya = ya + lambda*dx;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
end
end
elseif stack_solve_algo==3
@ -304,7 +304,7 @@ while ~(cvg==1 || iter>maxit_)
else
dx = za - ya;
ya = ya + lambda*dx;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
end
end
elseif stack_solve_algo==4
@ -316,7 +316,7 @@ while ~(cvg==1 || iter>maxit_)
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, T, periods, Blck_size, M);
dx = ya - yn;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(yn',length(y_index),periods)';
y(y_index, y_kmin+(1:periods))=reshape(yn',length(y_index),periods);
end
end
iter=iter+1;