From 7c1138eefb9040b70b1d17161489265c803426ea Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 7 Mar 2024 14:50:35 +0100
Subject: [PATCH] Bytecode: simplify FSTPG opcode

---
 mex/sources/bytecode/Evaluate.cc    | 15 +++++++--------
 mex/sources/bytecode/Evaluate.hh    |  8 ++++----
 mex/sources/bytecode/Interpreter.cc | 22 +++++-----------------
 mex/sources/bytecode/Interpreter.hh |  5 +++--
 preprocessor                        |  2 +-
 5 files changed, 20 insertions(+), 32 deletions(-)

diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc
index 50d7d9f720..01b7e616bd 100644
--- a/mex/sources/bytecode/Evaluate.cc
+++ b/mex/sources/bytecode/Evaluate.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2013-2023 Dynare Team
+ * Copyright © 2013-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -576,7 +576,7 @@ Evaluate::print_expression(const Evaluate::it_code_type& expr_begin,
           assign_lhs("residual(" + to_string(static_cast<FSTPR*>(*it_code)->pos + 1) + ")");
           break;
         case Tag::FSTPG:
-          assign_lhs("g1(" + to_string(static_cast<FSTPG*>(*it_code)->pos + 1) + ")");
+          assign_lhs("g1");
           break;
         case Tag::FSTPG2:
           {
@@ -991,10 +991,10 @@ Evaluate::print_expression(const Evaluate::it_code_type& expr_begin,
 void
 Evaluate::evaluateBlock(int it_, int y_kmin, double* __restrict__ y, int y_size,
                         double* __restrict__ x, int nb_row_x, double* __restrict__ params,
-                        const double* __restrict__ steady_y, double* __restrict__ u, int Per_u_,
-                        double* __restrict__ T, int T_nrows, map<int, double>& TEF,
+                        const double* __restrict__ steady_y, double& g1, double* __restrict__ u,
+                        int Per_u_, double* __restrict__ T, int T_nrows, map<int, double>& TEF,
                         map<pair<int, int>, double>& TEFD, map<tuple<int, int, int>, double>& TEFDD,
-                        double* __restrict__ r, double* __restrict__ g1, double* __restrict__ jacob,
+                        double* __restrict__ r, double* __restrict__ jacob,
                         double* __restrict__ jacob_exo, double* __restrict__ jacob_exo_det,
                         bool evaluate, bool no_derivatives)
 {
@@ -1446,11 +1446,10 @@ Evaluate::evaluateBlock(int it_, int y_kmin, double* __restrict__ y, int y_size,
           mexPrintf("FSTPG\n");
           mexEvalString("drawnow;");
 #endif
-          var = static_cast<FSTPG*>(*it_code)->pos;
-          g1[var] = Stack.top();
+          g1 = Stack.top();
 #ifdef DEBUG
           tmp_out << "=>";
-          mexPrintf(" g1[%d](%f)=%s\n", var, g1[var], tmp_out.str().c_str());
+          mexPrintf(" g1(%f)=%s\n", var, g1, tmp_out.str().c_str());
           tmp_out.str("");
 #endif
           Stack.pop();
diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh
index 282b357bda..d0291ef0a4 100644
--- a/mex/sources/bytecode/Evaluate.hh
+++ b/mex/sources/bytecode/Evaluate.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -98,10 +98,10 @@ public:
 
   void evaluateBlock(int it_, int y_kmin, double* __restrict__ y, int y_size,
                      double* __restrict__ x, int nb_row_x, double* __restrict__ params,
-                     const double* __restrict__ steady_y, double* __restrict__ u, int Per_u_,
-                     double* __restrict__ T, int T_nrows, map<int, double>& TEF,
+                     const double* __restrict__ steady_y, double& g1, double* __restrict__ u,
+                     int Per_u_, double* __restrict__ T, int T_nrows, map<int, double>& TEF,
                      map<pair<int, int>, double>& TEFD, map<tuple<int, int, int>, double>& TEFDD,
-                     double* __restrict__ r, double* __restrict__ g1, double* __restrict__ jacob,
+                     double* __restrict__ r, double* __restrict__ jacob,
                      double* __restrict__ jacob_exo, double* __restrict__ jacob_exo_det,
                      bool evaluate, bool no_derivatives);
 
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index 559e27e283..b61bebca59 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -149,7 +149,7 @@ Interpreter::solve_simple_one_periods()
                   if (verbosity >= 2)
                     mexPrintf("Reducing the path length in Newton step slowc=%f\n", slowc);
                   feclearexcept(FE_ALL_EXCEPT);
-                  y[Block_Contain[0].Variable + Per_y_] = ya - slowc * (r[0] / g1[0]);
+                  y[Block_Contain[0].Variable + Per_y_] = ya - slowc * (r[0] / g1);
                   if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
                     {
                       res1 = numeric_limits<double>::quiet_NaN();
@@ -165,7 +165,7 @@ Interpreter::solve_simple_one_periods()
       if (cvg)
         continue;
       feclearexcept(FE_ALL_EXCEPT);
-      y[Block_Contain[0].Variable + Per_y_] += -slowc * (rr / g1[0]);
+      y[Block_Contain[0].Variable + Per_y_] += -slowc * (rr / g1);
       if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
         {
           res1 = numeric_limits<double>::quiet_NaN();
@@ -183,8 +183,6 @@ Interpreter::solve_simple_one_periods()
 void
 Interpreter::solve_simple_over_periods(bool forward)
 {
-  g1 = static_cast<double*>(mxMalloc(sizeof(double)));
-  test_mxMalloc(g1, __LINE__, __FILE__, __func__, sizeof(double));
   r = static_cast<double*>(mxMalloc(sizeof(double)));
   test_mxMalloc(r, __LINE__, __FILE__, __func__, sizeof(double));
   if (steady_state)
@@ -207,7 +205,6 @@ Interpreter::solve_simple_over_periods(bool forward)
           it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
         }
     }
-  mxFree(g1);
   mxFree(r);
 }
 
@@ -280,8 +277,6 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
         }
       break;
     case BlockSimulationType::solveForwardSimple:
-      g1 = static_cast<double*>(mxMalloc(size * size * sizeof(double)));
-      test_mxMalloc(g1, __LINE__, __FILE__, __func__, size * size * sizeof(double));
       r = static_cast<double*>(mxMalloc(size * sizeof(double)));
       test_mxMalloc(r, __LINE__, __FILE__, __func__, size * sizeof(double));
       if (steady_state)
@@ -308,7 +303,6 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
                   residual[(it_ - y_kmin) * size + j] = r[j];
             }
         }
-      mxFree(g1);
       mxFree(r);
       break;
     case BlockSimulationType::solveForwardComplete:
@@ -380,8 +374,6 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
         }
       break;
     case BlockSimulationType::solveBackwardSimple:
-      g1 = static_cast<double*>(mxMalloc(size * size * sizeof(double)));
-      test_mxMalloc(g1, __LINE__, __FILE__, __func__, size * size * sizeof(double));
       r = static_cast<double*>(mxMalloc(size * sizeof(double)));
       test_mxMalloc(r, __LINE__, __FILE__, __func__, size * sizeof(double));
       if (steady_state)
@@ -408,7 +400,6 @@ Interpreter::evaluate_a_block(bool initialization, bool single_block, const stri
                   residual[(it_ - y_kmin) * size + j] = r[j];
             }
         }
