From a376d8c9fe9409d609d08e06df6f94d2afbdd96a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 17 Feb 2022 17:00:53 +0100
Subject: [PATCH] =?UTF-8?q?Fix=20steady=20state=20computation=20with=20byt?=
 =?UTF-8?q?ecode+block=20and=20solve=5Falgo=20=E2=A9=BD=204=20or=20?=
 =?UTF-8?q?=E2=A9=BE=209?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

– Temporary terms were not correctly passed between blocks
– solve_algo ⩾ 9 was incorrectly passed through bytecode own’s solver instead
  of through dynare_solve
---
 matlab/block_bytecode_mfs_steadystate.m |  6 +++---
 matlab/dynare_solve_block_or_bytecode.m | 27 +++++++++++++------------
 mex/sources/bytecode/SparseMatrix.cc    | 16 +--------------
 mex/sources/bytecode/bytecode.cc        |  5 ++++-
 tests/run_block_byte_tests_matlab.m     |  4 ++--
 tests/run_block_byte_tests_octave.m     |  4 ++--
 6 files changed, 26 insertions(+), 36 deletions(-)

diff --git a/matlab/block_bytecode_mfs_steadystate.m b/matlab/block_bytecode_mfs_steadystate.m
index 1747edd2a7..d9e275b76e 100644
--- a/matlab/block_bytecode_mfs_steadystate.m
+++ b/matlab/block_bytecode_mfs_steadystate.m
@@ -1,8 +1,8 @@
-function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, M)
+function [r, g1] = block_bytecode_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.
 
-% Copyright (C) 2009-2020 Dynare Team
+% Copyright (C) 2009-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -21,4 +21,4 @@ function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, M)
 
 indx = M.block_structure_stat.block(b).variable;
 y_all(indx) = y;
-[r, g1] = bytecode( y_all, exo, params, y_all, 1, y_all, 'evaluate', 'static', ['block = ' int2str(b) ]);
+[r, g1] = bytecode(y_all, exo, params, y_all, 1, y_all, T, 'evaluate', 'static', ['block = ' int2str(b) ]);
diff --git a/matlab/dynare_solve_block_or_bytecode.m b/matlab/dynare_solve_block_or_bytecode.m
index cb8df7a9be..ca5125ac3c 100644
--- a/matlab/dynare_solve_block_or_bytecode.m
+++ b/matlab/dynare_solve_block_or_bytecode.m
@@ -1,5 +1,5 @@
 function [x,info] = dynare_solve_block_or_bytecode(y, exo, params, options, M)
-% Copyright (C) 2010-2020 Dynare Team
+% Copyright (C) 2010-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -52,7 +52,7 @@ if options.block && ~options.bytecode
         [~, x, T, g1] = feval([M.fname '.static'], b, ss, exo, params, T);
     end
 elseif options.bytecode
-    if options.solve_algo > 4
+    if options.solve_algo >= 5 && options.solve_algo <= 8
         try
             x = bytecode('static', x, exo, params);
         catch ME
@@ -61,13 +61,14 @@ elseif options.bytecode
             return
         end
     elseif options.block
+        T = NaN(M.block_structure_stat.tmp_nbr, 1);
         for b = 1:length(M.block_structure_stat.block)
             if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
                     M.block_structure_stat.block(b).Simulation_Type ~= 2
                 [y, check] = dynare_solve('block_bytecode_mfs_steadystate', ...
                                           x(M.block_structure_stat ...
                                             .block(b).variable), ...
-                                          options, b, x, exo, params, M);
+                                          options, b, x, exo, params, T, M);
                 if check
                     %                    error(['STEADY: convergence problems in block '
                     %                    int2str(b)])
@@ -75,16 +76,16 @@ elseif options.bytecode
                     return
                 end
                 x(M.block_structure_stat.block(b).variable) = y;
