Block decomposition: in the dynamic file, y and T are now vectors as in non-block mode

matlab/dr_block.m

@@ -73,7 +73,7 @@ end
 if (options_.bytecode)
     [zz, data]= bytecode('dynamic','evaluate', z, zx, M_.params, dr.ys, 1, data);
 else
-    [r, data] = feval([M_.fname '.dynamic'], options_, M_, oo_, z', zx, M_.params, dr.ys, M_.maximum_lag+1, data);
+    [r, data] = feval([M_.fname '.dynamic'], options_, M_, oo_, dynvars_from_endo_simul(z, M_.maximum_lag+1, M_), zx, M_.params, dr.ys, M_.maximum_lag+1, data);
 end
 dr.full_rank = 1;
 dr.eigval = [];

matlab/dynvars_from_endo_simul.m

@@ -0,0 +1,26 @@
+function y2 = dynvars_from_endo_simul(y, it_, M_)
+% Given the matrix y of paths for all endogenous (same format as
+% oo_.endo_simul), and an iteration number (first simulation period corresponds
+% to it_=M_.maximum_lag+1), return a vector of endogenous values in the format
+% expected by the dynamic.m file (i.e. whose indices are described by
+% M_.lead_lag_incidence)
+
+% Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
+
+y2 = y(:, it_+(-M_.maximum_endo_lag:M_.maximum_endo_lead));
+y2 = y2(find(M_.lead_lag_incidence'));

matlab/lnsrch1_wrapper_one_boundary.m

@@ -1,18 +1,17 @@
-function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, steady_state, T, it_)
+function r = lnsrch1_wrapper_one_boundary(ya, ll_index, fname, y, x, params, steady_state, T, it_)
 % 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
+%   ll_index            [vector]        M_.lead_lag_incidence(M_.maximum_endo_lag+1, :)
-%                                       the block
 %   fname               [string]        name of the file containing the block
 %                                       to simulate
-%   y                   [matrix]        All the endogenous variables of the model
+%   y                   [vector]        Dynamic 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
+%   T                   [vector]        Temporary terms
 % OUTPUTS
 %   r                   [vector]        The residuals of the current block
 %
@@ -40,6 +39,5 @@ function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, stea
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licen
 
-%reshape the input arguments of the dynamic function
+y2(nonzeros(ll_index)) = ya(find(ll_index));
-y(it_, y_index) = ya;
+[r, ~, g1]=feval(fname, y, x, params, steady_state, T, it_, false);
-[r, T, g1]=feval(fname, y, x, params, steady_state, T, it_, false);

matlab/lnsrch1_wrapper_two_boundaries.m

@@ -1,5 +1,6 @@
 function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
-                                             params, steady_state, T, periods, y_kmin, y_size)
+                                             params, steady_state, T, periods, ...
+                                             y_size, M_)
 % wrapper for solve_one_boundary m-file when it is used with a dynamic
 % model
 %
@@ -15,9 +16,10 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
 %   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
 %                                       in the current block
+%   M_                                  Model description structure
+%
 % OUTPUTS
 %   ra                  [vector]        The residuals of the current block
 %
@@ -46,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(y_kmin+1:y_kmin+periods, y_index) = reshape(ya',length(y_index),periods)';
+y(M_.maximum_lag+(1:periods), y_index) = reshape(ya',length(y_index),periods)';
 ra = NaN(periods*y_size, 1);
-for it_ = y_kmin+(1:periods)
+for it_ = M_.maximum_lag+(1:periods)
-    [ra((it_-y_kmin-1)*y_size+(1:y_size)), T, g1]=feval(fname, y, x, params, steady_state, T, 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

matlab/perfect-foresight-models/solve_block_decomposed_problem.m

@@ -63,7 +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
-            [y, T] = feval(funcname, y, oo_.exo_simul, M_.params, oo_.steady_state, T, it_, false);
+            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, :)));
         end
     elseif M_.block_structure.block(blk).Simulation_Type == 3 || ... % solveForwardSimple
            M_.block_structure.block(blk).Simulation_Type == 4 || ... % solveBackwardSimple

matlab/solve_one_boundary.m

@@ -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, g1] = feval(fname, y, x, params, steady_state, T, 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
@@ -203,9 +203,9 @@ 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, y_index_eq, fname, y, x, params, steady_state, T, 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] = feval(fname, y, x, params, steady_state, T, it_, false);
+                [~, 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';
             elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0)) || (~is_dynamic && options.solve_algo==6)
@@ -262,8 +262,8 @@ for it_=start:incr:finish
                         y(y_index_eq) = phat;
                     end
                     if is_dynamic
-                        [r, T, g1] = feval(fname, y, x, params, ...
+                        [r, T(:, it_), g1] = feval(fname, dynvars_from_endo_simul(y', it_, M), x, params, ...
-                                           steady_state, T, it_, false);
+                                                   steady_state, T(:, it_), it_, false);
                     else
                         [r, T, g1] = feval(fname, y, x, params, T);
                     end

matlab/solve_two_boundaries.m

@@ -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, g1]=feval(fname, y, x, params, steady_state, T, 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
@@ -314,7 +314,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, T, periods, y_kmin, Blck_size);
+            [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)';
         end

preprocessor

@@ -1 +1 @@
-Subproject commit 379be6ccef55dc38af04edb881bf6aa57bdecd65
+Subproject commit 79763911b2af3cd3ce368afd8a7e604a0bb43715