diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 7efae9c9f3c6decfca1cd5e91a46cfbabda33761..8232d22a5a69253f51c05ab77a19a2cf00c50f02 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -665,87 +665,11 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, temporary_ter
int nze_exo_det, int nze_other_endo) const
{
BlockSimulationType simulation_type { blocks[blk].simulation_type };
- int block_size { blocks[blk].size };
int block_mfs_size { blocks[blk].mfs_size };
int block_recursive_size { blocks[blk].getRecursiveSize() };
- deriv_node_temp_terms_t tef_terms;
-
- auto write_eq_tt = [&](int eq)
- {
- for (auto it : blocks_temporary_terms[blk][eq])
- {
- if (dynamic_cast<AbstractExternalFunctionNode *>(it))
- it->writeExternalFunctionOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
-
- output << " ";
- it->writeOutput(output, output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms);
- output << '=';
- it->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
- temporary_terms.insert(it);
- output << ';' << endl;
- }
- };
-
- // The equations
- for (int eq {0}; eq < block_size; eq++)
- {
- write_eq_tt(eq);
-
- EquationType equ_type { getBlockEquationType(blk, eq) };
- BinaryOpNode *e { getBlockEquationExpr(blk, eq) };
- expr_t lhs { e->arg1 }, rhs { e->arg2 };
- switch (simulation_type)
- {
- case BlockSimulationType::evaluateBackward:
- case BlockSimulationType::evaluateForward:
- evaluation:
- if (equ_type == EquationType::evaluateRenormalized)
- {
- e = getBlockEquationRenormalizedExpr(blk, eq);
- lhs = e->arg1;
- rhs = e->arg2;
- }
- else if (equ_type != EquationType::evaluate)
- {
- cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl;
- exit(EXIT_FAILURE);
- }
- output << " ";
- lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- output << '=';
- rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- output << ';' << endl;
- break;
- case BlockSimulationType::solveBackwardSimple:
- case BlockSimulationType::solveForwardSimple:
- case BlockSimulationType::solveBackwardComplete:
- case BlockSimulationType::solveForwardComplete:
- case BlockSimulationType::solveTwoBoundariesComplete:
- case BlockSimulationType::solveTwoBoundariesSimple:
- if (eq < block_recursive_size)
- goto evaluation;
- output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
- << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type)
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=(";
- goto end;
- default:
- end:
- lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- output << ")-(";
- rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- output << ");" << endl;
- }
- }
-
- // The Jacobian if we have to solve the block
-
- /* Write temporary terms for derivatives.
- Also note that in the case of “evaluate” blocks, derivatives are not
- computed in deterministic mode; still their temporary terms must be
- computed even in that mode, because they may be needed in subsequent
- blocks. */
- write_eq_tt(blocks[blk].size);
+ // Write residuals and temporary terms (incl. for derivatives)
+ writePerBlockHelper<output_type>(blk, output, temporary_terms);
if constexpr(isCOutput(output_type))
output << " if (stochastic_mode) {" << endl;
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index bcb4a6d89e26b809856040ed42fdfafe677bfd05..a03e96597f35b7d0dbdb7ad8b789d5591d423a6f 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -271,6 +271,10 @@ protected:
template<ExprNodeOutputType output_type>
pair<vector<ostringstream>, vector<ostringstream>> writeModelFileHelper() const;
+ // Writes per-block residuals and temporary terms (incl. for derivatives)
+ template<ExprNodeOutputType output_type>
+ void writePerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
+
/* Helper for writing derivatives w.r.t. parameters.
Returns { tt, rp, gp, rpp, gpp, hp, g3p }.
g3p is empty if requesting a static output type. */
@@ -812,6 +816,89 @@ ModelTree::writeModelFileHelper() const
return { move(d_output), move(tt_output) };
}
+template<ExprNodeOutputType output_type>
+void
+ModelTree::writePerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
+{
+ int block_recursive_size { blocks[blk].getRecursiveSize() };
+
+ // The equations
+ deriv_node_temp_terms_t tef_terms;
+
+ auto write_eq_tt = [&](int eq)
+ {
+ for (auto it : blocks_temporary_terms[blk][eq])
+ {
+ if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+ it->writeExternalFunctionOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
+
+ output << " ";
+ it->writeOutput(output, output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms);
+ output << '=';
+ it->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
+ temporary_terms.insert(it);
+ output << ';' << endl;
+ }
+ };
+
+ for (int eq {0}; eq < blocks[blk].size; eq++)
+ {
+ write_eq_tt(eq);
+
+ EquationType equ_type { getBlockEquationType(blk, eq) };
+ BinaryOpNode *e { getBlockEquationExpr(blk, eq) };
+ expr_t lhs { e->arg1 }, rhs { e->arg2 };
+ switch (blocks[blk].simulation_type)
+ {
+ case BlockSimulationType::evaluateBackward:
+ case BlockSimulationType::evaluateForward:
+ evaluation:
+ if (equ_type == EquationType::evaluateRenormalized)
+ {
+ e = getBlockEquationRenormalizedExpr(blk, eq);
+ lhs = e->arg1;
+ rhs = e->arg2;
+ }
+ else if (equ_type != EquationType::evaluate)
+ {
+ cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl;
+ exit(EXIT_FAILURE);
+ }
+ output << " ";
+ lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ output << '=';
+ rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ output << ';' << endl;
+ break;
+ case BlockSimulationType::solveBackwardSimple:
+ case BlockSimulationType::solveForwardSimple:
+ case BlockSimulationType::solveBackwardComplete:
+ case BlockSimulationType::solveForwardComplete:
+ case BlockSimulationType::solveTwoBoundariesComplete:
+ case BlockSimulationType::solveTwoBoundariesSimple:
+ if (eq < block_recursive_size)
+ goto evaluation;
+ output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=(";
+ lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ output << ")-(";
+ rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ output << ");" << endl;
+ break;
+ default:
+ cerr << "Incorrect type for block " << blk+1 << endl;
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ /* Write temporary terms for derivatives.
+ This is done even for “evaluate” blocks, whose derivatives are not
+ always computed at runtime; still those temporary terms may be needed by
+ subsequent blocks. */
+ write_eq_tt(blocks[blk].size);
+}
+
template<ExprNodeOutputType output_type>
tuple<ostringstream, ostringstream, ostringstream, ostringstream,
ostringstream, ostringstream, ostringstream>
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index a67637d7a5f0745d2080ecbe08b2ef22347ab25f..092ebdc046af0fb1ef5feebe99504963d6f088e6 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -176,82 +176,8 @@ StaticModel::writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms
BlockSimulationType simulation_type { blocks[blk].simulation_type };
int block_recursive_size { blocks[blk].getRecursiveSize() };
- // The equations
- deriv_node_temp_terms_t tef_terms;
-
- auto write_eq_tt = [&](int eq)
- {
- for (auto it : blocks_temporary_terms[blk][eq])
- {
- if (dynamic_cast<AbstractExternalFunctionNode *>(it))
- it->writeExternalFunctionOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
-
- output << " ";
- it->writeOutput(output, output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms);
- output << '=';
- it->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
- temporary_terms.insert(it);
- output << ';' << endl;
- }
- };
-
- for (int eq {0}; eq < blocks[blk].size; eq++)
- {
- write_eq_tt(eq);
-
- EquationType equ_type { getBlockEquationType(blk, eq) };
- BinaryOpNode *e { getBlockEquationExpr(blk, eq) };
- expr_t lhs { e->arg1 }, rhs { e->arg2 };
- switch (simulation_type)
- {
- case BlockSimulationType::evaluateBackward:
- case BlockSimulationType::evaluateForward:
- evaluation:
- if (equ_type == EquationType::evaluateRenormalized)
- {
- e = getBlockEquationRenormalizedExpr(blk, eq);
- lhs = e->arg1;
- rhs = e->arg2;
- }
- else if (equ_type != EquationType::evaluate)
- {
- cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl;
- exit(EXIT_FAILURE);
- }
- output << " ";
- lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- output << '=';
- rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- output << ';' << endl;
- break;
- case BlockSimulationType::solveBackwardSimple:
- case BlockSimulationType::solveForwardSimple:
- case BlockSimulationType::solveBackwardComplete:
- case BlockSimulationType::solveForwardComplete:
- if (eq < block_recursive_size)
- goto evaluation;
- output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
- << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type)
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=(";
- lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- output << ")-(";
- rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- output << ");" << endl;
- break;
- default:
- cerr << "Incorrect type for block " << blk+1 << endl;
- exit(EXIT_FAILURE);
- }
- }
-
- /* Write temporary terms for derivatives.
- This is done even for “evaluate” blocks, whose derivatives are not
- computed at runtime; still those temporary terms may be needed by
- subsequent blocks (not calling write_eq_tt() would not be a bug though,
- because those terms would not be added to temporary_terms_union and would
- therefore not be used; still, it’s probably better performance-wise to
- use those temporary terms). */
- write_eq_tt(blocks[blk].size);
+ // Write residuals and temporary terms (incl. for derivatives)
+ writePerBlockHelper<output_type>(blk, output, temporary_terms);
// The Jacobian if we have to solve the block
if (simulation_type != BlockSimulationType::evaluateBackward