diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 45a406f1d6ecb197030b786f48435ef2db9729cb..40b3d30cb13e1e1791dcddeb45a25b782cc75ecc 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1461,6 +1461,8 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
<< "M_.block_structure.block(" << blk+1 << ").NNZDerivatives = " << blocks_derivatives[blk].size() << ';' << endl;
}
+ writeBlockDriverSparseIndicesHelper<true>(output);
+
output << "M_.block_structure.variable_reordered = [";
for (int i = 0; i < symbol_table.endo_nbr(); i++)
output << " " << endo_idx_block2orig[i]+1;
@@ -3101,10 +3103,14 @@ DynamicModel::computeChainRuleJacobian()
size_t nb_blocks { blocks.size() };
blocks_derivatives.resize(nb_blocks);
+ blocks_jacobian_sparse_column_major_order.resize(nb_blocks);
+ blocks_jacobian_sparse_colptr.resize(nb_blocks);
for (int blk {0}; blk < static_cast<int>(nb_blocks); blk++)
{
int nb_recursives = blocks[blk].getRecursiveSize();
+ int mfs_size { blocks[blk].mfs_size };
+ BlockSimulationType simulation_type { blocks[blk].simulation_type };
// Create a map from recursive vars to their defining (normalized) equation
map<int, BinaryOpNode *> recursive_vars;
@@ -3144,6 +3150,25 @@ DynamicModel::computeChainRuleJacobian()
if (d != Zero)
blocks_derivatives[blk][{ eq, var, lag }] = d;
}
+
+ // Compute the sparse representation of the Jacobian
+ if (simulation_type != BlockSimulationType::evaluateForward
+ && simulation_type != BlockSimulationType::evaluateBackward)
+ {
+ const bool one_boundary {simulation_type == BlockSimulationType::solveBackwardSimple
+ || simulation_type == BlockSimulationType::solveForwardSimple
+ || simulation_type == BlockSimulationType::solveBackwardComplete
+ || simulation_type == BlockSimulationType::solveForwardComplete};
+ 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))
+ blocks_jacobian_sparse_column_major_order[blk].emplace(pair{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);
+ }
}
/* Also store information and derivatives w.r.t. other types of variables
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 8bdca37458efc0fedd187f4df1f09fcf37e9a03e..2936fdbd84a66cc5abf4e0fb45fed085b0c23ade 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -128,7 +128,7 @@ private:
/* Computes the number of nonzero elements in deterministic Jacobian of
block-decomposed model */
int nzeDeterministicJacobianForBlock(int blk) const;
- //! Helper for writing the per-block dynamic files of block decomposed models
+ // Helper for writing the per-block dynamic files of block decomposed models (legacy representation)
template<ExprNodeOutputType output_type>
void writeDynamicPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const;
//! Writes the per-block dynamic files of block decomposed model (MATLAB version)
@@ -695,6 +695,8 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, temporary_ter
int nze_stochastic, int nze_deterministic, int nze_exo,
int nze_exo_det, int nze_other_endo) const
{
+ static_assert(!isSparseModelOutput(output_type));
+
BlockSimulationType simulation_type { blocks[blk].simulation_type };
int block_mfs_size { blocks[blk].mfs_size };
int block_recursive_size { blocks[blk].getRecursiveSize() };
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index b31cac943a228c52838065543baa6e0a4c2075f6..ef94feaabb12ee9496c83241b3c7f1ad7df33d7b 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -125,6 +125,14 @@ ModelTree::copyHelper(const ModelTree &m)
blocks_temporary_terms.push_back(convert_vector_tt(it));
for (const auto &it : m.blocks_temporary_terms_idxs)
blocks_temporary_terms_idxs.emplace(f(it.first), it.second);
+
+ for (const auto &it : m.blocks_jacobian_sparse_column_major_order)
+ {
+ map<pair<int, int>, expr_t, columnMajorOrderLess> v;
+ for (const auto &it2 : it)
+ v.emplace(it2.first, f(it2.second));
+ blocks_jacobian_sparse_column_major_order.push_back(v);
+ }
}
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
@@ -160,6 +168,7 @@ ModelTree::ModelTree(const ModelTree &m) :
blocks{m.blocks},
endo2block{m.endo2block},
eq2block{m.eq2block},
+ blocks_jacobian_sparse_colptr{m.blocks_jacobian_sparse_colptr},
endo2eq{m.endo2eq},
cutoff{m.cutoff},
mfs{m.mfs}
@@ -204,6 +213,8 @@ ModelTree::operator=(const ModelTree &m)
eq2block = m.eq2block;
blocks_temporary_terms.clear();
blocks_temporary_terms_idxs.clear();
+ blocks_jacobian_sparse_column_major_order.clear();
+ blocks_jacobian_sparse_colptr = m.blocks_jacobian_sparse_colptr;
endo2eq = m.endo2eq;
cutoff = m.cutoff;
mfs = m.mfs;
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 8947d052239273956ad69fec9bc3185b92620e8f..226b0e8d6e2e37f3758e67881db392bb1cc2d912 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -238,6 +238,18 @@ protected:
the vector of all temporary terms */
temporary_terms_idxs_t blocks_temporary_terms_idxs;
+ /* The nonzero values of the sparse block Jacobians in column-major order
+ (which is the natural order for Compressed Sparse Column (CSC) storage).
+ The pair of indices is (row, column). Nothing is stored for blocks of
+ type “evaluate”, since their Jacobian is not used in the sparse
+ representation (since the latter does not support the stochastic mode,
+ which only exists in the legacy representation). */
+ vector<SparseColumnMajorOrderMatrix> blocks_jacobian_sparse_column_major_order;
+ /* Column indices for the sparse block Jacobian in Compressed Sparse Column
+ (CSC) storage (corresponds to the “jc” vector in MATLAB terminology).
+ Same remark as above regarding blocks of type “evaluate”. */
+ vector<vector<int>> blocks_jacobian_sparse_colptr;
+
//! Computes derivatives
/*! \param order the derivation order
\param vars the derivation IDs w.r.t. which compute the derivatives */
@@ -296,9 +308,15 @@ protected:
void writeModelCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
// Writes per-block residuals and temporary terms (incl. for derivatives)
+ // This part is common to the legacy and sparse representations
template<ExprNodeOutputType output_type>
void writePerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
+ // Writes per-block Jacobian (sparse representation)
+ // Assumes temporary terms for derivatives are already set.
+ template<ExprNodeOutputType output_type>
+ void writeSparsePerBlockJacobianHelper(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. */
@@ -329,6 +347,10 @@ protected:
template<bool dynamic>
void writeJsonSparseIndicesHelper(ostream &output) const;
+ // Helper for writing sparse block derivatives indices in MATLAB/Octave driver file
+ template<bool dynamic>
+ void writeBlockDriverSparseIndicesHelper(ostream &output) const;
+
/* Helper for writing JSON output for residuals and derivatives.
Returns mlv and derivatives output at each derivation order. */
template<bool dynamic>
@@ -2254,6 +2276,56 @@ ModelTree::writeJsonSparseIndicesHelper(ostream &output) const
}
}
+template<bool dynamic>
+void
+ModelTree::writeBlockDriverSparseIndicesHelper(ostream &output) const
+{
+ for (int blk {0}; blk < static_cast<int>(blocks.size()); blk++)
+ {
+ const string struct_name { "M_.block_structure"s + (dynamic ? "" : "_stat") + ".block(" + to_string(blk+1) + ")." };
+
+ // Write indices for the sparse Jacobian (both naive and CSC storage)
+ output << struct_name << "g1_sparse_rowval = int32([";
+ for (const auto &[indices, d1] : blocks_jacobian_sparse_column_major_order.at(blk))
+ output << indices.first+1 << ' ';
+ output << "]);" << endl
+ << struct_name << "g1_sparse_colval = int32([";
+ for (const auto &[indices, d1] : blocks_jacobian_sparse_column_major_order.at(blk))
+ output << indices.second+1 << ' ';
+ output << "]);" << endl
+ << struct_name << "g1_sparse_colptr = int32([";
+ for (int it : blocks_jacobian_sparse_colptr.at(blk))
+ output << it+1 << ' ';
+ output << "]);" << endl;
+ }
+}
+
+template<ExprNodeOutputType output_type>
+void
+ModelTree::writeSparsePerBlockJacobianHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
+{
+ static_assert(isSparseModelOutput(output_type));
+
+ // NB: stochastic mode is currently unsupported by sparse representation
+ /* See also the comment above the definition of
+ blocks_jacobian_sparse_column_major_order and
+ blocks_jacobian_sparse_column_major_colptr */
+
+ assert(blocks[blk].simulation_type != BlockSimulationType::evaluateForward
+ && blocks[blk].simulation_type != BlockSimulationType::evaluateBackward);
+
+ int k {0};
+ for (const auto &[row_col, d1] : blocks_jacobian_sparse_column_major_order[blk])
+ {
+ output << "g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+ d1->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ output << ";" << endl;
+ k++;
+ }
+}
+
template<bool dynamic>
void
ModelTree::writeSparseModelJuliaFiles(const string &basename) const
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index c005f124cd1b29b15e86f69713a1e816bf564f79..acb299d3ec7f4f41205b4291a96ef5d86c482b75 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -846,6 +846,8 @@ StaticModel::writeBlockDriverOutput(ostream &output) const
output << "];" << endl
<< "M_.block_structure_stat.tmp_nbr = " << blocks_temporary_terms_idxs.size()
<< ";" << endl;
+
+ writeBlockDriverSparseIndicesHelper<false>(output);
}
SymbolType
@@ -911,10 +913,15 @@ void
StaticModel::computeChainRuleJacobian()
{
int nb_blocks = blocks.size();
+
blocks_derivatives.resize(nb_blocks);
+ blocks_jacobian_sparse_column_major_order.resize(nb_blocks);
+ blocks_jacobian_sparse_colptr.resize(nb_blocks);
+
for (int blk = 0; blk < nb_blocks; blk++)
{
int nb_recursives = blocks[blk].getRecursiveSize();
+ BlockSimulationType simulation_type { blocks[blk].simulation_type };
map<int, BinaryOpNode *> recursive_vars;
for (int i = 0; i < nb_recursives; i++)
@@ -926,8 +933,8 @@ StaticModel::computeChainRuleJacobian()
recursive_vars[deriv_id] = getBlockEquationExpr(blk, i);
}
- assert(blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesSimple
- && blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesComplete);
+ assert(simulation_type != BlockSimulationType::solveTwoBoundariesSimple
+ && simulation_type != BlockSimulationType::solveTwoBoundariesComplete);
int size = blocks[blk].size;
@@ -942,6 +949,19 @@ StaticModel::computeChainRuleJacobian()
blocks_derivatives[blk][{ eq, var, 0 }] = d1;
}
}
+
+ // Compute the sparse representation of the Jacobian
+ if (simulation_type != BlockSimulationType::evaluateForward
+ && simulation_type != BlockSimulationType::evaluateBackward)
+ {
+ for (const auto &[indices, d1] : blocks_derivatives[blk])
+ {
+ auto &[eq, var, lag] { indices };
+ assert(lag == 0);
+ blocks_jacobian_sparse_column_major_order[blk].emplace(pair{eq, var}, d1);
+ }
+ blocks_jacobian_sparse_colptr[blk] = computeCSCColPtr(blocks_jacobian_sparse_column_major_order[blk], size);
+ }
}
}
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index e50dc69610b11a6e0425baefc17f8a7d2f1b64a6..ae893c7aa9e1c7f0f893c97062043fa8eef663aa 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -44,7 +44,7 @@ private:
then compiles it with the per-block functions into a single MEX */
void writeStaticBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
- //! Helper for writing a per-block static file of block decomposed model
+ // Helper for writing a per-block static file of block decomposed model (legacy representation)
template<ExprNodeOutputType output_type>
void writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
@@ -178,6 +178,8 @@ template<ExprNodeOutputType output_type>
void
StaticModel::writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
{
+ static_assert(!isSparseModelOutput(output_type));
+
BlockSimulationType simulation_type { blocks[blk].simulation_type };
int block_recursive_size { blocks[blk].getRecursiveSize() };