From 3c7b02e87ef717764f06b7430140bb672d515368 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 3 Nov 2023 11:52:54 +0100
Subject: [PATCH] =?UTF-8?q?=F0=9F=90=9B=20Bytecode+block+mfs>0:=20fix=20Ja?=
 =?UTF-8?q?cobian=20computation=20in=20=E2=80=9Cevaluate=E2=80=9D=20mode?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

When parsing the FSTPG{2,3} opcodes, the bytecode MEX internally constructs a
Jacobian with as many lines as feedback variables. But the preprocessor would
also output instructions for filling derivatives of equations corresponding to
recursive variables in the “evaluate” mode, thus leading to memory corruption
in the bytecode MEX.
---
 src/DynamicModel.cc | 7 +++----
 src/ModelTree.hh    | 8 +++++---
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 3752d019..8eb6d38e 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -2354,7 +2354,7 @@ DynamicModel::computeChainRuleJacobian()
               break;
             }
 
-          if (d != Zero)
+          if (d != Zero && eq >= nb_recursives)
             blocks_derivatives[blk][{ eq, var, lag }] = d;
         }
 
@@ -2369,9 +2369,8 @@ DynamicModel::computeChainRuleJacobian()
           for (const auto &[indices, d1] : blocks_derivatives[blk])
             {
               auto &[eq, var, lag] { indices };
-              assert(lag >= -1 && lag <= 1);
-              if (eq >= nb_recursives && var >= nb_recursives
-                  && !(one_boundary && lag != 0))
+              assert(lag >= -1 && lag <= 1 && eq >= nb_recursives);
+              if (var >= nb_recursives && !(one_boundary && lag != 0))
                 blocks_jacobian_sparse_column_major_order[blk].try_emplace({eq-nb_recursives, var-nb_recursives+static_cast<int>(!one_boundary)*(lag+1)*mfs_size}, d1);
             }
           blocks_jacobian_sparse_colptr[blk] = computeCSCColPtr(blocks_jacobian_sparse_column_major_order[blk], (one_boundary ? 1 : 3)*mfs_size);
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index aa0f6fd8..d7cdea22 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -1776,7 +1776,8 @@ ModelTree::writeBlockBytecodeHelper(BytecodeWriter &code_file, int block, tempor
                 const auto &[eq, var, lag] {indices};
                 int eqr {getBlockEquationID(block, eq)};
                 int varr {getBlockVariableID(block, var)};
-                if (eq >= block_recursive && var >= block_recursive)
+                assert(eq >= block_recursive);
+                if (var >= block_recursive)
                   {
                     if constexpr(dynamic)
                       if (lag != 0
@@ -1837,10 +1838,11 @@ ModelTree::writeBlockBytecodeHelper(BytecodeWriter &code_file, int block, tempor
       int varr {getBlockVariableID(block, var)};
       code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag};
       d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+      assert(eq >= block_recursive);
       if constexpr(dynamic)
-        code_file << FSTPG3_{eq, var, lag, getBlockJacobianEndoCol(block, var, lag)};
+        code_file << FSTPG3_{eq-block_recursive, var, lag, getBlockJacobianEndoCol(block, var, lag)};
       else
-        code_file << FSTPG2_{eq, getBlockJacobianEndoCol(block, var, lag)};
+        code_file << FSTPG2_{eq-block_recursive, getBlockJacobianEndoCol(block, var, lag)};
     }
 
   // Update jump offset for previous JMP
-- 
GitLab