From 0305b427a657fccb3e16198f1b4f15f8ae75dd23 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 17 Jun 2020 17:28:59 +0200
Subject: [PATCH] Block decomposition: use same orientation convention as
 oo_.endo_simul for matrix of simulated paths

Improves consistency, but also efficiency (less transpose needed).
---
 matlab/lnsrch1_wrapper_two_boundaries.m       |  4 +--
 .../solve_block_decomposed_problem.m          | 11 ++++----
 matlab/solve_one_boundary.m                   | 28 +++++++++----------
 matlab/solve_two_boundaries.m                 | 18 ++++++------
 4 files changed, 30 insertions(+), 31 deletions(-)

diff --git a/matlab/lnsrch1_wrapper_two_boundaries.m b/matlab/lnsrch1_wrapper_two_boundaries.m
index b26018b5e..4f6d903ce 100644
--- a/matlab/lnsrch1_wrapper_two_boundaries.m
+++ b/matlab/lnsrch1_wrapper_two_boundaries.m
@@ -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
diff --git a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
index 8f66bbf5c..fe0963331 100644
--- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
+++ b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
@@ -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;
diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m
index 3c69330d2..0b3530b1a 100644
--- a/matlab/solve_one_boundary.m
+++ b/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(:, 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
diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m
index 7c221fbe4..c46ac5c22 100644
--- a/matlab/solve_two_boundaries.m
+++ b/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(:, 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;