diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index b7dbe15b..fe4d40b1 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -1538,6 +1538,8 @@ DynamicModel::writeDynamicCFile(const string &basename, const int order) const string filename_mex = basename + "/model/src/dynamic_mex.c"; ofstream mDynamicModelFile, mDynamicMexFile; + int ntt = temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size(); + mDynamicModelFile.open(filename, ios::out | ios::binary); if (!mDynamicModelFile.is_open()) { @@ -1570,9 +1572,13 @@ DynamicModel::writeDynamicCFile(const string &basename, const int order) const writePowerDerivCHeader(mDynamicModelFile); writeNormcdfCHeader(mDynamicModelFile); + mDynamicModelFile << endl; + // Writing the function body writeDynamicModel(mDynamicModelFile, true, false); + mDynamicModelFile << endl; + writePowerDeriv(mDynamicModelFile); writeNormcdf(mDynamicModelFile); mDynamicModelFile.close(); @@ -1592,73 +1598,83 @@ DynamicModel::writeDynamicCFile(const string &basename, const int order) const << " * Warning : this file is generated automatically by Dynare" << endl << " * from model file (.mod)" << endl << endl - << " */" << endl << endl - << "#include \"mex.h\"" << endl << endl - << "void Dynamic(double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, double *g1, double *v2, double *v3);" << endl + << " */" << endl + << endl + << "#include " << endl + << "#include \"mex.h\"" << endl + << endl + << "const int ntt = " << ntt << ";" << endl + << "void dynamic_resid_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T);" << endl + << "void dynamic_resid(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *residual);" << endl + << "void dynamic_g1_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T);" << endl + << "void dynamic_g1(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *g1);" << endl + << "void dynamic_g2_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T);" << endl + << "void dynamic_g2(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *v2);" << endl + << "void dynamic_g3_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T);" << endl + << "void dynamic_g3(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *v3);" << endl << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl << "{" << endl - << " double *y, *x, *params, *steady_state;" << endl - << " double *residual, *g1, *v2, *v3;" << endl - << " int nb_row_x, it_;" << endl - << endl << " /* Check that no derivatives of higher order than computed are being requested */" << endl << " if (nlhs > " << order + 1 << ")" << endl << " mexErrMsgTxt(\"Derivatives of higher order than computed have been requested\");" << endl << " /* Create a pointer to the input matrix y. */" << endl - << " y = mxGetPr(prhs[0]);" << endl + << " double *y = mxGetPr(prhs[0]);" << endl << endl << " /* Create a pointer to the input matrix x. */" << endl - << " x = mxGetPr(prhs[1]);" << endl + << " double *x = mxGetPr(prhs[1]);" << endl << endl << " /* Create a pointer to the input matrix params. */" << endl - << " params = mxGetPr(prhs[2]);" << endl + << " double *params = mxGetPr(prhs[2]);" << endl << endl << " /* Create a pointer to the input matrix steady_state. */" << endl - << " steady_state = mxGetPr(prhs[3]);" << endl + << " double *steady_state = mxGetPr(prhs[3]);" << endl << endl << " /* Fetch time index */" << endl - << " it_ = (int) mxGetScalar(prhs[4]) - 1;" << endl + << " int it_ = (int) mxGetScalar(prhs[4]) - 1;" << endl << endl << " /* Gets number of rows of matrix x. */" << endl - << " nb_row_x = mxGetM(prhs[1]);" << endl + << " int nb_row_x = mxGetM(prhs[1]);" << endl + << endl + << " double *T = (double *) malloc(sizeof(double)*ntt);" << endl - << " residual = NULL;" << endl << " if (nlhs >= 1)" << endl << " {" << endl << " /* Set the output pointer to the output matrix residual. */" << endl << " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl - << " /* Create a C pointer to a copy of the output matrix residual. */" << endl - << " residual = mxGetPr(plhs[0]);" << endl + << " double *residual = mxGetPr(plhs[0]);" << endl + << " dynamic_resid_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl + << " dynamic_resid(y, x, nb_row_x, params, steady_state, it_, T, residual);" << endl << " }" << endl << endl - << " g1 = NULL;" << endl << " if (nlhs >= 2)" << endl << " {" << endl << " /* Set the output pointer to the output matrix g1. */" << endl << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr << ", mxREAL);" << endl - << " /* Create a C pointer to a copy of the output matrix g1. */" << endl - << " g1 = mxGetPr(plhs[1]);" << endl + << " double *g1 = mxGetPr(plhs[1]);" << endl + << " dynamic_g1_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl + << " dynamic_g1(y, x, nb_row_x, params, steady_state, it_, T, g1);" << endl << " }" << endl << endl - << " v2 = NULL;" << endl << " if (nlhs >= 3)" << endl << " {" << endl << " /* Set the output pointer to the output matrix v2. */" << endl << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3 << ", mxREAL);" << endl - << " v2 = mxGetPr(plhs[2]);" << endl + << " double *v2 = mxGetPr(plhs[2]);" << 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 << " }" << endl << endl - << " v3 = NULL;" << endl << " if (nlhs >= 4)" << endl << " {" << endl << " /* Set the output pointer to the output matrix v3. */" << endl << " plhs[3] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3 << ", mxREAL);" << endl - << " v3 = mxGetPr(plhs[3]);" << endl + << " double *v3 = mxGetPr(plhs[3]);" << 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 << " }" << endl << endl - << " /* Call the C subroutines. */" << endl - << " Dynamic(y, x, nb_row_x, params, steady_state, it_, residual, g1, v2, v3);" << endl + << " free(T);" << "}" << endl; mDynamicMexFile.close(); } @@ -2579,40 +2595,46 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput, } else if (output_type == ExprNodeOutputType::CDynamicModel) { - DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, double *g1, double *v2, double *v3)" << endl + DynamicOutput << "void dynamic_resid_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T)" << endl + << "{" << endl + << model_tt_output.str() + << "}" << endl + << endl + << "void dynamic_resid(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *residual)" << endl << "{" << endl << " double lhs, rhs;" << endl - << endl - << " /* Residual equations */" << endl - << model_tt_output.str() << model_output.str() - << " /* Jacobian */" << endl - << " if (g1 == NULL)" << endl - << " return;" << endl + << "}" << endl << endl + << "void dynamic_g1_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T)" << endl + << "{" << endl << jacobian_tt_output.str() + << "}" << endl + << endl + << "void dynamic_g1(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *g1)" << endl + << "{" << endl << jacobian_output.str() - << endl; - - if (second_derivatives.size()) - DynamicOutput << " /* Hessian for endogenous and exogenous variables */" << endl - << " if (v2 == NULL)" << endl - << " return;" << endl - << endl - << hessian_tt_output.str() - << hessian_output.str() - << endl; - - if (third_derivatives.size()) - DynamicOutput << " /* Third derivatives for endogenous and exogenous variables */" << endl - << " if (v3 == NULL)" << endl - << " return;" << endl - << endl - << third_derivatives_tt_output.str() - << third_derivatives_output.str() - << endl; - - DynamicOutput << "}" << endl << endl; + << "}" << endl + << endl + << "void dynamic_g2_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T)" << endl + << "{" << endl + << hessian_tt_output.str() + << "}" << endl + << endl + << "void dynamic_g2(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *v2)" << endl + << "{" << endl + << hessian_output.str() + << "}" << endl + << endl + << "void dynamic_g3_tt(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T)" << endl + << "{" << endl + << third_derivatives_tt_output.str() + << "}" << endl + << endl + << "void dynamic_g3(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *v3)" << endl + << "{" << endl + << third_derivatives_output.str() + << "}" << endl; } else { @@ -6021,7 +6043,7 @@ DynamicModel::dynamicOnlyEquationsNbr() const } bool -DynamicModel::isChecksumMatching(const string &basename) const +DynamicModel::isChecksumMatching(const string &basename, bool block) const { boost::crc_32_type result; @@ -6033,7 +6055,7 @@ DynamicModel::isChecksumMatching(const string &basename) const << equation_tag.second.first << equation_tag.second.second; - ExprNodeOutputType buffer_type = ExprNodeOutputType::CDynamicModel; + ExprNodeOutputType buffer_type = block ? ExprNodeOutputType::matlabDynamicModelSparse : ExprNodeOutputType::CDynamicModel; for (int eq = 0; eq < (int) equations.size(); eq++) { diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index 1e4094a0..b021b969 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -642,7 +642,7 @@ public: //! Returns true if a parameter was used in the model block with a lead or lag bool ParamUsedWithLeadLag() const; - bool isChecksumMatching(const string &basename) const; + bool isChecksumMatching(const string &basename, bool block) const; }; //! Classes to re-order derivatives for various sparse storage formats diff --git a/src/ExprNode.cc b/src/ExprNode.cc index d26f2ff1..8bc9d961 100644 --- a/src/ExprNode.cc +++ b/src/ExprNode.cc @@ -107,7 +107,7 @@ ExprNode::checkIfTemporaryTermThenWrite(ostream &output, ExprNodeOutputType outp if (output_type == ExprNodeOutputType::matlabDynamicModelSparse) output << "T" << idx << "(it_)"; else - if (output_type == ExprNodeOutputType::matlabStaticModelSparse || isCOutput(output_type)) + if (output_type == ExprNodeOutputType::matlabStaticModelSparse) output << "T" << idx; else { diff --git a/src/ModFile.cc b/src/ModFile.cc index ad7da6a2..3908276c 100644 --- a/src/ModFile.cc +++ b/src/ModFile.cc @@ -796,7 +796,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo , const bool nopreprocessoroutput ) const { - bool hasModelChanged = !dynamic_model.isChecksumMatching(basename); + bool hasModelChanged = !dynamic_model.isChecksumMatching(basename, block); if (!check_model_changes) hasModelChanged = true; diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 1dde43ac..a5a67e3c 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -1325,9 +1325,7 @@ ModelTree::writeModelLocalVariableTemporaryTerms(const temporary_terms_t &tto, c temporary_terms_t tt2; for (auto it : tt) { - if (isCOutput(output_type)) - output << "double "; - else if (isJuliaOutput(output_type)) + if (isJuliaOutput(output_type)) output << " @inbounds const "; it.first->writeOutput(output, output_type, tto, temporary_terms_idxs, tef_terms); @@ -1357,9 +1355,7 @@ ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, if (dynamic_cast(*it) != nullptr) (*it)->writeExternalFunctionOutput(output, output_type, tt2, tt_idxs, tef_terms); - if (isCOutput(output_type)) - output << "double "; - else if (isJuliaOutput(output_type)) + if (isJuliaOutput(output_type)) output << " @inbounds "; (*it)->writeOutput(output, output_type, tt, tt_idxs, tef_terms); diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 788f53e4..fd47cb75 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -1615,37 +1615,36 @@ StaticModel::writeStaticModel(const string &basename, } else if (output_type == ExprNodeOutputType::CStaticModel) { - StaticOutput << "void Static(double *y, double *x, int nb_row_x, double *params, double *residual, double *g1, double *v2)" << endl + StaticOutput << "void static_resid_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T)" << endl + << "{" << endl + << model_tt_output.str() + << "}" << endl + << endl + << "void static_resid(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *residual)" << endl << "{" << endl << " double lhs, rhs;" << endl - << endl - << " /* Residual equations */" << endl - << model_tt_output.str() << model_output.str() - << " /* Jacobian */" << endl - << " if (g1 == NULL)" << endl - << " return;" << endl + << "}" << endl << endl + << "void static_g1_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T)" << endl + << "{" << endl << jacobian_tt_output.str() + << "}" << endl + << endl + << "void static_g1(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *g1)" << endl + << "{" << endl << jacobian_output.str() - << endl; - - if (second_derivatives.size()) - StaticOutput << " /* Hessian for endogenous and exogenous variables */" << endl - << " if (v2 == NULL)" << endl - << " return;" << endl - << endl - << hessian_tt_output.str() - << hessian_output.str() - << endl; - if (third_derivatives.size()) - StaticOutput << " /* Third derivatives for endogenous and exogenous variables */" << endl - << " if (v3 == NULL)" << endl - << " return;" << endl - << endl - << third_derivatives_tt_output.str() - << third_derivatives_output.str() - << endl; + << "}" << endl + << endl + << "void static_g2_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T)" << endl + << "{" << endl + << hessian_tt_output.str() + << "}" << endl + << endl + << "void static_g2(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *v2)" << endl + << "{" << endl + << hessian_output.str() + << "}" << endl; } else { @@ -1866,6 +1865,8 @@ StaticModel::writeStaticCFile(const string &basename) const string filename = basename + "/model/src/static.c"; string filename_mex = basename + "/model/src/static_mex.c"; + int ntt = temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size(); + ofstream output; output.open(filename, ios::out | ios::binary); if (!output.is_open()) @@ -1900,9 +1901,12 @@ StaticModel::writeStaticCFile(const string &basename) const writePowerDerivCHeader(output); writeNormcdfCHeader(output); + output << endl; + // Writing the function body writeStaticModel(output, true, false); - output << "}" << endl << endl; + + output << endl; writePowerDeriv(output); writeNormcdf(output); @@ -1922,57 +1926,64 @@ StaticModel::writeStaticCFile(const string &basename) const << " *" << endl << " * Warning : this file is generated automatically by Dynare" << endl << " * from model file (.mod)" << endl << endl - << " */" << endl << endl - << "#include \"mex.h\"" << endl << endl - << "void Static(double *y, double *x, int nb_row_x, double *params, double *residual, double *g1, double *v2);" << endl + << " */" << endl + << endl + << "#include " << endl + << "#include \"mex.h\"" << endl + << endl + << "const int ntt = " << ntt << ";" << endl + << "void static_resid_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T);" << endl + << "void static_resid(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *residual);" << endl + << "void static_g1_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T);" << endl + << "void static_g1(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *g1);" << endl + << "void static_g2_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T);" << endl + << "void static_g2(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *v2);" << endl << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl << "{" << endl - << " double *y, *x, *params;" << endl - << " double *residual, *g1, *v2;" << endl - << " int nb_row_x;" << endl - << endl << " /* Create a pointer to the input matrix y. */" << endl - << " y = mxGetPr(prhs[0]);" << endl + << " double *y = mxGetPr(prhs[0]);" << endl << endl << " /* Create a pointer to the input matrix x. */" << endl - << " x = mxGetPr(prhs[1]);" << endl + << " double *x = mxGetPr(prhs[1]);" << endl << endl << " /* Create a pointer to the input matrix params. */" << endl - << " params = mxGetPr(prhs[2]);" << endl + << " double *params = mxGetPr(prhs[2]);" << endl << endl << " /* Gets number of rows of matrix x. */" << endl - << " nb_row_x = mxGetM(prhs[1]);" << endl + << " int nb_row_x = mxGetM(prhs[1]);" << endl + << endl + << " double *T = (double *) malloc(sizeof(double)*ntt);" << endl - << " residual = NULL;" << endl << " if (nlhs >= 1)" << endl << " {" << endl << " /* Set the output pointer to the output matrix residual. */" << endl << " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl - << " /* Create a C pointer to a copy of the output matrix residual. */" << endl - << " residual = mxGetPr(plhs[0]);" << endl + << " double *residual = mxGetPr(plhs[0]);" << endl + << " static_resid_tt(y, x, nb_row_x, params, T);" << endl + << " static_resid(y, x, nb_row_x, params, T, residual);" << endl << " }" << endl << endl - << " g1 = NULL;" << endl << " if (nlhs >= 2)" << endl << " {" << endl << " /* Set the output pointer to the output matrix g1. */" << endl << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl - << " /* Create a C pointer to a copy of the output matrix g1. */" << endl - << " g1 = mxGetPr(plhs[1]);" << endl + << " double *g1 = mxGetPr(plhs[1]);" << endl + << " static_g1_tt(y, x, nb_row_x, params, T);" << endl + << " static_g1(y, x, nb_row_x, params, T, g1);" << endl << " }" << endl << endl - << " v2 = NULL;" << endl << " if (nlhs >= 3)" << endl << " {" << endl << " /* Set the output pointer to the output matrix v2. */" << endl << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3 << ", mxREAL);" << endl - << " v2 = mxGetPr(plhs[2]);" << endl + << " double *v2 = mxGetPr(plhs[2]);" << endl + << " static_g2_tt(y, x, nb_row_x, params, T);" << endl + << " static_g2(y, x, nb_row_x, params, T, v2);" << endl << " }" << endl << endl - << " /* Call the C subroutines. */" << endl - << " Static(y, x, nb_row_x, params, residual, g1, v2);" << endl - << "}" << endl << endl; + << " free(T);" << endl + << "}" << endl; output.close(); }