-      mxFree(g1);
       mxFree(r);
       break;
     case BlockSimulationType::solveBackwardComplete:
@@ -2634,9 +2625,9 @@ Interpreter::compute_block_time(int my_Per_u_, bool evaluate, bool no_derivative
 
   try
     {
-      evaluator.evaluateBlock(it_, y_kmin, y, y_size, x, nb_row_x, params, steady_y, u, my_Per_u_,
-                              T, periods, TEF, TEFD, TEFDD, r, g1, jacob, jacob_exo, jacob_exo_det,
-                              evaluate, no_derivatives);
+      evaluator.evaluateBlock(it_, y_kmin, y, y_size, x, nb_row_x, params, steady_y, g1, u,
+                              my_Per_u_, T, periods, TEF, TEFD, TEFDD, r, jacob, jacob_exo,
+                              jacob_exo_det, evaluate, no_derivatives);
     }
   catch (FloatingPointException& e)
     {
@@ -4464,8 +4455,6 @@ Interpreter::solve_non_linear()
 void
 Interpreter::Simulate_Newton_One_Boundary(bool forward)
 {
-  g1 = static_cast<double*>(mxMalloc(size * size * sizeof(double)));
-  test_mxMalloc(g1, __LINE__, __FILE__, __func__, size * size * sizeof(double));
   r = static_cast<double*>(mxMalloc(size * sizeof(double)));
   test_mxMalloc(r, __LINE__, __FILE__, __func__, size * sizeof(double));
   iter = 0;
@@ -4520,7 +4509,6 @@ Interpreter::Simulate_Newton_One_Boundary(bool forward)
       mxFree(Ax_save);
       mxFree(b_save);
     }
-  mxFree(g1);
   mxFree(r);
 }
 
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index 8adabb5d76..26c2fc6aeb 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -132,9 +132,10 @@ private:
   int nb_row_x;
   int y_kmin, y_kmax, periods;
   double *x, *params;
+  double g1; // The derivative of a simple (single-equation) block in simulate mode
   double* u;
   double* steady_y;
-  double *g1, *r, *res;
+  double *r, *res;
   vector<mxArray*> jacobian_block, jacobian_exo_block, jacobian_det_exo_block;
   mxArray* GlobalTemporaryTerms;
   int it_;
diff --git a/preprocessor b/preprocessor
index 585dc63680..ab64435e8c 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 585dc63680bfdb5dcdc02ba1418648e7a0420b2c
+Subproject commit ab64435e8c05073b9af974fbfdc2dedfe18dcc4d
-- 
GitLab