From ad5e196d3064fd3b431b0a0dd34c66f42b52dbc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org> Date: Tue, 23 Jun 2020 15:13:04 +0200 Subject: [PATCH] Block decomposition now compatible with 'use_dll' option --- src/DataTree.cc | 7 + src/DataTree.hh | 4 +- src/DynamicModel.cc | 431 +++++++++++++++++++++++++++++++++++++++----- src/DynamicModel.hh | 7 + src/ModFile.cc | 4 +- src/ModelTree.cc | 8 +- src/ModelTree.hh | 2 +- src/StaticModel.cc | 236 ++++++++++++++++++++++-- src/StaticModel.hh | 6 + 9 files changed, 635 insertions(+), 70 deletions(-) diff --git a/src/DataTree.cc b/src/DataTree.cc index 37d429f5..3e279896 100644 --- a/src/DataTree.cc +++ b/src/DataTree.cc @@ -929,6 +929,13 @@ DataTree::writePowerDeriv(ostream &output) const << "}" << endl; } +void +DataTree::writePowerDerivHeader(ostream &output) const +{ + if (isBinaryOpUsed(BinaryOpcode::powerDeriv)) + output << "double getPowerDeriv(double x, double p, int k);" << endl; +} + string DataTree::packageDir(const string &package) { diff --git a/src/DataTree.hh b/src/DataTree.hh index 23ce9c11..4a82c379 100644 --- a/src/DataTree.hh +++ b/src/DataTree.hh @@ -278,8 +278,10 @@ public: //! Returns the minimum lag (as a negative number) of the given symbol in the whole data tree (and not only in the equations !!) /*! Returns 0 if the symbol is not used */ int minLagForSymbol(int symb_id) const; - //! Write getPowerDeriv in C + //! Write getPowerDeriv in C (function body) void writePowerDeriv(ostream &output) const; + //! Write getPowerDeriv in C (prototype) + void writePowerDerivHeader(ostream &output) const; //! Thrown when trying to access an unknown variable by deriv_id class UnknownDerivIDException { diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 5abc8f43..db45cf67 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -305,7 +305,10 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutpu // Write temporary terms for derivatives write_eq_tt(blocks[blk].size); - output << " if stochastic_mode" << endl; + if (isCOutput(output_type)) + output << " if (stochastic_mode) {" << endl; + else + output << " if stochastic_mode" << endl; ostringstream i_output, j_output, v_output; int line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type); @@ -323,7 +326,7 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutpu v_output << ';' << endl; line_counter++; } - assert(line_counter == nze_stochastic+1); + assert(line_counter == nze_stochastic+ARRAY_SUBSCRIPT_OFFSET(output_type)); output << i_output.str() << j_output.str() << v_output.str(); i_output.str(""); @@ -344,7 +347,7 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutpu v_output << ';' << endl; line_counter++; } - assert(line_counter == nze_exo+1); + assert(line_counter == nze_exo+ARRAY_SUBSCRIPT_OFFSET(output_type)); output << i_output.str() << j_output.str() << v_output.str(); i_output.str(""); @@ -365,7 +368,7 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutpu v_output << ';' << endl; line_counter++; } - assert(line_counter == nze_exo_det+1); + assert(line_counter == nze_exo_det+ARRAY_SUBSCRIPT_OFFSET(output_type)); output << i_output.str() << j_output.str() << v_output.str(); i_output.str(""); @@ -386,14 +389,17 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutpu v_output << ';' << endl; line_counter++; } - assert(line_counter == nze_other_endo+1); + assert(line_counter == nze_other_endo+ARRAY_SUBSCRIPT_OFFSET(output_type)); output << i_output.str() << j_output.str() << v_output.str(); // Deterministic mode if (simulation_type != BlockSimulationType::evaluateForward && simulation_type != BlockSimulationType::evaluateBackward) { - output << " else" << endl; + if (isCOutput(output_type)) + output << " } else {" << endl; + else + output << " else" << endl; i_output.str(""); j_output.str(""); v_output.str(""); @@ -439,10 +445,39 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutpu line_counter++; } } - assert(line_counter == nze_deterministic+1); + assert(line_counter == nze_deterministic+ARRAY_SUBSCRIPT_OFFSET(output_type)); output << i_output.str() << j_output.str() << v_output.str(); } - output << " end" << endl; + if (isCOutput(output_type)) + output << " }" << endl; + else + output << " end" << endl; +} + +int +DynamicModel::nzeDeterministicJacobianForBlock(int blk) const +{ + BlockSimulationType simulation_type = blocks[blk].simulation_type; + int block_recursive_size = blocks[blk].getRecursiveSize(); + + int nze_deterministic = 0; + if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete + || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) + nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(), + [=](const auto &kv) { + auto [eq, var, lag] = kv.first; + return eq >= block_recursive_size && var >= block_recursive_size; + }); + else if (simulation_type == BlockSimulationType::solveBackwardSimple + || simulation_type == BlockSimulationType::solveForwardSimple + || simulation_type == BlockSimulationType::solveBackwardComplete + || simulation_type == BlockSimulationType::solveForwardComplete) + nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(), + [](const auto &kv) { + auto [eq, var, lag] = kv.first; + return lag == 0; + }); + return nze_deterministic; } void @@ -455,27 +490,10 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) 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(); // Number of nonzero derivatives for the various Jacobians int nze_stochastic = blocks_derivatives[blk].size(); - int nze_deterministic = 0; - if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete - || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) - nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(), - [=](const auto &kv) { - auto [eq, var, lag] = kv.first; - return eq >= block_recursive_size && var >= block_recursive_size; - }); - else if (simulation_type == BlockSimulationType::solveBackwardSimple - || simulation_type == BlockSimulationType::solveForwardSimple - || simulation_type == BlockSimulationType::solveBackwardComplete - || simulation_type == BlockSimulationType::solveForwardComplete) - nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(), - [](const auto &kv) { - auto [eq, var, lag] = kv.first; - return lag == 0; - }); + int nze_deterministic = nzeDeterministicJacobianForBlock(blk); int nze_other_endo = blocks_derivatives_other_endo[blk].size(); int nze_exo = blocks_derivatives_exo[blk].size(); int nze_exo_det = blocks_derivatives_exo_det[blk].size(); @@ -483,6 +501,11 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const string filename = packageDir(basename + ".block") + "/dynamic_" + to_string(blk+1) + ".m"; ofstream output; output.open(filename, ios::out | ios::binary); + if (!output.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } output << "%" << endl << "% " << filename << " : Computes dynamic version of one block" << endl @@ -566,6 +589,206 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const } } +void +DynamicModel::writeDynamicPerBlockCFiles(const string &basename) const +{ + temporary_terms_t temporary_terms; // Temp terms written so far + + for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++) + { + BlockSimulationType simulation_type = blocks[blk].simulation_type; + int block_size = blocks[blk].size; + int block_mfs_size = blocks[blk].mfs_size; + + // Number of nonzero derivatives for the various Jacobians + int nze_stochastic = blocks_derivatives[blk].size(); + int nze_deterministic = nzeDeterministicJacobianForBlock(blk); + int nze_other_endo = blocks_derivatives_other_endo[blk].size(); + int nze_exo = blocks_derivatives_exo[blk].size(); + int nze_exo_det = blocks_derivatives_exo_det[blk].size(); + + string filename = basename + "/model/src/dynamic_" + to_string(blk+1) + ".c"; + ofstream output; + output.open(filename, ios::out | ios::binary); + if (!output.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + + output << "/* Block " << blk+1 << endl + << " " << BlockSim(simulation_type) << " */" << endl + << endl + << "#include <math.h>" << endl + << "#include <stdlib.h>" << endl + << "#include <stdbool.h>" << endl + << R"(#include "mex.h")" << endl + << endl + << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl + << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl + << endl; + + // Write function definition if BinaryOpcode::powerDeriv is used + writePowerDerivHeader(output); + + output << endl; + + if (simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + output << "void dynamic_" << blk+1 << "(double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, double *T, int it_, bool stochastic_mode, double *g1_i, double *g1_j, double *g1_v, double *g1_x_i, double *g1_x_j, double *g1_x_v, double *g1_xd_i, double *g1_xd_j, double *g1_xd_v, double *g1_o_i, double *g1_o_j, double *g1_o_v)" << endl; + else + output << "void dynamic_" << blk+1 << "(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, double *T, int it_, bool stochastic_mode, double *residual, double *g1_i, double *g1_j, double *g1_v, double *g1_x_i, double *g1_x_j, double *g1_x_v, double *g1_xd_i, double *g1_xd_j, double *g1_xd_v, double *g1_o_i, double *g1_o_j, double *g1_o_v)" << endl; + output << '{' << endl; + + writeDynamicPerBlockHelper(blk, output, ExprNodeOutputType::CDynamicModel, temporary_terms, + nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo); + + output << '}' << endl + << endl; + + ostringstream header; + if (simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + header << "void dynamic_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, const mxArray *steady_state, mxArray *T, const mxArray *it_, const mxArray *stochastic_mode, mxArray **g1, mxArray **g1_x, mxArray **g1_xd, mxArray **g1_o)"; + else + header << "void dynamic_" << blk+1 << "_mx(const mxArray *y, const mxArray *x, const mxArray *params, const mxArray *steady_state, mxArray *T, const mxArray *it_, const mxArray *stochastic_mode, mxArray **residual, mxArray **g1, mxArray **g1_x, mxArray **g1_xd, mxArray **g1_o)"; + output << header.str() << endl + << '{' << endl + << " int nb_row_x = mxGetM(x);" << endl; + + if (simulation_type != BlockSimulationType::evaluateForward + && simulation_type != BlockSimulationType::evaluateBackward) + output << " *residual = mxCreateDoubleMatrix(" << block_mfs_size << ",1,mxREAL);" << endl; + + output << " mxArray *g1_i = NULL, *g1_j = NULL, *g1_v = NULL;" << endl + << " mxArray *g1_x_i = NULL, *g1_x_j = NULL, *g1_x_v = NULL;" << endl + << " mxArray *g1_xd_i = NULL, *g1_xd_j = NULL, *g1_xd_v = NULL;" << endl + << " mxArray *g1_o_i = NULL, *g1_o_j = NULL, *g1_o_v = NULL;" << endl + << " if (mxGetScalar(stochastic_mode)) {" << endl + << " g1_i=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl + << " g1_j=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl + << " g1_v=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl + << " g1_x_i=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl + << " g1_x_j=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl + << " g1_x_v=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl + << " g1_xd_i=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl + << " g1_xd_j=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl + << " g1_xd_v=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl + << " g1_o_i=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl + << " g1_o_j=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl + << " g1_o_v=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl; + if (simulation_type != BlockSimulationType::evaluateForward + && simulation_type != BlockSimulationType::evaluateBackward) + output << " } else {" << endl + << " g1_i=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl + << " g1_j=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl + << " g1_v=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl; + output << " }" << endl + << endl; + + if (simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + output << " dynamic_" << blk+1 << "(mxGetPr(y), mxGetPr(x), nb_row_x, mxGetPr(params), mxGetPr(steady_state), mxGetPr(T), mxGetScalar(it_), mxGetScalar(stochastic_mode), g1_i ? mxGetPr(g1_i) : NULL, g1_j ? mxGetPr(g1_j) : NULL, g1_v ? mxGetPr(g1_v) : NULL, g1_x_i ? mxGetPr(g1_x_i) : NULL, g1_x_j ? mxGetPr(g1_x_j) : NULL, g1_x_v ? mxGetPr(g1_x_v) : NULL, g1_xd_i ? mxGetPr(g1_xd_i) : NULL, g1_xd_j ? mxGetPr(g1_xd_j) : NULL, g1_xd_v ? mxGetPr(g1_xd_v) : NULL, g1_o_i ? mxGetPr(g1_o_i) : NULL, g1_o_j ? mxGetPr(g1_o_j) : NULL, g1_o_v ? mxGetPr(g1_o_v) : NULL);" << endl; + else + output << " dynamic_" << blk+1 << "(mxGetPr(y), mxGetPr(x), nb_row_x, mxGetPr(params), mxGetPr(steady_state), mxGetPr(T), mxGetScalar(it_), mxGetScalar(stochastic_mode), mxGetPr(*residual), g1_i ? mxGetPr(g1_i) : NULL, g1_j ? mxGetPr(g1_j) : NULL, g1_v ? mxGetPr(g1_v) : NULL, g1_x_i ? mxGetPr(g1_x_i) : NULL, g1_x_j ? mxGetPr(g1_x_j) : NULL, g1_x_v ? mxGetPr(g1_x_v) : NULL, g1_xd_i ? mxGetPr(g1_xd_i) : NULL, g1_xd_j ? mxGetPr(g1_xd_j) : NULL, g1_xd_v ? mxGetPr(g1_xd_v) : NULL, g1_o_i ? mxGetPr(g1_o_i) : NULL, g1_o_j ? mxGetPr(g1_o_j) : NULL, g1_o_v ? mxGetPr(g1_o_v) : NULL);" << endl; + + output << endl + << " if (mxGetScalar(stochastic_mode)) {" << endl + << " mxArray *m = mxCreateDoubleScalar(" << block_size << ");" << endl + << " mxArray *n = mxCreateDoubleScalar(" << blocks_jacob_cols_endo[blk].size() << ");" << endl + << " mxArray *plhs[1];" << endl + << " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl + << R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl + << " *g1=plhs[0];" << endl + << " mxDestroyArray(g1_i);" << endl + << " mxDestroyArray(g1_j);" << endl + << " mxDestroyArray(g1_v);" << endl + << " mxDestroyArray(n);" << endl + << " n = mxCreateDoubleScalar(" << blocks_jacob_cols_exo[blk].size() << ");" << endl + << " mxArray *prhs_x[5] = { g1_x_i, g1_x_j, g1_x_v, m, n };" << endl + << R"( mexCallMATLAB(1, plhs, 5, prhs_x, "sparse");)" << endl + << " *g1_x=plhs[0];" << endl + << " mxDestroyArray(g1_x_i);" << endl + << " mxDestroyArray(g1_x_j);" << endl + << " mxDestroyArray(g1_x_v);" << endl + << " mxDestroyArray(n);" << endl + << " n = mxCreateDoubleScalar(" << blocks_jacob_cols_exo_det[blk].size() << ");" << endl + << " mxArray *prhs_xd[5] = { g1_xd_i, g1_xd_j, g1_xd_v, m, n };" << endl + << R"( mexCallMATLAB(1, plhs, 5, prhs_xd, "sparse");)" << endl + << " *g1_xd=plhs[0];" << endl + << " mxDestroyArray(g1_xd_i);" << endl + << " mxDestroyArray(g1_xd_j);" << endl + << " mxDestroyArray(g1_xd_v);" << endl + << " mxDestroyArray(n);" << endl + << " n = mxCreateDoubleScalar(" << blocks_jacob_cols_other_endo[blk].size() << ");" << endl + << " mxArray *prhs_o[5] = { g1_o_i, g1_o_j, g1_o_v, m, n };" << endl + << R"( mexCallMATLAB(1, plhs, 5, prhs_o, "sparse");)" << endl + << " *g1_o=plhs[0];" << endl + << " mxDestroyArray(g1_o_i);" << endl + << " mxDestroyArray(g1_o_j);" << endl + << " mxDestroyArray(g1_o_v);" << endl + << " mxDestroyArray(n);" << endl + << " mxDestroyArray(m);" << endl + << " } else {" << endl; + switch (simulation_type) + { + case BlockSimulationType::evaluateForward: + case BlockSimulationType::evaluateBackward: + output << " *g1=mxCreateDoubleMatrix(0,0,mxREAL);" << endl; + break; + case BlockSimulationType::solveBackwardSimple: + case BlockSimulationType::solveForwardSimple: + case BlockSimulationType::solveBackwardComplete: + case BlockSimulationType::solveForwardComplete: + output << " mxArray *m = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl + << " mxArray *n = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl + << " mxArray *plhs[1];" << endl + << " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl + << R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl + << " *g1=plhs[0];" << endl + << " mxDestroyArray(g1_i);" << endl + << " mxDestroyArray(g1_j);" << endl + << " mxDestroyArray(g1_v);" << endl + << " mxDestroyArray(n);" << endl + << " mxDestroyArray(m);" << endl; + break; + case BlockSimulationType::solveTwoBoundariesSimple: + case BlockSimulationType::solveTwoBoundariesComplete: + output << " mxArray *m = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl + << " mxArray *n = mxCreateDoubleScalar(" << 3*block_mfs_size << ");" << endl + << " mxArray *plhs[1];" << endl + << " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl + << R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl + << " *g1=plhs[0];" << endl + << " mxDestroyArray(g1_i);" << endl + << " mxDestroyArray(g1_j);" << endl + << " mxDestroyArray(g1_v);" << endl + << " mxDestroyArray(n);" << endl + << " mxDestroyArray(m);" << endl; + break; + default: + break; + } + output << " *g1_x=mxCreateDoubleMatrix(0,0,mxREAL);" << endl + << " *g1_xd=mxCreateDoubleMatrix(0,0,mxREAL);" << endl + << " *g1_o=mxCreateDoubleMatrix(0,0,mxREAL);" << endl + << " }" << endl + << "}" << endl; + output.close(); + + filename = basename + "/model/src/dynamic_" + to_string(blk+1) + ".h"; + ofstream header_output; + header_output.open(filename, ios::out | ios::binary); + if (!header_output.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + header_output << header.str() << ';' << endl; + header_output.close(); + } +} + void DynamicModel::writeDynamicBytecode(const string &basename) const { @@ -574,8 +797,6 @@ DynamicModel::writeDynamicBytecode(const string &basename) const unsigned int instruction_number = 0; bool file_open = false; - filesystem::create_directories(basename + "/model/bytecode"); - string main_name = basename + "/model/bytecode/dynamic.cod"; code_file.open(main_name, ios::out | ios::binary | ios::ate); if (!code_file.is_open()) @@ -837,7 +1058,6 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename, bool linear_deco vector<int> feedback_variables; bool file_open = false; string main_name; - filesystem::create_directories(basename + "/model/bytecode"); if (linear_decomposition) main_name = basename + "/model/bytecode/non_linear.cod"; else @@ -1195,7 +1415,6 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const void DynamicModel::writeDynamicCFile(const string &basename) const { - filesystem::create_directories(basename + "/model/src"); string filename = basename + "/model/src/dynamic.c"; int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size(); @@ -1408,6 +1627,92 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const output.close(); } +void +DynamicModel::writeDynamicBlockCFile(const string &basename) const +{ + string filename = basename + "/model/src/dynamic.c"; + + ofstream output; + output.open(filename, ios::out | ios::binary); + if (!output.is_open()) + { + cerr << "Error: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + + output << "#include <math.h>" << endl + << R"(#include "mex.h")" << endl; + + for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++) + output << R"(#include "dynamic_)" << blk+1 << R"(.h")" << endl; + + output << endl; + writePowerDeriv(output); + + output << endl + << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl + << "{" << endl + << " if (nrhs != 8)" << endl + << R"( mexErrMsgTxt("Requires exactly 8 input arguments");)" << endl + << " if (nlhs > 7)" << endl + << R"( mexErrMsgTxt("Accepts at most 7 output arguments");)" << endl + << " int nblock = (int) mxGetScalar(prhs[0]);" << endl + << " const mxArray *y = prhs[1], *x = prhs[2], *params = prhs[3], *steady_state = prhs[4], *T = prhs[5], *it_ = prhs[6], *stochastic_mode = prhs[7];" << endl + << " mxArray *T_new = mxDuplicateArray(T);" << endl + << " mxArray *y_new = mxDuplicateArray(y);" << endl + << " mxArray *residual, *g1, *g1_x, *g1_xd, *g1_o;" << endl + << " switch (nblock)" << endl + << " {" << endl; + + for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++) + { + output << " case " << blk+1 << ':' << endl; + + BlockSimulationType simulation_type = blocks[blk].simulation_type; + + if (simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + output << " dynamic_" << blk+1 << "_mx(y_new, x, params, steady_state, T_new, it_, stochastic_mode, &g1, &g1_x, &g1_xd, &g1_o);" << endl + << " residual = mxCreateDoubleMatrix(0,0,mxREAL);" << endl; + else + output << " dynamic_" << blk+1 << "_mx(y, x, params, steady_state, T_new, it_, stochastic_mode, &residual, &g1, &g1_x, &g1_xd, &g1_o);" << endl; + output << " break;" << endl; + } + output << " }" << endl + << endl + << " if (nlhs >= 1)" << endl + << " plhs[0] = residual;" << endl + << " else" << endl + << " mxDestroyArray(residual);" << endl + << " if (nlhs >= 2)" << endl + << " plhs[1] = y_new;" << endl + << " else" << endl + << " mxDestroyArray(y_new);" << endl + << " if (nlhs >= 3)" << endl + << " plhs[2] = T_new;" << endl + << " else" << endl + << " mxDestroyArray(T_new);" << endl + << " if (nlhs >= 4)" << endl + << " plhs[3] = g1;" << endl + << " else" << endl + << " mxDestroyArray(g1);" << endl + << " if (nlhs >= 5)" << endl + << " plhs[4] = g1_x;" << endl + << " else" << endl + << " mxDestroyArray(g1_x);" << endl + << " if (nlhs >= 6)" << endl + << " plhs[5] = g1_xd;" << endl + << " else" << endl + << " mxDestroyArray(g1_xd);" << endl + << " if (nlhs >= 7)" << endl + << " plhs[6] = g1_o;" << endl + << " else" << endl + << " mxDestroyArray(g1_o);" << endl + << "}" << endl; + + output.close(); +} + void DynamicModel::writeWrapperFunctions(const string &basename, const string &ending) const { @@ -4256,24 +4561,64 @@ DynamicModel::computeBlockDynJacobianCols() void DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const { - if ((block && bytecode) || linear_decomposition) - writeDynamicBlockBytecode(basename, linear_decomposition); - else if (!block && bytecode) - writeDynamicBytecode(basename); - else if (block && !bytecode) + filesystem::path model_dir{basename}; + model_dir /= "model"; + if (use_dll) + filesystem::create_directories(model_dir / "src"); + if (bytecode) + filesystem::create_directories(model_dir / "bytecode"); + + if (linear_decomposition) { - writeDynamicPerBlockMFiles(basename); - writeDynamicBlockMFile(basename); + if (bytecode) + writeDynamicBlockBytecode(basename, linear_decomposition); + else + { + cerr << "'linear_decomposition' option requires the 'bytecode' option" << endl; + exit(EXIT_FAILURE); + } } - else if (use_dll) + else if (block) { - writeDynamicCFile(basename); - compileMEX(basename, "dynamic", mexext, matlabroot, dynareroot); + if (bytecode) + writeDynamicBlockBytecode(basename, linear_decomposition); + else if (use_dll) + { + writeDynamicPerBlockCFiles(basename); + writeDynamicBlockCFile(basename); + vector<filesystem::path> src_files{blocks.size() + 1}; + for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++) + src_files[blk] = model_dir / "src" / ("dynamic_" + to_string(blk+1) + ".c"); + src_files[blocks.size()] = model_dir / "src" / "dynamic.c"; + compileMEX(basename, "dynamic", mexext, src_files, matlabroot, dynareroot); + } + else if (julia) + { + cerr << "'block' option is not available with Julia" << endl; + exit(EXIT_FAILURE); + } + else // M-files + { + writeDynamicPerBlockMFiles(basename); + writeDynamicBlockMFile(basename); + } } - else if (julia) - writeDynamicJuliaFile(basename); else - writeDynamicMFile(basename); + { + if (bytecode) + writeDynamicBytecode(basename); + else if (use_dll) + { + writeDynamicCFile(basename); + compileMEX(basename, "dynamic", mexext, { model_dir / "src" / "dynamic.c" }, + matlabroot, dynareroot); + } + else if (julia) + writeDynamicJuliaFile(basename); + else + writeDynamicMFile(basename); + } + writeSetAuxiliaryVariables(basename, julia); } diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index 42548b74..17064c1c 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -133,10 +133,17 @@ private: void writeDynamicModel(const string &basename, ostream &DynamicOutput, bool use_dll, bool julia) const; //! Writes the main dynamic function of block decomposed model (MATLAB version) void writeDynamicBlockMFile(const string &basename) const; + //! Writes the main dynamic function of block decomposed model (C version) + void writeDynamicBlockCFile(const string &basename) const; + /* 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 void writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, 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) void writeDynamicPerBlockMFiles(const string &basename) const; + //! Writes the per-block dynamic files of block decomposed model (C version) + void writeDynamicPerBlockCFiles(const string &basename) const; //! Writes the code of the block-decomposed model in virtual machine bytecode void writeDynamicBlockBytecode(const string &basename, bool linear_decomposition) const; //! Writes the code of the model in virtual machine bytecode diff --git a/src/ModFile.cc b/src/ModFile.cc index eeff5fc4..7670e2da 100644 --- a/src/ModFile.cc +++ b/src/ModFile.cc @@ -187,9 +187,9 @@ ModFile::checkPass(bool nostrict, bool stochastic) exit(EXIT_FAILURE); } - if (use_dll && (block || bytecode || linear_decomposition)) + if (use_dll && bytecode) { - cerr << "ERROR: In 'model' block, 'use_dll' option is not compatible with 'block', 'bytecode' or 'linear_decomposition'" << endl; + cerr << "ERROR: In 'model' block, 'use_dll' option is not compatible with 'bytecode'" << endl; exit(EXIT_FAILURE); } if (block && linear_decomposition) diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 4180f5db..0774c321 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -1977,7 +1977,7 @@ ModelTree::matlab_arch(const string &mexext) } void -ModelTree::compileMEX(const string &basename, const string &funcname, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const +ModelTree::compileMEX(const string &basename, const string &funcname, const string &mexext, const vector<filesystem::path> &src_files, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const { const string opt_flags = "-O3 -g0 --param ira-max-conflict-table-size=1 -fno-forward-propagate -fno-gcse -fno-dce -fno-dse -fno-tree-fre -fno-tree-pre -fno-tree-cselim -fno-tree-dse -fno-tree-dce -fno-tree-pta -fno-gcse-after-reload"; @@ -2056,8 +2056,6 @@ ModelTree::compileMEX(const string &basename, const string &funcname, const stri } } - auto model_dir = filesystem::path{basename} / "model" / "src"; - filesystem::path src{model_dir / (funcname + ".c")}; filesystem::path mex_dir{"+" + basename}; filesystem::path binary{mex_dir / (funcname + "." + mexext)}; @@ -2089,7 +2087,9 @@ ModelTree::compileMEX(const string &basename, const string &funcname, const stri if (!user_set_add_flags.empty()) cmd << user_set_add_flags << " "; - cmd << src << " -o " << binary << " "; + for (auto &src : src_files) + cmd << src << " "; + cmd << "-o " << binary << " "; if (user_set_subst_libs.empty()) cmd << libs; diff --git a/src/ModelTree.hh b/src/ModelTree.hh index 1b3eb22d..518d3bdb 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -408,7 +408,7 @@ private: //! Returns the name of the MATLAB architecture given the extension used for MEX files static string matlab_arch(const string &mexext); //! Compiles a MEX file - void compileMEX(const string &basename, const string &funcname, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const; + void compileMEX(const string &basename, const string &funcname, const string &mexext, const vector<filesystem::path> &src_files, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const; public: ModelTree(SymbolTable &symbol_table_arg, diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 87e09b09..1fad42f9 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -240,6 +240,11 @@ StaticModel::writeStaticPerBlockMFiles(const string &basename) const string filename = packageDir(basename + ".block") + "/static_" + to_string(blk+1) + ".m"; ofstream output; output.open(filename, ios::out | ios::binary); + if (!output.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } output << "%" << endl << "% " << filename << " : Computes static version of one block" << endl << "%" << endl @@ -279,6 +284,100 @@ StaticModel::writeStaticPerBlockMFiles(const string &basename) const } } +void +StaticModel::writeStaticPerBlockCFiles(const string &basename) const +{ + temporary_terms_t temporary_terms; // Temp terms written so far + + for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++) + { + BlockSimulationType simulation_type = blocks[blk].simulation_type; + + string filename = basename + "/model/src/static_" + to_string(blk+1) + ".c"; + ofstream output; + output.open(filename, ios::out | ios::binary); + if (!output.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + output << "/* Block " << blk+1 << endl + << " " << BlockSim(simulation_type) << " */" << endl + << endl + << "#include <math.h>" << endl + << "#include <stdlib.h>" << endl + << R"(#include "mex.h")" << endl + << endl + << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl + << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl + << endl; + + // Write function definition if BinaryOpcode::powerDeriv is used + writePowerDerivHeader(output); + + output << endl; + + if (simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + output << "void static_" << blk+1 << "(double *y, const double *x, const double *params, double *T)" << endl; + else + output << "void static_" << blk+1 << "(const double *y, const double *x, const double *params, double *T, double *residual, double *g1_i, double *g1_j, double *g1_v)" << endl; + output << '{' << endl; + + writeStaticPerBlockHelper(blk, output, ExprNodeOutputType::CStaticModel, temporary_terms); + + output << '}' << endl + << endl; + + ostringstream header; + if (simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + { + header << "void static_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, mxArray *T)"; + output << header.str() << endl + << '{' << endl + << " static_" << blk+1 << "(mxGetPr(y), mxGetPr(x), mxGetPr(params), mxGetPr(T));" << endl + << '}' << endl; + } + else + { + header << "void static_" << blk+1 << "_mx(const mxArray *y, const mxArray *x, const mxArray *params, mxArray *T, mxArray **residual, mxArray **g1)"; + output << header.str() << endl + << '{' << endl + << " *residual = mxCreateDoubleMatrix(" << blocks[blk].mfs_size << ",1,mxREAL);" << endl + << " mxArray *g1_i = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl + << " mxArray *g1_j = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl + << " mxArray *g1_v = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl + << " static_" << blk+1 << "(mxGetPr(y), mxGetPr(x), mxGetPr(params), mxGetPr(T), mxGetPr(*residual), mxGetPr(g1_i), mxGetPr(g1_j), mxGetPr(g1_v));" << endl + << " mxArray *plhs[1];" << endl + << " mxArray *m = mxCreateDoubleScalar(" << blocks[blk].mfs_size << ");" << endl + << " mxArray *n = mxCreateDoubleScalar(" << blocks[blk].mfs_size << ");" << endl + << " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl + << R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl + << " *g1 = plhs[0];" << endl + << " mxDestroyArray(g1_i);" << endl + << " mxDestroyArray(g1_j);" << endl + << " mxDestroyArray(g1_v);" << endl + << " mxDestroyArray(m);" << endl + << " mxDestroyArray(n);" << endl + << '}' << endl; + } + + output.close(); + + filename = basename + "/model/src/static_" + to_string(blk+1) + ".h"; + ofstream header_output; + header_output.open(filename, ios::out | ios::binary); + if (!header_output.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + header_output << header.str() << ';' << endl; + header_output.close(); + } +} + void StaticModel::writeStaticBytecode(const string &basename) const { @@ -287,8 +386,6 @@ StaticModel::writeStaticBytecode(const string &basename) const unsigned int instruction_number = 0; bool file_open = false; - filesystem::create_directories(basename + "/model/bytecode"); - string main_name = basename + "/model/bytecode/static.cod"; code_file.open(main_name, ios::out | ios::binary | ios::ate); if (!code_file.is_open()) @@ -462,8 +559,6 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const vector<int> feedback_variables; bool file_open = false; - filesystem::create_directories(basename + "/model/bytecode"); - string main_name = basename + "/model/bytecode/static.cod"; code_file.open(main_name, ios::out | ios::binary | ios::ate); if (!code_file.is_open()) @@ -1578,7 +1673,6 @@ void StaticModel::writeStaticCFile(const string &basename) const { // Writing comments and function definition command - filesystem::create_directories(basename + "/model/src"); string filename = basename + "/model/src/static.c"; int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size(); @@ -1661,24 +1755,54 @@ StaticModel::writeStaticJuliaFile(const string &basename) const void StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const { - if (block && bytecode) - writeStaticBlockBytecode(basename); - else if (!block && bytecode) - writeStaticBytecode(basename); - else if (block && !bytecode) + filesystem::path model_dir{basename}; + model_dir /= "model"; + if (use_dll) + filesystem::create_directories(model_dir / "src"); + if (bytecode) + filesystem::create_directories(model_dir / "bytecode"); + + if (block) { - writeStaticPerBlockMFiles(basename); - writeStaticBlockMFile(basename); + if (bytecode) + writeStaticBlockBytecode(basename); + else if (use_dll) + { + writeStaticPerBlockCFiles(basename); + writeStaticBlockCFile(basename); + vector<filesystem::path> src_files{blocks.size() + 1}; + for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++) + src_files[blk] = model_dir / "src" / ("static_" + to_string(blk+1) + ".c"); + src_files[blocks.size()] = model_dir / "src" / "static.c"; + compileMEX(basename, "static", mexext, src_files, matlabroot, dynareroot); + } + else if (julia) + { + cerr << "'block' option is not available with Julia" << endl; + exit(EXIT_FAILURE); + } + else // M-files + { + writeStaticPerBlockMFiles(basename); + writeStaticBlockMFile(basename); + } } - else if (use_dll) + else { - writeStaticCFile(basename); - compileMEX(basename, "static", mexext, matlabroot, dynareroot); + if (bytecode) + writeStaticBytecode(basename); + else if (use_dll) + { + writeStaticCFile(basename); + compileMEX(basename, "static", mexext, { model_dir / "src" / "static.c" }, + matlabroot, dynareroot); + } + else if (julia) + writeStaticJuliaFile(basename); + else // M-files + writeStaticMFile(basename); } - else if (julia) - writeStaticJuliaFile(basename); - else - writeStaticMFile(basename); + writeSetAuxiliaryVariables(basename, julia); } @@ -1727,6 +1851,80 @@ StaticModel::writeStaticBlockMFile(const string &basename) const output.close(); } +void +StaticModel::writeStaticBlockCFile(const string &basename) const +{ + string filename = basename + "/model/src/static.c"; + + ofstream output; + output.open(filename, ios::out | ios::binary); + if (!output.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + + output << "#include <math.h>" << endl + << R"(#include "mex.h")" << endl; + + for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++) + output << R"(#include "static_)" << blk+1 << R"(.h")" << endl; + + output << endl; + writePowerDeriv(output); + + output << endl + << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl + << "{" << endl + << " if (nrhs != 5)" << endl + << R"( mexErrMsgTxt("Requires exactly 5 input arguments");)" << endl + << " if (nlhs > 4)" << endl + << R"( mexErrMsgTxt("Accepts at most 4 output arguments");)" << endl + << " int nblock = (int) mxGetScalar(prhs[0]);" << endl + << " const mxArray *y = prhs[1], *x = prhs[2], *params = prhs[3], *T = prhs[4];" << endl + << " mxArray *T_new = mxDuplicateArray(T);" << endl + << " mxArray *y_new = mxDuplicateArray(y);" << endl + << " mxArray *residual, *g1;" << endl + << " switch (nblock)" << endl + << " {" << endl; + + for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++) + { + output << " case " << blk+1 << ':' << endl; + + BlockSimulationType simulation_type = blocks[blk].simulation_type; + + if (simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + output << " static_" << blk+1 << "_mx(y_new, x, params, T_new);" << endl + << " residual = mxCreateDoubleMatrix(0,0,mxREAL);" << endl + << " g1 = mxCreateDoubleMatrix(0,0,mxREAL);" << endl; + else + output << " static_" << blk+1 << "_mx(y, x, params, T_new, &residual, &g1);" << endl; + output << " break;" << endl; + } + output << " }" << endl + << endl + << " if (nlhs >= 1)" << endl + << " plhs[0] = residual;" << endl + << " else" << endl + << " mxDestroyArray(residual);" << endl + << " if (nlhs >= 2)" << endl + << " plhs[1] = y_new;" << endl + << " else" << endl + << " mxDestroyArray(y_new);" << endl + << " if (nlhs >= 3)" << endl + << " plhs[2] = T_new;" << endl + << " else" << endl + << " mxDestroyArray(T_new);" << endl + << " if (nlhs >= 4)" << endl + << " plhs[3] = g1;" << endl + << " else" << endl + << " mxDestroyArray(g1);" << endl + << "}" << endl; + output.close(); +} + void StaticModel::writeDriverOutput(ostream &output, bool block) const { diff --git a/src/StaticModel.hh b/src/StaticModel.hh index 12327ca7..3ca11f69 100644 --- a/src/StaticModel.hh +++ b/src/StaticModel.hh @@ -48,12 +48,18 @@ private: //! Writes the main static function of block decomposed model (MATLAB version) void writeStaticBlockMFile(const string &basename) const; + //! Writes the main static function of block decomposed model (C version) + void writeStaticBlockCFile(const string &basename) const; + //! Helper for writing a per-block static file of block decomposed model void writeStaticPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms) const; //! Writes the per-block static files of block decomposed model (MATLAB version) void writeStaticPerBlockMFiles(const string &basename) const; + //! Writes the per-block static files of block decomposed model (C version) + void writeStaticPerBlockCFiles(const string &basename) const; + //! Writes the code of the block-decomposed model in virtual machine bytecode void writeStaticBlockBytecode(const string &basename) const; -- GitLab