-            else
-                try
-                    [nulldev, nulldev1, x] = bytecode(x, exo, params, ...
-                                                      x, 1, x, 'evaluate', 'static', ...
-                                                      ['block = ' int2str(b)]);
-                catch ME
-                    disp(ME.message);
-                    info = 1;
-                    return
-                end
+            end
+            % Compute endogenous if the block is of type evaluate forward/backward
+            % Also update the temporary terms vector (needed for the dynare_solve case)
+            try
+                [~, ~, x, T] = bytecode(x, exo, params, x, 1, x, T, 'evaluate', 'static', ...
+                                        ['block = ' int2str(b)]);
+            catch ME
+                disp(ME.message);
+                info = 1;
+                return
             end
         end
     else
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 76bc87b1be..4214d1c679 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -3787,25 +3787,13 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in
                                          + "%d)\n");
         }
     }
+
   if (print_it)
     {
       if (steady_state)
         {
           switch (solve_algo)
             {
-            case 0:
-              mexPrintf("MODEL STEADY STATE: MATLAB fsolve\n");
-              break;
-            case 1:
-              mexPrintf("MODEL STEADY STATE: MATLAB solve1\n");
-              break;
-            case 2:
-            case 4:
-              mexPrintf("MODEL STEADY STATE: block decomposition + MATLAB solve1\n");
-              break;
-            case 3:
-              mexPrintf("MODEL STEADY STATE: MATLAB csolve\n");
-              break;
             case 5:
               mexPrintf("MODEL STEADY STATE: (method=ByteCode own solver)\n");
               break;
@@ -3818,8 +3806,6 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in
             case 8:
               mexPrintf(preconditioner_print_out("MODEL STEADY STATE: (method=BiCGStab)\n", preconditioner, true).c_str());
               break;
-            default:
-              mexPrintf("MODEL STEADY STATE: (method=Unknown - %d - )\n", stack_solve_algo);
             }
         }
 
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index 2ab41b4561..d2f0f5d8e1 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2021 Dynare Team
+ * Copyright © 2007-2022 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -686,6 +686,9 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   if (stack_solve_algo == 7 && !steady_state)
     mexErrMsgTxt("Bytecode: Can't use option stack_solve_algo=7\n");
 
+  if (steady_state && !evaluate && (solve_algo < 5 || solve_algo > 8))
+    mexErrMsgTxt("Bytecode: solve_algo must be between 5 and 8 when using the internal steady state solver\n");
+
   size_t size_of_direction = col_y*row_y*sizeof(double);
   auto *y = static_cast<double *>(mxMalloc(size_of_direction));
   error_msg.test_mxMalloc(y, __LINE__, __FILE__, __func__, size_of_direction);
diff --git a/tests/run_block_byte_tests_matlab.m b/tests/run_block_byte_tests_matlab.m
index 5fe94b3d2a..de922c0322 100644
--- a/tests/run_block_byte_tests_matlab.m
+++ b/tests/run_block_byte_tests_matlab.m
@@ -1,4 +1,4 @@
-% Copyright (C) 2011-2020 Dynare Team
+% Copyright (C) 2011-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -50,7 +50,7 @@ for blockFlag = 0:1
             solve_algos = [1:4 6:9];
             stack_solve_algos = 0:4;
         else
-            solve_algos = 1:8;
+            solve_algos = 1:9;
             stack_solve_algos = 0:5;
         end
         if has_optimization_toolbox
diff --git a/tests/run_block_byte_tests_octave.m b/tests/run_block_byte_tests_octave.m
index 9c52eebbd1..b6319ca272 100644
--- a/tests/run_block_byte_tests_octave.m
+++ b/tests/run_block_byte_tests_octave.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2009-2020 Dynare Team
+## Copyright (C) 2009-2022 Dynare Team
 ##
 ## This file is part of Dynare.
 ##
@@ -52,7 +52,7 @@ for blockFlag = 0:1
             solve_algos = [0:4 6:9];
             stack_solve_algos = 0:4;
         else
-            solve_algos = 0:8;
+            solve_algos = 0:9;
             stack_solve_algos = 0:5;
         endif
 
-- 
GitLab