From db8c5439f7d0e82837d645bf3d52b8489b66563c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 23 Jun 2020 17:50:50 +0200
Subject: [PATCH] use_dll: higher order derivatives are now returned as sparse
matrices by static/dynamic files
Previously they were returned as 3-column matrices. But this was inconsistent
with the M-file mode.
---
src/DynamicModel.cc | 104 +++++++++++++++++++++++++++++---------------
src/ModelTree.cc | 11 -----
src/ModelTree.hh | 4 --
src/StaticModel.cc | 87 ++++++++++++++++++++++--------------
4 files changed, 124 insertions(+), 82 deletions(-)
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index fb809507..ea237b54 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1476,18 +1476,40 @@ DynamicModel::writeDynamicCFile(const string &basename) const
<< endl
<< " if (nlhs >= 3)" << endl
<< " {" << endl
- << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3 << ", mxREAL);" << endl
- << " double *v2 = mxGetPr(plhs[2]);" << endl
+ << " mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+ << " mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+ << " mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
<< " dynamic_g2_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl
- << " dynamic_g2(y, x, nb_row_x, params, steady_state, it_, T, v2);" << endl
+ << " dynamic_g2(y, x, nb_row_x, params, steady_state, it_, T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
+ << " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
+ << " mxArray *n = mxCreateDoubleScalar(" << dynJacobianColsNbr*dynJacobianColsNbr << ");" << endl
+ << " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
+ << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
+ << " plhs[2] = plhs_sparse[0];" << endl
+ << " mxDestroyArray(g2_i);" << endl
+ << " mxDestroyArray(g2_j);" << endl
+ << " mxDestroyArray(g2_v);" << endl
+ << " mxDestroyArray(m);" << endl
+ << " mxDestroyArray(n);" << endl
<< " }" << endl
<< endl
<< " if (nlhs >= 4)" << endl
<< " {" << endl
- << " plhs[3] = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 3 << ", mxREAL);" << endl
- << " double *v3 = mxGetPr(plhs[3]);" << endl
+ << " mxArray *g3_i = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
+ << " mxArray *g3_j = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
+ << " mxArray *g3_v = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
<< " dynamic_g3_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl
- << " dynamic_g3(y, x, nb_row_x, params, steady_state, it_, T, v3);" << endl
+ << " dynamic_g3(y, x, nb_row_x, params, steady_state, it_, T, mxGetPr(g3_i), mxGetPr(g3_j), mxGetPr(g3_v));" << endl
+ << " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
+ << " mxArray *n = mxCreateDoubleScalar(" << dynJacobianColsNbr*dynJacobianColsNbr*dynJacobianColsNbr << ");" << endl
+ << " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g3_i, g3_j, g3_v, m, n };" << endl
+ << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
+ << " plhs[3] = plhs_sparse[0];" << endl
+ << " mxDestroyArray(g3_i);" << endl
+ << " mxDestroyArray(g3_j);" << endl
+ << " mxDestroyArray(g3_v);" << endl
+ << " mxDestroyArray(m);" << endl
+ << " mxDestroyArray(n);" << endl
<< " }" << endl
<< endl
<< " free(T);" << endl
@@ -1951,7 +1973,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
then the third. This gives a significant performance boost in use_dll
mode (at both compilation and runtime), because it facilitates memory
accesses and expression reusage. */
- ostringstream col0_output, col1_output, col2_output;
+ ostringstream i_output, j_output, v_output;
int k = 0; // Current line index in the 3-column matrix
for (const auto &[vidx, d] : derivatives[i])
@@ -1973,16 +1995,19 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
}
else
{
- sparseHelper(i, col0_output, k, 0, output_type);
- col0_output << "=" << eq + 1 << ";" << endl;
-
- sparseHelper(i, col1_output, k, 1, output_type);
- col1_output << "=" << col_idx + 1 << ";" << endl;
-
- sparseHelper(i, col2_output, k, 2, output_type);
- col2_output << "=";
- d->writeOutput(col2_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms);
- col2_output << ";" << endl;
+ i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << "=" << eq + 1 << ";" << endl;
+ j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << "=" << col_idx + 1 << ";" << endl;
+ v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+ d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms);
+ v_output << ";" << endl;
k++;
}
@@ -1997,23 +2022,27 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
<< "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
else
{
- sparseHelper(2, col0_output, k, 0, output_type);
- col0_output << "=" << eq + 1 << ";" << endl;
-
- sparseHelper(2, col1_output, k, 1, output_type);
- col1_output << "=" << col_idx_sym + 1 << ";" << endl;
-
- sparseHelper(2, col2_output, k, 2, output_type);
- col2_output << "=";
- sparseHelper(2, col2_output, k-1, 2, output_type);
- col2_output << ";" << endl;
+ i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << "=" << eq + 1 << ";" << endl;
+ j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << "=" << col_idx_sym + 1 << ";" << endl;
+ v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="
+ << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
k++;
}
}
}
if (output_type != ExprNodeOutputType::juliaDynamicModel)
- d_output[i] << col0_output.str() << col1_output.str() << col2_output.str();
+ d_output[i] << i_output.str() << j_output.str() << v_output.str();
}
if (output_type == ExprNodeOutputType::matlabDynamicModel)
@@ -2052,18 +2081,17 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
ncols *= dynJacobianColsNbr;
ntt += temporary_terms_derivatives[i].size();
string gname = "g" + to_string(i);
- string vname = "v" + to_string(i);
string gprevname = "g" + to_string(i-1);
init_output.str("");
end_output.str("");
if (derivatives[i].size())
{
- init_output << vname << " = zeros(" << NNZDerivatives[i] << ",3);";
+ init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl
+ << gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl
+ << gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl;
end_output << gname << " = sparse("
- << vname << "(:,1),"
- << vname << "(:,2),"
- << vname << "(:,3),"
+ << gname << "_i," << gname << "_j," << gname << "_v,"
<< nrows << "," << ncols << ");";
}
else
@@ -2085,13 +2113,19 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
for (size_t i = 0; i < d_output.size(); i++)
{
string funcname = i == 0 ? "resid" : "g" + to_string(i);
- string argname = i == 0 ? "residual" : i == 1 ? "g1" : "v" + to_string(i);
DynamicOutput << "void dynamic_" << funcname << "_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T)" << endl
<< "{" << endl
<< tt_output[i].str()
<< "}" << endl
<< endl
- << "void dynamic_" << funcname << "(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *" << argname << ")" << endl
+ << "void dynamic_" << funcname << "(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, ";
+ if (i == 0)
+ DynamicOutput << "double *residual";
+ else if (i == 1)
+ DynamicOutput << "double *g1";
+ else
+ DynamicOutput << "double *" << funcname << "_i, double *" << funcname << "_j, double *" << funcname << "_v";
+ DynamicOutput << ")" << endl
<< "{" << endl;
if (i == 0)
DynamicOutput << " double lhs, rhs;" << endl;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 0774c321..856c2260 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1792,17 +1792,6 @@ ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutput
output << RIGHT_ARRAY_SUBSCRIPT(output_type);
}
-void
-ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
-{
- output << "v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
- if (isMatlabOutput(output_type) || isJuliaOutput(output_type))
- output << row_nb + 1 << "," << col_nb + 1;
- else
- output << row_nb + col_nb * NNZDerivatives[order];
- output << RIGHT_ARRAY_SUBSCRIPT(output_type);
-}
-
void
ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
{
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 518d3bdb..0ea39d5c 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -456,10 +456,6 @@ public:
//! Helper for writing the Jacobian elements in MATLAB and C
/*! Writes either (i+1,j+1) or [i+j*no_eq] */
void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
- //! Helper for writing the sparse Hessian or third derivatives in MATLAB and C
- /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[2]]
- If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[3]] */
- void sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Returns all the equation tags associated to an equation
inline map<string, string>
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 0cfab15e..ae5ff71c 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1301,7 +1301,7 @@ StaticModel::writeStaticModel(const string &basename,
then the third. This gives a significant performance boost in use_dll
mode (at both compilation and runtime), because it facilitates memory
accesses and expression reusage. */
- ostringstream col0_output, col1_output, col2_output;
+ ostringstream i_output, j_output, v_output;
int k = 0; // Current line index in the 3-column matrix
for (const auto &[vidx, d] : derivatives[i])
@@ -1323,16 +1323,19 @@ StaticModel::writeStaticModel(const string &basename,
}
else
{
- sparseHelper(i, col0_output, k, 0, output_type);
- col0_output << "=" << eq + 1 << ";" << endl;
-
- sparseHelper(i, col1_output, k, 1, output_type);
- col1_output << "=" << col_idx + 1 << ";" << endl;
-
- sparseHelper(i, col2_output, k, 2, output_type);
- col2_output << "=";
- d->writeOutput(col2_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms);
- col2_output << ";" << endl;
+ i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << "=" << eq + 1 << ";" << endl;
+ j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << "=" << col_idx + 1 << ";" << endl;
+ v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+ d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms);
+ v_output << ";" << endl;
k++;
}
@@ -1347,23 +1350,27 @@ StaticModel::writeStaticModel(const string &basename,
<< "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
else
{
- sparseHelper(2, col0_output, k, 0, output_type);
- col0_output << "=" << eq + 1 << ";" << endl;
-
- sparseHelper(2, col1_output, k, 1, output_type);
- col1_output << "=" << col_idx_sym + 1 << ";" << endl;
-
- sparseHelper(2, col2_output, k, 2, output_type);
- col2_output << "=";
- sparseHelper(2, col2_output, k-1, 2, output_type);
- col2_output << ";" << endl;
+ i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << "=" << eq + 1 << ";" << endl;
+ j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << "=" << col_idx_sym + 1 << ";" << endl;
+ v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="
+ << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
k++;
}
}
}
if (output_type != ExprNodeOutputType::juliaStaticModel)
- d_output[i] << col0_output.str() << col1_output.str() << col2_output.str();
+ d_output[i] << i_output.str() << j_output.str() << v_output.str();
}
if (output_type == ExprNodeOutputType::matlabStaticModel)
@@ -1407,18 +1414,17 @@ StaticModel::writeStaticModel(const string &basename,
ncols *= JacobianColsNbr;
ntt += temporary_terms_derivatives[i].size();
string gname = "g" + to_string(i);
- string vname = "v" + to_string(i);
string gprevname = "g" + to_string(i-1);
init_output.str("");
end_output.str("");
if (derivatives[i].size())
{
- init_output << vname << " = zeros(" << NNZDerivatives[i] << ",3);";
+ init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl
+ << gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl
+ << gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl;
end_output << gname << " = sparse("
- << vname << "(:,1),"
- << vname << "(:,2),"
- << vname << "(:,3),"
+ << gname << "_i," << gname << "_j," << gname << "_v,"
<< nrows << "," << ncols << ");";
}
else
@@ -1440,13 +1446,19 @@ StaticModel::writeStaticModel(const string &basename,
for (size_t i = 0; i < d_output.size(); i++)
{
string funcname = i == 0 ? "resid" : "g" + to_string(i);
- string argname = i == 0 ? "residual" : i == 1 ? "g1" : "v" + to_string(i);
StaticOutput << "void static_" << funcname << "_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T)" << endl
<< "{" << endl
<< tt_output[i].str()
<< "}" << endl
<< endl
- << "void static_" << funcname << "(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *" << argname << ")" << endl
+ << "void static_" << funcname << "(const double *y, const double *x, int nb_row_x, const double *params, const double *T, ";
+ if (i == 0)
+ StaticOutput << "double *residual";
+ else if (i == 1)
+ StaticOutput << "double *g1";
+ else
+ StaticOutput << "double *" << funcname << "_i, double *" << funcname << "_j, double *" << funcname << "_v";
+ StaticOutput << ")" << endl
<< "{" << endl;
if (i == 0)
StaticOutput << " double lhs, rhs;" << endl;
@@ -1732,10 +1744,21 @@ StaticModel::writeStaticCFile(const string &basename) const
<< endl
<< " if (nlhs >= 3)" << endl
<< " {" << endl
- << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3 << ", mxREAL);" << endl
- << " double *v2 = mxGetPr(plhs[2]);" << endl
+ << " mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+ << " mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+ << " mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
<< " static_g2_tt(y, x, nb_row_x, params, T);" << endl
- << " static_g2(y, x, nb_row_x, params, T, v2);" << endl
+ << " static_g2(y, x, nb_row_x, params, T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
+ << " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
+ << " mxArray *n = mxCreateDoubleScalar(" << symbol_table.endo_nbr()*symbol_table.endo_nbr() << ");" << endl
+ << " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
+ << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
+ << " plhs[2] = plhs_sparse[0];" << endl
+ << " mxDestroyArray(g2_i);" << endl
+ << " mxDestroyArray(g2_j);" << endl
+ << " mxDestroyArray(g2_v);" << endl
+ << " mxDestroyArray(m);" << endl
+ << " mxDestroyArray(n);" << endl
<< " }" << endl
<< endl
<< " free(T);" << endl
--
GitLab