stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Bytecode: various simplifications and modernizations

time-shift
Sébastien Villemot 2021-02-03 18:10:01 +01:00
parent 8e9085eea4
commit 4893f0e82c
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
11 changed files with 830 additions and 1465 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

124
mex/sources/bytecode/ErrorHandling.hh
View File

@ -38,11 +38,9 @@
# define CHAR_LENGTH 2
#endif
//#define DEBUG
using namespace std;
const int NO_ERROR_ON_EXIT = 0;
const int ERROR_ON_EXIT = 1;
constexpr int NO_ERROR_ON_EXIT = 0, ERROR_ON_EXIT = 1;
using code_liste_type = vector<pair<Tags, void *>>;
using it_code_type = code_liste_type::const_iterator;
@ -60,7 +58,7 @@ public:
return ErrorMsg;
}
inline void
completeErrorMsg(string ErrorMsg_arg)
completeErrorMsg(const string &ErrorMsg_arg)
{
ErrorMsg += ErrorMsg_arg;
}
@ -69,9 +67,9 @@ public:
class FloatingPointExceptionHandling : public GeneralExceptionHandling
{
public:
FloatingPointExceptionHandling(string value) : GeneralExceptionHandling(string("Floating point error in bytecode: " + value))
FloatingPointExceptionHandling(const string &value) : GeneralExceptionHandling("Floating point error in bytecode: " + value)
{
};
}
};
class LogExceptionHandling : public FloatingPointExceptionHandling
@ -81,10 +79,8 @@ public:
LogExceptionHandling(double value_arg) : FloatingPointExceptionHandling("log(X)"),
value(value_arg)
{
ostringstream tmp;
tmp << " with X=" << value << "\n";
completeErrorMsg(tmp.str());
};
completeErrorMsg(" with X=" + to_string(value) + "\n");
}
};
class Log10ExceptionHandling : public FloatingPointExceptionHandling
@ -94,10 +90,8 @@ public:
Log10ExceptionHandling(double value_arg) : FloatingPointExceptionHandling("log10(X)"),
value(value_arg)
{
ostringstream tmp;
tmp << " with X=" << value << "\n";
completeErrorMsg(tmp.str());
};
completeErrorMsg(" with X=" + to_string(value) + "\n");
}
};
class DivideExceptionHandling : public FloatingPointExceptionHandling
@ -108,10 +102,8 @@ public:
value1(value1_arg),
value2(value2_arg)
{
ostringstream tmp;
tmp << " with X=" << value2 << "\n";
completeErrorMsg(tmp.str());
};
completeErrorMsg(" with X=" + to_string(value2) + "\n");
}
};
class PowExceptionHandling : public FloatingPointExceptionHandling
@ -122,12 +114,10 @@ public:
value1(value1_arg),
value2(value2_arg)
{
ostringstream tmp;
if (fabs(value1) > 1e-10)
tmp << " with X=" << value1 << "\n";
completeErrorMsg(" with X=" + to_string(value1) + "\n");
else
tmp << " with X=" << value1 << " and a=" << value2 << "\n";
completeErrorMsg(tmp.str());
completeErrorMsg(" with X=" + to_string(value1) + " and a=" + to_string(value2) + "\n");
};
};
@ -144,7 +134,8 @@ public:
class FatalExceptionHandling : public GeneralExceptionHandling
{
public:
FatalExceptionHandling(string ErrorMsg_arg) : GeneralExceptionHandling("Fatal error in bytecode:")
FatalExceptionHandling(const string &ErrorMsg_arg)
: GeneralExceptionHandling("Fatal error in bytecode:")
{
completeErrorMsg(ErrorMsg_arg);
};
@ -191,8 +182,8 @@ public:
vector<s_plan> splan, spfplan;
vector<mxArray *> jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block;
map<unsigned int, double> TEF;
map<pair<unsigned int, unsigned int>, double > TEFD;
map<pair<unsigned int, pair<unsigned int, unsigned int>>, double > TEFDD;
map<pair<unsigned int, unsigned int>, double> TEFD;
map<tuple<unsigned int, unsigned int, unsigned int>, double> TEFDD;
ExpressionType EQN_type;
it_code_type it_code_expr;
@ -201,7 +192,7 @@ public:
size_t /*unsigned int*/ endo_name_length, exo_name_length, param_name_length;
unsigned int EQN_equation, EQN_block, EQN_block_number;
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
vector<pair<string, pair<SymbolType, unsigned int>>> Variable_list;
vector<tuple<string, SymbolType, unsigned int>> Variable_list;
inline
ErrorMsg()
@ -263,9 +254,9 @@ public:
else
{
if (dollar.compare(&i) == 0)
pos1 = int (temp.length());
pos1 = static_cast<int>(temp.length());
else
pos2 = int (temp.length());
pos2 = static_cast<int>(temp.length());
if (pos1 >= 0 && pos2 >= 0)
{
tmp_n.erase(pos1, pos2-pos1+1);
@ -287,19 +278,19 @@ public:
for (unsigned int i = 0; i < endo_name_length; i++)
if (P_endo_names[CHAR_LENGTH*(variable_num+i*nb_endo)] != ' ')
res << P_endo_names[CHAR_LENGTH*(variable_num+i*nb_endo)];
Variable_list.emplace_back(res.str(), make_pair(SymbolType::endogenous, variable_num));
Variable_list.emplace_back(res.str(), SymbolType::endogenous, variable_num);
}
for (unsigned int variable_num = 0; variable_num < static_cast<unsigned int>(nb_exo); variable_num++)
{
for (unsigned int i = 0; i < exo_name_length; i++)
if (P_exo_names[CHAR_LENGTH*(variable_num+i*nb_exo)] != ' ')
res << P_exo_names[CHAR_LENGTH*(variable_num+i*nb_exo)];
Variable_list.emplace_back(res.str(), make_pair(SymbolType::exogenous, variable_num));
Variable_list.emplace_back(res.str(), SymbolType::exogenous, variable_num);
}
}
inline int
get_ID(const string variable_name, SymbolType *variable_type)
get_ID(const string &variable_name, SymbolType *variable_type)
{
if (!is_load_variable_list)
{
@ -311,19 +302,19 @@ public:
bool notfound = true;
while (notfound && i < static_cast<int>(n))
{
if (variable_name == Variable_list[i].first)
if (variable_name == get<0>(Variable_list[i]))
{
notfound = false;
*variable_type = Variable_list[i].second.first;
return Variable_list[i].second.second;
*variable_type = get<1>(Variable_list[i]);
return get<2>(Variable_list[i]);
}
i++;
}
return (-1);
return -1;
}
inline string
get_variable(const SymbolType variable_type, const unsigned int variable_num) const
get_variable(SymbolType variable_type, unsigned int variable_num) const
{
ostringstream res;
switch (variable_type)
@ -362,13 +353,13 @@ public:
default:
break;
}
return (res.str());
return res.str();
}
inline string
error_location(bool evaluate, bool steady_state, int size, int block_num, int it_, int Per_u_)
{
stringstream Error_loc("");
stringstream Error_loc;
if (!steady_state)
switch (EQN_type)
{
@ -467,7 +458,7 @@ public:
}
it_code_type it_code_ret;
Error_loc << endl << add_underscore_to_fpe(" " + print_expression(it_code_expr, evaluate, size, block_num, steady_state, Per_u_, it_, it_code_ret, true));
return (Error_loc.str());
return Error_loc.str();
}
inline string
@ -1214,7 +1205,7 @@ public:
mexPrintf("<");
#endif
if (compute)
Stackf.push(double (v1f < v2f));
Stackf.push(static_cast<double>(v1f < v2f));
tmp_out.str("");
found = v1.find(" ");
if (found != string::npos)
@ -1236,7 +1227,7 @@ public:
mexPrintf(">");
#endif
if (compute)
Stackf.push(double (v1f > v2f));
Stackf.push(static_cast<double>(v1f > v2f));
tmp_out.str("");
found = v1.find(" ");
if (found != string::npos)
@ -1258,7 +1249,7 @@ public:
mexPrintf("<=");
#endif
if (compute)
Stackf.push(double (v1f <= v2f));
Stackf.push(static_cast<double>(v1f <= v2f));
tmp_out.str("");
found = v1.find(" ");
if (found != string::npos)
@ -1280,7 +1271,7 @@ public:
mexPrintf(">=");
#endif
if (compute)
Stackf.push(double (v1f >= v2f));
Stackf.push(static_cast<double>(v1f >= v2f));
tmp_out.str("");
found = v1.find(" ");
if (found != string::npos)
@ -1302,7 +1293,7 @@ public:
mexPrintf("==");
#endif
if (compute)
Stackf.push(double (v1f == v2f));
Stackf.push(static_cast<double>(v1f == v2f));
tmp_out.str("");
found = v1.find(" ");
if (found != string::npos)
@ -1324,7 +1315,7 @@ public:
mexPrintf("!=");
#endif
if (compute)
Stackf.push(double (v1f != v2f));
Stackf.push(static_cast<double>(v1f != v2f));
tmp_out.str("");
found = v1.find(" ");
if (found != string::npos)
@ -1380,7 +1371,7 @@ public:
Stack.pop();
if (compute)
{
int derivOrder = int (nearbyint(Stackf.top()));
int derivOrder = static_cast<int>(nearbyint(Stackf.top()));
Stackf.pop();
if (fabs(v1f) < near_zero && v2f > 0
&& derivOrder > v2f
@ -1967,7 +1958,7 @@ public:
{
#ifdef DEBUG
double rr = Stackf.top();
mexPrintf("FSTP TEFD[make_pair(indx, row)]=%f done\n", rr);
mexPrintf("FSTP TEFD[{ indx, row }]=%f done\n", rr);
mexEvalString("drawnow;");
#endif
Stackf.pop();
@ -1987,13 +1978,13 @@ public:
#ifdef DEBUG
mexPrintf("FLDTEFD\n");
mexPrintf("indx=%d row=%d Stack.size()=%d\n", indx, row, Stack.size());
auto it = TEFD.find(make_pair(indx, row-1));
mexPrintf("FLD TEFD[make_pair(indx, row)]=%f done\n", it->second);
auto it = TEFD.find({ indx, row-1 });
mexPrintf("FLD TEFD[{ indx, row }]=%f done\n", it->second);
mexEvalString("drawnow;");
#endif
if (compute)
{
auto it = TEFD.find(make_pair(indx, row-1));
auto it = TEFD.find({ indx, row-1 });
Stackf.push(it->second);
}
tmp_out.str("");
@ -2016,8 +2007,8 @@ public:
#ifdef DEBUG
double rr = Stackf.top();
mexPrintf("rr=%f\n", rr);
auto it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1)));
mexPrintf("FSTP TEFDD[make_pair(indx, make_pair(row, col))]=%f done\n", it->second);
auto it = TEFDD.find({ indx, row-1, col-1 });
mexPrintf("FSTP TEFDD[{ indx, row, col }]=%f done\n", it->second);
mexEvalString("drawnow;");
#endif
Stackf.pop();
@ -2039,13 +2030,13 @@ public:
#ifdef DEBUG
mexPrintf("FLDTEFD\n");
mexPrintf("indx=%d Stack.size()=%d\n", indx, Stack.size());
auto it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1)));
mexPrintf("FLD TEFD[make_pair(indx, make_pair(row, col))]=%f done\n", it->second);
auto it = TEFDD.find({ indx, row-1, col-1 });
mexPrintf("FLD TEFD[{ indx, row, col }]=%f done\n", it->second);
mexEvalString("drawnow;");
#endif
if (compute)
{
auto it = TEFDD.find(make_pair(indx, make_pair(row-1, col-1)));
auto it = TEFDD.find({ indx, row-1, col-1 });
Stackf.push(it->second);
}
tmp_out.str("");
@ -2087,10 +2078,11 @@ public:
case Tags::FOK:
break;
default:
ostringstream tmp;
mexPrintf("Error it_code->first=%d unknown\n", it_code->first); mexEvalString("drawnow;");
tmp << " in print_expression, unknown opcode " << static_cast<int>(it_code->first) << "!! FENDEQU=" << static_cast<int>(Tags::FENDEQU) << "\n";
throw FatalExceptionHandling(tmp.str());
throw FatalExceptionHandling(" in print_expression, unknown opcode "
+ to_string(static_cast<int>(it_code->first))
+ "!! FENDEQU="
+ to_string(static_cast<int>(Tags::FENDEQU)) + "\n");
}
it_code++;
}
@ -2098,21 +2090,15 @@ public:
mexPrintf("print_expression end tmp_out.str().c_str()=%s\n", tmp_out.str().c_str()); mexEvalString("drawnow;");
#endif
it_code_ret = it_code;
return (tmp_out.str());
return tmp_out.str();
}
void
inline
test_mxMalloc(void *z, int line, string file, string func, int amount)
static inline void
test_mxMalloc(void *z, int line, const string &file, const string &func, int amount)
{
if (z == nullptr && (amount > 0))
{
ostringstream tmp;
tmp << " mxMalloc: out of memory " << amount << " bytes required at line " << line << " in function " << func << " (file " << file;
throw FatalExceptionHandling(tmp.str());
}
if (!z && amount > 0)
throw FatalExceptionHandling(" mxMalloc: out of memory " + to_string(amount) + " bytes required at line " + to_string(line) + " in function " + func + " (file " + file);
}
};
#endif

228
mex/sources/bytecode/Evaluate.cc
View File

@ -20,6 +20,8 @@
#include <cstring>
#include <sstream>
#include <cmath>
#include <limits>
#include "Evaluate.hh"
#ifdef MATLAB_MEX_FILE
@ -35,7 +37,7 @@ Evaluate::Evaluate()
block = -1;
}
Evaluate::Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg, const double slowc_arg) :
Evaluate::Evaluate(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, int periods_arg, int minimal_solving_periods_arg, double slowc_arg) :
print_it(print_it_arg), minimal_solving_periods(minimal_solving_periods_arg)
{
symbol_table_endo_nbr = 0;
@ -54,10 +56,10 @@ Evaluate::Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_
double
Evaluate::pow1(double a, double b)
{
double r = pow_(a, b);
double r = pow(a, b);
if (isnan(r) || isinf(r))
{
res1 = NAN;
res1 = std::numeric_limits<double>::quiet_NaN();
r = 0.0000000000000000000000001;
if (print_error)
throw PowExceptionHandling(a, b);
@ -71,7 +73,7 @@ Evaluate::divide(double a, double b)
double r = a / b;
if (isnan(r) || isinf(r))
{
res1 = NAN;
res1 = std::numeric_limits<double>::quiet_NaN();
r = 1e70;
if (print_error)
throw DivideExceptionHandling(a, b);
@ -85,7 +87,7 @@ Evaluate::log1(double a)
double r = log(a);
if (isnan(r) || isinf(r))
{
res1 = NAN;
res1 = std::numeric_limits<double>::quiet_NaN();
r = -1e70;
if (print_error)
throw LogExceptionHandling(a);
@ -99,7 +101,7 @@ Evaluate::log10_1(double a)
double r = log(a);
if (isnan(r) || isinf(r))
{
res1 = NAN;
res1 = std::numeric_limits<double>::quiet_NaN();
r = -1e70;
if (print_error)
throw Log10ExceptionHandling(a);
@ -108,7 +110,7 @@ Evaluate::log10_1(double a)
}
void
Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int block_num, const int size, const bool steady_state,*/ const bool no_derivative)
Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
{
int var = 0, lag = 0, op;
unsigned int eq, pos_col;
@ -702,11 +704,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
//store in the jacobian matrix
rr = Stack.top();
if (EQN_type != ExpressionType::FirstEndoDerivative)
{
ostringstream tmp;
tmp << " in compute_block_time, impossible case " << static_cast<int>(EQN_type) << " not implement in static jacobian\n";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" in compute_block_time, impossible case " + to_string(static_cast<int>(EQN_type)) + " not implement in static jacobian\n");
eq = static_cast<FSTPG2_ *>(it_code->second)->get_row();
var = static_cast<FSTPG2_ *>(it_code->second)->get_col();
#ifdef DEBUG
@ -779,9 +777,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
jacob_exo_det[eq + size*pos_col] = rr;
break;
default:
ostringstream tmp;
tmp << " in compute_block_time, variable " << static_cast<int>(EQN_type) << " not used yet\n";
throw FatalExceptionHandling(tmp.str());
throw FatalExceptionHandling(" in compute_block_time, variable " + to_string(static_cast<int>(EQN_type)) + " not used yet\n");
}
// #ifdef DEBUG
// tmp_out << "=>";
@ -840,37 +836,37 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
#endif
break;
case BinaryOpcode::less:
Stack.push(double (v1 < v2));
Stack.push(static_cast<double>(v1 < v2));
#ifdef DEBUG
mexPrintf("v1=%f v2=%f v1 < v2 = %f\n", v1, v2, double (v1 < v2));
mexPrintf("v1=%f v2=%f v1 < v2 = %f\n", v1, v2, static_cast<double>(v1 < v2));
#endif
break;
case BinaryOpcode::greater:
Stack.push(double (v1 > v2));
Stack.push(static_cast<double>(v1 > v2));
#ifdef DEBUG
tmp_out << " |" << v1 << ">" << v2 << "|";
#endif
break;
case BinaryOpcode::lessEqual:
Stack.push(double (v1 <= v2));
Stack.push(static_cast<double>(v1 <= v2));
#ifdef DEBUG
tmp_out << " |" << v1 << "<=" << v2 << "|";
#endif
break;
case BinaryOpcode::greaterEqual:
Stack.push(double (v1 >= v2));
Stack.push(static_cast<double>(v1 >= v2));
#ifdef DEBUG
tmp_out << " |" << v1 << ">=" << v2 << "|";
#endif
break;
case BinaryOpcode::equalEqual:
Stack.push(double (v1 == v2));
Stack.push(static_cast<double>(v1 == v2));
#ifdef DEBUG
tmp_out << " |" << v1 << "==" << v2 << "|";
#endif
break;
case BinaryOpcode::different:
Stack.push(double (v1 != v2));
Stack.push(static_cast<double>(v1 != v2));
#ifdef DEBUG
tmp_out << " |" << v1 << "!=" << v2 << "|";
#endif
@ -896,7 +892,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
break;
case BinaryOpcode::powerDeriv:
{
int derivOrder = int (nearbyint(Stack.top()));
int derivOrder = static_cast<int>(nearbyint(Stack.top()));
Stack.pop();
try
{
@ -941,9 +937,8 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
default:
{
mexPrintf("Error\n");
ostringstream tmp;
tmp << " in compute_block_time, unknown binary operator " << op << "\n";
throw FatalExceptionHandling(tmp.str());
throw FatalExceptionHandling(" in compute_block_time, unknown binary operator "
+ to_string(op) + "\n");
}
}
break;
@ -981,7 +976,6 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
go_on = false;
}
Stack.push(tmp);
//if (isnan(res1))
#ifdef DEBUG
tmp_out << " |log(" << v1 << ")|";
@ -1090,9 +1084,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
default:
{
mexPrintf("Error\n");
ostringstream tmp;
tmp << " in compute_block_time, unknown unary operator " << op << "\n";
throw FatalExceptionHandling(tmp.str());
throw FatalExceptionHandling(" in compute_block_time, unknown unary operator " + to_string(op) + "\n");
}
}
break;
@ -1121,9 +1113,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
default:
{
mexPrintf("Error\n");
ostringstream tmp;
tmp << " in compute_block_time, unknown trinary operator " << op << "\n";
throw FatalExceptionHandling(tmp.str());
throw FatalExceptionHandling(" in compute_block_time, unknown trinary operator " + to_string(op) + "\n");
}
}
break;
@ -1183,11 +1173,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
Stack.pop();
}
if (mexCallMATLAB(nb_output_arguments, output_arguments, nb_input_arguments, input_arguments, function_name.c_str()))
{
ostringstream tmp;
tmp << " external function: " << function_name << " not found";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" external function: " + function_name + " not found");
double *rr = mxGetPr(output_arguments[0]);
Stack.push(*rr);
@ -1197,7 +1183,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
double *FD1 = mxGetPr(output_arguments[1]);
size_t rows = mxGetN(output_arguments[1]);
for (unsigned int i = 0; i < rows; i++)
TEFD[make_pair(indx, i)] = FD1[i];
TEFD[{ indx, i }] = FD1[i];
}
if (function_type == ExternalFunctionType::withFirstAndSecondDerivative)
{
@ -1208,7 +1194,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
unsigned int k = 0;
for (unsigned int j = 0; j < cols; j++)
for (unsigned int i = 0; i < rows; i++)
TEFDD[make_pair(indx, make_pair(i, j))] = FD2[k++];
TEFDD[{ indx, i, j }] = FD2[k++];
}
}
break;
@ -1240,11 +1226,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
#endif
nb_input_arguments = 3;
if (mexCallMATLAB(nb_output_arguments, output_arguments, nb_input_arguments, input_arguments, function_name.c_str()))
{
ostringstream tmp;
tmp << " external function: " << function_name << " not found";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" external function: " + function_name + " not found");
double *rr = mxGetPr(output_arguments[0]);
#ifdef DEBUG
mexPrintf("*rr=%f\n", *rr);
@ -1263,17 +1245,12 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
Stack.pop();
}
if (mexCallMATLAB(nb_output_arguments, output_arguments, nb_input_arguments, input_arguments, function_name.c_str()))
{
ostringstream tmp;
tmp << " external function: " << function_name << " not found";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" external function: " + function_name + " not found");
unsigned int indx = fc->get_indx();
double *FD1 = mxGetPr(output_arguments[0]);
//mexPrint
size_t rows = mxGetN(output_arguments[0]);
for (unsigned int i = 0; i < rows; i++)
TEFD[make_pair(indx, i)] = FD1[i];
TEFD[{ indx, i }] = FD1[i];
}
break;
case ExternalFunctionType::numericalSecondDerivative:
@ -1306,11 +1283,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
#endif
nb_input_arguments = 3;
if (mexCallMATLAB(nb_output_arguments, output_arguments, nb_input_arguments, input_arguments, function_name.c_str()))
{
ostringstream tmp;
tmp << " external function: " << function_name << " not found";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" external function: " + function_name + " not found");
double *rr = mxGetPr(output_arguments[0]);
Stack.push(*rr);
}
@ -1326,11 +1299,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
Stack.pop();
}
if (mexCallMATLAB(nb_output_arguments, output_arguments, nb_input_arguments, input_arguments, function_name.c_str()))
{
ostringstream tmp;
tmp << " external function: " << function_name << " not found";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" external function: " + function_name + " not found");
unsigned int indx = fc->get_indx();
double *FD2 = mxGetPr(output_arguments[2]);
size_t rows = mxGetM(output_arguments[0]);
@ -1338,7 +1307,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
unsigned int k = 0;
for (unsigned int j = 0; j < cols; j++)
for (unsigned int i = 0; i < rows; i++)
TEFDD[make_pair(indx, make_pair(i, j))] = FD2[k++];
TEFDD[{ indx, i, j }] = FD2[k++];
}
break;
}
@ -1377,9 +1346,9 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
#endif
if (function_type == ExternalFunctionType::numericalFirstDerivative)
{
TEFD[make_pair(indx, row-1)] = Stack.top();
TEFD[{ indx, row-1 }] = Stack.top();
#ifdef DEBUG
mexPrintf("FSTP TEFD[make_pair(indx, row)]=%f done\n", TEFD[make_pair(indx, row-1)]);
mexPrintf("FSTP TEFD[{ indx, row }]=%f done\n", TEFD[{ indx, row-1 }]);
mexEvalString("drawnow;");
#endif
Stack.pop();
@ -1394,10 +1363,10 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
#ifdef DEBUG
mexPrintf("FLDTEFD\n");
mexPrintf("indx=%d row=%d Stack.size()=%d\n", indx, row, Stack.size());
mexPrintf("FLD TEFD[make_pair(indx, row)]=%f done\n", TEFD[make_pair(indx, row-1)]);
mexPrintf("FLD TEFD[{ indx, row }]=%f done\n", TEFD[{ indx, row-1 }]);
mexEvalString("drawnow;");
#endif
Stack.push(TEFD[make_pair(indx, row-1)]);
Stack.push(TEFD[{ indx, row-1 }]);
}
break;
case Tags::FSTPTEFDD:
@ -1411,9 +1380,9 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
#endif
if (function_type == ExternalFunctionType::numericalSecondDerivative)
{
TEFDD[make_pair(indx, make_pair(row-1, col-1))] = Stack.top();
TEFDD[{ indx, row-1, col-1 }] = Stack.top();
#ifdef DEBUG
mexPrintf("FSTP TEFDD[make_pair(indx, make_pair(row, col))]=%f done\n", TEFDD[make_pair(indx, make_pair(row, col))]);
mexPrintf("FSTP TEFDD[{ indx, row, col }]=%f done\n", TEFDD[{ indx, row, col }]);
mexEvalString("drawnow;");
#endif
Stack.pop();
@ -1429,10 +1398,10 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
#ifdef DEBUG
mexPrintf("FLDTEFD\n");
mexPrintf("indx=%d Stack.size()=%d\n", indx, Stack.size());
mexPrintf("FLD TEFD[make_pair(indx, make_pair(row, col))]=%f done\n", TEFDD[make_pair(indx, make_pair(row, col))]);
mexPrintf("FLD TEFD[{ indx, row, col }]=%f done\n", TEFDD[{ indx, row, col }]);
mexEvalString("drawnow;");
#endif
Stack.push(TEFDD[make_pair(indx, make_pair(row-1, col-1))]);
Stack.push(TEFDD[{ indx, row-1, col-1 }]);
}
break;
case Tags::FCUML:
@ -1476,16 +1445,10 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
case Tags::FOK:
op = static_cast<FOK_ *>(it_code->second)->get_arg();
if (Stack.size() > 0)
{
ostringstream tmp;
tmp << " in compute_block_time, stack not empty\n";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" in compute_block_time, stack not empty\n");
break;
default:
ostringstream tmp;
tmp << " in compute_block_time, unknown opcode " << static_cast<int>(it_code->first) << "\n";
throw FatalExceptionHandling(tmp.str());
throw FatalExceptionHandling(" in compute_block_time, unknown opcode " + to_string(static_cast<int>(it_code->first)) + "\n");
}
#ifdef DEBUG
mexPrintf("it_code++=%d\n", it_code);
@ -1499,7 +1462,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
}
void
Evaluate::evaluate_over_periods(const bool forward)
Evaluate::evaluate_over_periods(bool forward)
{
if (steady_state)
compute_block_time(0, false, false);
@ -1535,7 +1498,7 @@ Evaluate::solve_simple_one_periods()
double ya;
double slowc_save = slowc;
res1 = 0;
while (!(cvg || (iter > maxit_)))
while (!(cvg || iter > maxit_))
{
it_code = start_code;
Per_y_ = it_*y_size;
@ -1543,7 +1506,7 @@ Evaluate::solve_simple_one_periods()
compute_block_time(0, false, false);
if (!isfinite(res1))
{
res1 = NAN;
res1 = std::numeric_limits<double>::quiet_NaN();
while ((isinf(res1) || isnan(res1)) && (slowc > 1e-9))
{
it_code = start_code;
@ -1559,7 +1522,6 @@ Evaluate::solve_simple_one_periods()
double rr;
rr = r[0];
cvg = (fabs(rr) < solve_tolf);
//mexPrintf("g1=%x, g1[0]=%f, type=%d, block_num=%d, it_=%d y=%x\n", g1, g1[0], type, block_num, it_, y);
if (cvg)
continue;
try
@ -1575,15 +1537,12 @@ Evaluate::solve_simple_one_periods()
}
slowc = slowc_save;
if (!cvg)
{
ostringstream tmp;
tmp << " in Solve Forward simple, convergence not achieved in block " << block_num+1 << ", after " << iter << " iterations\n";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" in Solve Forward simple, convergence not achieved in block "
+ to_string(block_num+1) + ", after " + to_string(iter) + " iterations\n");
}
void
Evaluate::solve_simple_over_periods(const bool forward)
Evaluate::solve_simple_over_periods(bool forward)
{
g1 = static_cast<double *>(mxMalloc(sizeof(double)));
test_mxMalloc(g1, __LINE__, __FILE__, __func__, sizeof(double));
@ -1615,13 +1574,13 @@ Evaluate::solve_simple_over_periods(const bool forward)
}
void
Evaluate::set_block(const int size_arg, const int type_arg, string file_name_arg, string bin_base_name_arg, const int block_num_arg,
const bool is_linear_arg, const int symbol_table_endo_nbr_arg, const int Block_List_Max_Lag_arg, const int Block_List_Max_Lead_arg, const int u_count_int_arg, const int block_arg)
Evaluate::set_block(int size_arg, int type_arg, string file_name_arg, string bin_base_name_arg, int block_num_arg,
bool is_linear_arg, int symbol_table_endo_nbr_arg, int Block_List_Max_Lag_arg, int Block_List_Max_Lead_arg, int u_count_int_arg, int block_arg)
{
size = size_arg;
type = type_arg;
file_name = file_name_arg;
bin_base_name = bin_base_name_arg;
file_name = move(file_name_arg);
bin_base_name = move(bin_base_name_arg);
block_num = block_num_arg;
is_linear = is_linear_arg;
symbol_table_endo_nbr = symbol_table_endo_nbr_arg;
@ -1632,14 +1591,14 @@ Evaluate::set_block(const int size_arg, const int type_arg, string file_name_arg
}
void
Evaluate::evaluate_complete(const bool no_derivatives)
Evaluate::evaluate_complete(bool no_derivatives)
{
it_code = start_code;
compute_block_time(0, false, no_derivatives);
}
void
Evaluate::compute_complete_2b(const bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx)
Evaluate::compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx)
{
res1 = 0;
*_res1 = 0;
@ -1653,23 +1612,19 @@ Evaluate::compute_complete_2b(const bool no_derivatives, double *_res1, double *
int shift = (it_-y_kmin) * size;
compute_block_time(Per_u_, false, no_derivatives);
if (!(isnan(res1) || isinf(res1)))
{
for (int i = 0; i < size; i++)
{
for (int i = 0; i < size; i++)
double rr;
rr = r[i];
res[i+shift] = rr;
if (max_res < fabs(rr))
{
double rr;
rr = r[i];
res[i+shift] = rr;
if (max_res < fabs(rr))
{
*_max_res = fabs(rr);
*_max_res_idx = i;
}
*_res2 += rr*rr;
*_res1 += fabs(rr);
*_max_res = fabs(rr);
*_max_res_idx = i;
}
*_res2 += rr*rr;
*_res1 += fabs(rr);
}
}
else
return;
}
@ -1678,7 +1633,7 @@ Evaluate::compute_complete_2b(const bool no_derivatives, double *_res1, double *
}
bool
Evaluate::compute_complete(const bool no_derivatives, double &_res1, double &_res2, double &_max_res, int &_max_res_idx)
Evaluate::compute_complete(bool no_derivatives, double &_res1, double &_res2, double &_max_res, int &_max_res_idx)
{
bool result;
res1 = 0;
@ -1686,23 +1641,21 @@ Evaluate::compute_complete(const bool no_derivatives, double &_res1, double &_re
compute_block_time(0, false, no_derivatives);
if (!(isnan(res1) || isinf(res1)))
{
{
_res1 = 0;
_res2 = 0;
_max_res = 0;
for (int i = 0; i < size; i++)
{
double rr;
rr = r[i];
if (max_res < fabs(rr))
{
_max_res = fabs(rr);
_max_res_idx = i;
}
_res2 += rr*rr;
_res1 += fabs(rr);
}
}
_res1 = 0;
_res2 = 0;
_max_res = 0;
for (int i = 0; i < size; i++)
{
double rr;
rr = r[i];
if (max_res < fabs(rr))
{
_max_res = fabs(rr);
_max_res_idx = i;
}
_res2 += rr*rr;
_res1 += fabs(rr);
}
result = true;
}
else
@ -1714,7 +1667,6 @@ bool
Evaluate::compute_complete(double lambda, double *crit)
{
double res1_ = 0, res2_ = 0, max_res_ = 0;
//double res1 = 0, res2, max_res;
int max_res_idx_ = 0;
if (steady_state)
{
@ -1727,28 +1679,18 @@ Evaluate::compute_complete(double lambda, double *crit)
Per_u_ = 0;
Per_y_ = 0;
if (compute_complete(true, res1, res2, max_res, max_res_idx))
{
res2_ = res2;
/*res1_ = res1;
if (max_res > max_res_)
{
max_res = max_res_;
max_res_idx = max_res_idx_;
}*/
}
res2_ = res2;
else
return false;
}
else
{
for (int it = y_kmin; it < periods+y_kmin; it++)
{
for (int i = 0; i < size; i++)
{
int eq = index_vara[i];
y[eq+it*y_size] = ya[eq+it*y_size] + lambda * direction[eq+it*y_size];
}
}
for (int i = 0; i < size; i++)
{
int eq = index_vara[i];
y[eq+it*y_size] = ya[eq+it*y_size] + lambda * direction[eq+it*y_size];
}
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
{
Per_u_ = (it_-y_kmin)*u_count_int;

24
mex/sources/bytecode/Evaluate.hh
View File

@ -1,5 +1,5 @@
/*
* Copyright © 2007-2017 Dynare Team
* Copyright © 2007-2021 Dynare Team
*
* This file is part of Dynare.
*
@ -29,8 +29,6 @@
#include <dynmex.h>
#include "ErrorHandling.hh"
#define pow_ pow
class Evaluate : public ErrorMsg
{
private:
@ -43,10 +41,10 @@ protected:
double divide(double a, double b);
double log1(double a);
double log10_1(double a);
void evaluate_over_periods(const bool forward);
void evaluate_over_periods(bool forward);
void solve_simple_one_periods();
void solve_simple_over_periods(const bool forward);
void compute_block_time(const int Per_u_, const bool evaluate, const bool no_derivatives);
void solve_simple_over_periods(bool forward);
void compute_block_time(int Per_u_, bool evaluate, bool no_derivatives);
code_liste_type code_liste;
it_code_type it_code;
int Block_Count, Per_u_, Per_y_;
@ -55,7 +53,6 @@ protected:
double *direction;
double solve_tolf;
bool GaussSeidel;
map<pair<pair<int, int>, int>, int> IM_i;
int equation, derivative_equation, derivative_variable;
string filename;
int stack_solve_algo, solve_algo;
@ -77,13 +74,12 @@ public:
bool steady_state;
double slowc;
Evaluate();
Evaluate(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg, const double slowc);
//typedef void (Interpreter::*InterfpreterMemFn)(const int block_num, const int size, const bool steady_state, int it);
void set_block(const int size_arg, const int type_arg, string file_name_arg, string bin_base_name_arg, const int block_num_arg,
const bool is_linear_arg, const int symbol_table_endo_nbr_arg, const int Block_List_Max_Lag_arg, const int Block_List_Max_Lead_arg, const int u_count_int_arg, const int block_arg);
void evaluate_complete(const bool no_derivatives);
bool compute_complete(const bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
void compute_complete_2b(const bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx);
Evaluate(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, int periods_arg, int minimal_solving_periods_arg, double slowc);
void set_block(int size_arg, int type_arg, string file_name_arg, string bin_base_name_arg, int block_num_arg,
bool is_linear_arg, int symbol_table_endo_nbr_arg, int Block_List_Max_Lag_arg, int Block_List_Max_Lead_arg, int u_count_int_arg, int block_arg);
void evaluate_complete(bool no_derivatives);
bool compute_complete(bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
void compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx);
bool compute_complete(double lambda, double *crit);
};

156
mex/sources/bytecode/Interpreter.cc
View File

@ -21,11 +21,8 @@
#include <sstream>
#include <algorithm>
#include "Interpreter.hh"
#define BIG 1.0e+8;
#define SMALL 1.0e-5;
///#define DEBUG
Interpreter::~Interpreter() = default;
constexpr double BIG = 1.0e+8, SMALL = 1.0e-5;
Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
double *direction_arg, size_t y_size_arg,
@ -43,12 +40,9 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
steady_y = steady_y_arg;
steady_x = steady_x_arg;
direction = direction_arg;
//y_size = y_size_arg;
nb_row_x = nb_row_x_arg;
nb_row_xd = nb_row_xd_arg;
periods = periods_arg;
//y_kmax = y_kmax_arg;
//y_kmin = y_kmin_arg;
maxit_ = maxit_arg_;
solve_tolf = solve_tolf_arg;
size_of_direction = size_of_direction_arg;
@ -67,9 +61,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
col_y = col_y_arg;
GlobalTemporaryTerms = GlobalTemporaryTerms_arg;
print_error = print_error_arg;
//steady_state = steady_state_arg;
print_it = print_it_arg;
}
void
@ -82,7 +74,7 @@ Interpreter::evaluate_a_block(bool initialization)
case BlockSimulationType::evaluateForward:
if (steady_state)
{
compute_block_time(0, true, /*block_num, size, steady_state, */ false);
compute_block_time(0, true, false);
if (block >= 0)
for (int j = 0; j < size; j++)
residual[j] = y[Block_Contain[j].Variable] - ya[Block_Contain[j].Variable];
@ -97,7 +89,7 @@ Interpreter::evaluate_a_block(bool initialization)
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, /*block_num, size, steady_state, */ false);
compute_block_time(0, true, false);
if (block >= 0)
for (int j = 0; j < size; j++)
residual[it_*size+j] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
@ -114,7 +106,7 @@ Interpreter::evaluate_a_block(bool initialization)
test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double));
if (steady_state)
{
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[Block_Contain[j].Equation] = r[j];
@ -129,7 +121,7 @@ Interpreter::evaluate_a_block(bool initialization)
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[Per_y_+Block_Contain[j].Equation] = r[j];
@ -154,7 +146,7 @@ Interpreter::evaluate_a_block(bool initialization)
test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double));
if (steady_state)
{
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[Block_Contain[j].Equation] = r[j];
@ -169,7 +161,7 @@ Interpreter::evaluate_a_block(bool initialization)
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[it_*y_size+Block_Contain[j].Equation] = r[j];
@ -183,7 +175,7 @@ Interpreter::evaluate_a_block(bool initialization)
case BlockSimulationType::evaluateBackward:
if (steady_state)
{
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block >= 0)
for (int j = 0; j < size; j++)
residual[j] = y[Block_Contain[j].Variable] - ya[Block_Contain[j].Variable];
@ -198,7 +190,7 @@ Interpreter::evaluate_a_block(bool initialization)
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block >= 0)
for (int j = 0; j < size; j++)
residual[it_*size+j] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
@ -215,7 +207,7 @@ Interpreter::evaluate_a_block(bool initialization)
test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double));
if (steady_state)
{
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[Block_Contain[j].Equation] = r[j];
@ -230,7 +222,7 @@ Interpreter::evaluate_a_block(bool initialization)
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[Per_y_+Block_Contain[j].Equation] = r[j];
@ -252,7 +244,7 @@ Interpreter::evaluate_a_block(bool initialization)
test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double));
if (steady_state)
{
compute_block_time(0, true, /*block_num, size, steady_state,*/ false);
compute_block_time(0, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[Block_Contain[j].Equation] = r[j];
@ -267,7 +259,7 @@ Interpreter::evaluate_a_block(bool initialization)
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, /*block_num, size, steady_state, */ false);
compute_block_time(0, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[Per_y_+Block_Contain[j].Equation] = r[j];
@ -294,7 +286,7 @@ Interpreter::evaluate_a_block(bool initialization)
Per_u_ = (it_-y_kmin)*u_count_int;
Per_y_ = it_*y_size;
it_code = begining;
compute_block_time(Per_u_, true, /*block_num, size, steady_state,*/ false);
compute_block_time(Per_u_, true, false);
if (block < 0)
for (int j = 0; j < size; j++)
residual[it_*y_size+Block_Contain[j].Equation] = r[j];
@ -310,7 +302,7 @@ Interpreter::evaluate_a_block(bool initialization)
}
int
Interpreter::simulate_a_block(vector_table_conditional_local_type vector_table_conditional_local)
Interpreter::simulate_a_block(const vector_table_conditional_local_type &vector_table_conditional_local)
{
it_code_type begining;
max_res = 0;
@ -450,17 +442,9 @@ Interpreter::simulate_a_block(vector_table_conditional_local_type vector_table_c
max_res_idx = 0;
memcpy(y_save, y, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
if (vector_table_conditional_local.size())
{
for (auto & it1 : vector_table_conditional_local)
{
if (it1.is_cond)
{
//mexPrintf("y[%d] = %f\n", it1->var_endo + y_kmin * size, y[it1->var_endo + y_kmin * size]);
y[it1.var_endo + y_kmin * size] = it1.constrained_value;
}
}
}
for (auto & it1 : vector_table_conditional_local)
if (it1.is_cond)
y[it1.var_endo + y_kmin * size] = it1.constrained_value;
compute_complete_2b(false, &res1, &res2, &max_res, &max_res_idx);
end_code = it_code;
if (!(isnan(res1) || isinf(res1)))
@ -480,11 +464,9 @@ Interpreter::simulate_a_block(vector_table_conditional_local_type vector_table_c
}
}
if (!cvg)
{
ostringstream tmp;
tmp << " in Solve two boundaries, convergence not achieved in block " << block_num+1 << ", after " << iter << " iterations\n";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" in Solve two boundaries, convergence not achieved in block "
+ to_string(block_num+1) + ", after "
+ to_string(iter) + " iterations\n");
}
else
{
@ -516,9 +498,7 @@ Interpreter::simulate_a_block(vector_table_conditional_local_type vector_table_c
End_Solver();
break;
default:
ostringstream tmp;
tmp << " in simulate_a_block, Unknown type = " << type << "\n";
throw FatalExceptionHandling(tmp.str());
throw FatalExceptionHandling(" in simulate_a_block, Unknown type = " + to_string(type) + "\n");
return ERROR_ON_EXIT;
}
return NO_ERROR_ON_EXIT;
@ -577,52 +557,46 @@ Interpreter::ReadCodeFile(string file_name, CodeLoad &code)
code_liste = code.get_op_code(file_name);
EQN_block_number = code.get_block_number();
if (!code_liste.size())
{
ostringstream tmp;
tmp << " in compute_blocks, " << file_name << ".cod cannot be opened\n";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" in compute_blocks, " + file_name + ".cod cannot be opened\n");
if (block >= static_cast<int>(code.get_block_number()))
{
ostringstream tmp;
tmp << " in compute_blocks, input argument block = " << block+1 << " is greater than the number of blocks in the model (" << code.get_block_number() << " see M_.block_structure_stat.block)\n";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" in compute_blocks, input argument block = " + to_string(block+1)
+ " is greater than the number of blocks in the model ("
+ to_string(code.get_block_number()) + " see M_.block_structure_stat.block)\n");
}
void
Interpreter::check_for_controlled_exo_validity(FBEGINBLOCK_ *fb, vector<s_plan> sconstrained_extended_path)
Interpreter::check_for_controlled_exo_validity(FBEGINBLOCK_ *fb, const vector<s_plan> &sconstrained_extended_path)
{
vector<int> exogenous = fb->get_exogenous();
vector<int> endogenous = fb->get_endogenous();
for (auto & it : sconstrained_extended_path)
{
if ((find(endogenous.begin(), endogenous.end(), it.exo_num) != endogenous.end()) && (find(exogenous.begin(), exogenous.end(), it.var_num) == exogenous.end()))
{
ostringstream tmp;
tmp << "\n the conditional forecast involving as constrained variable " << get_variable(SymbolType::endogenous, it.exo_num) << " and as endogenized exogenous " << get_variable(SymbolType::exogenous, it.var_num) << " that do not appear in block=" << Block_Count+1 << ")\n You should not use block in model options\n";
throw FatalExceptionHandling(tmp.str());
}
else if ((find(endogenous.begin(), endogenous.end(), it.exo_num) != endogenous.end()) && (find(exogenous.begin(), exogenous.end(), it.var_num) != exogenous.end()) && ((fb->get_type() == static_cast<uint8_t>(BlockSimulationType::evaluateForward)) || (fb->get_type() == static_cast<uint8_t>(BlockSimulationType::evaluateBackward))))
{
ostringstream tmp;
tmp << "\n the conditional forecast cannot be implemented for the block=" << Block_Count+1 << ") that has to be evaluated instead to be solved\n You should not use block in model options\n";
throw FatalExceptionHandling(tmp.str());
}
else if (find(previous_block_exogenous.begin(), previous_block_exogenous.end(), it.var_num) != previous_block_exogenous.end())
{
ostringstream tmp;
tmp << "\n the conditional forecast involves in the block " << Block_Count+1 << " the endogenized exogenous " << get_variable(SymbolType::exogenous, it.var_num) << " that appear also in a previous block\n You should not use block in model options\n";
throw FatalExceptionHandling(tmp.str());
}
if (find(endogenous.begin(), endogenous.end(), it.exo_num) != endogenous.end()
&& find(exogenous.begin(), exogenous.end(), it.var_num) == exogenous.end())
throw FatalExceptionHandling("\n the conditional forecast involving as constrained variable "
+ get_variable(SymbolType::endogenous, it.exo_num)
+ " and as endogenized exogenous " + get_variable(SymbolType::exogenous, it.var_num)
+ " that do not appear in block=" + to_string(Block_Count+1)
+ ")\n You should not use block in model options\n");
else if (find(endogenous.begin(), endogenous.end(), it.exo_num) != endogenous.end()
&& find(exogenous.begin(), exogenous.end(), it.var_num) != exogenous.end()
&& (fb->get_type() == static_cast<uint8_t>(BlockSimulationType::evaluateForward)
|| fb->get_type() == static_cast<uint8_t>(BlockSimulationType::evaluateBackward)))
throw FatalExceptionHandling("\n the conditional forecast cannot be implemented for the block="
+ to_string(Block_Count+1) + ") that has to be evaluated instead to be solved\n You should not use block in model options\n");
else if (find(previous_block_exogenous.begin(), previous_block_exogenous.end(), it.var_num)
!= previous_block_exogenous.end())
throw FatalExceptionHandling("\n the conditional forecast involves in the block "
+ to_string(Block_Count+1) + " the endogenized exogenous "
+ get_variable(SymbolType::exogenous, it.var_num)
+ " that appear also in a previous block\n You should not use block in model options\n");
}
for (auto it : exogenous)
previous_block_exogenous.push_back(it);
}
bool
Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int block, bool last_call, bool constrained, vector<s_plan> sconstrained_extended_path, vector_table_conditional_local_type vector_table_conditional_local)
Interpreter::MainLoop(const string &bin_basename, const CodeLoad &code, bool evaluate, int block, bool last_call, bool constrained, const vector<s_plan> &sconstrained_extended_path, const vector_table_conditional_local_type &vector_table_conditional_local)
{
int var;
Block_Count = -1;
@ -729,10 +703,7 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo
delete fb;
}
if (block >= 0)
{
go_on = false;
}
go_on = false;
break;
case Tags::FEND:
@ -752,9 +723,7 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo
T = static_cast<double *>(mxMalloc(var*(periods+y_kmin+y_kmax)*sizeof(double)));
test_mxMalloc(T, __LINE__, __FILE__, __func__, var*(periods+y_kmin+y_kmax)*sizeof(double));
if (block >= 0)
{
it_code = code_liste.begin() + code.get_begin_block(block);
}
it_code = code_liste.begin() + code.get_begin_block(block);
else
it_code++;
break;
@ -769,7 +738,8 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo
{
if (GlobalTemporaryTerms == nullptr)
{
mexPrintf("GlobalTemporaryTerms is NULL\n"); mexEvalString("drawnow;");
mexPrintf("GlobalTemporaryTerms is NULL\n");
mexEvalString("drawnow;");
}
if (var != static_cast<int>(mxGetNumberOfElements(GlobalTemporaryTerms)))
GlobalTemporaryTerms = mxCreateDoubleMatrix(var, 1, mxREAL);
@ -787,9 +757,9 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo
it_code++;
break;
default:
ostringstream tmp;
tmp << " in compute_blocks, unknown command " << static_cast<int>(it_code->first) << " (block=" << Block_Count << ")\n";
throw FatalExceptionHandling(tmp.str());
throw FatalExceptionHandling(" in compute_blocks, unknown command "
+ to_string(static_cast<int>(it_code->first)) + " (block="
+ to_string(Block_Count) + ")\n");
}
}
max_res = max_res_local;
@ -839,16 +809,13 @@ Interpreter::elastic(string str, unsigned int len, bool left)
}
bool
Interpreter::extended_path(string file_name, string bin_basename, bool evaluate, int block, int &nb_blocks, int nb_periods, vector<s_plan> sextended_path, vector<s_plan> sconstrained_extended_path, vector<string> dates, table_conditional_global_type table_conditional_global)
Interpreter::extended_path(const string &file_name, const string &bin_basename, bool evaluate, int block, int &nb_blocks, int nb_periods, const vector<s_plan> &sextended_path, const vector<s_plan> &sconstrained_extended_path, const vector<string> &dates, const table_conditional_global_type &table_conditional_global)
{
CodeLoad code;
ReadCodeFile(file_name, code);
it_code = code_liste.begin();
it_code_type Init_Code = code_liste.begin();
/*size_t size_of_direction = y_size*(periods + y_kmax + y_kmin)*sizeof(double);
double *y_save = (double *) mxMalloc(size_of_direction);
double *x_save = (double *) mxMalloc((periods + y_kmax + y_kmin) * col_x *sizeof(double));*/
size_t size_of_direction = y_size*col_y*sizeof(double);
auto *y_save = static_cast<double *>(mxMalloc(size_of_direction));
test_mxMalloc(y_save, __LINE__, __FILE__, __func__, size_of_direction);
@ -899,14 +866,12 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
mexEvalString("drawnow;");
}
for (const auto & it : sextended_path)
{
x[y_kmin + (it.exo_num - 1) * /*(nb_periods + y_kmax + y_kmin)*/ nb_row_x] = it.value[t];
}
x[y_kmin + (it.exo_num - 1) * nb_row_x] = it.value[t];
it_code = Init_Code;
vector_table_conditional_local.clear();
if (table_conditional_global.size())
vector_table_conditional_local = table_conditional_global[t];
vector_table_conditional_local = table_conditional_global.at(t);
if (t < nb_periods)
MainLoop(bin_basename, code, evaluate, block, false, true, sconstrained_extended_path, vector_table_conditional_local);
else
@ -937,11 +902,6 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
}
}
print_it = old_print_it;
/*for (int j = 0; j < y_size; j++)
{
for(int k = nb_periods; k < periods; k++)
y_save[j + (k + y_kmin) * y_size] = y[ j + ( k - (nb_periods-1) + y_kmin) * y_size];
}*/
for (int i = 0; i < y_size * col_y; i++)
y[i] = y_save[i];
for (int j = 0; j < col_x * nb_row_x; j++)
@ -959,7 +919,7 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
}
bool
Interpreter::compute_blocks(string file_name, string bin_basename, bool evaluate, int block, int &nb_blocks)
Interpreter::compute_blocks(const string &file_name, const string &bin_basename, bool evaluate, int block, int &nb_blocks)
{
CodeLoad code;
ReadCodeFile(file_name, code);

37
mex/sources/bytecode/Interpreter.hh
View File

@ -30,8 +30,6 @@
#include "Evaluate.hh"
#include <dynmex.h>
//#define DEBUGC
using namespace std;
class Interpreter : public dynSparseMatrix
@ -40,11 +38,10 @@ private:
vector<int> previous_block_exogenous;
protected:
void evaluate_a_block(bool initialization);
int simulate_a_block(vector_table_conditional_local_type vector_table_conditional_local);
int simulate_a_block(const vector_table_conditional_local_type &vector_table_conditional_local);
void print_a_block();
string elastic(string str, unsigned int len, bool left);
public:
~Interpreter();
Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
double *direction_arg, size_t y_size_arg,
size_t nb_row_x_arg, size_t nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
@ -52,42 +49,42 @@ public:
string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg,
bool steady_state_arg, bool print_it_arg, int col_x_arg, int col_y_arg);
bool extended_path(string file_name, string bin_basename, bool evaluate, int block, int &nb_blocks, int nb_periods, vector<s_plan> sextended_path, vector<s_plan> sconstrained_extended_path, vector<string> dates, table_conditional_global_type table_conditional_global);
bool compute_blocks(string file_name, string bin_basename, bool evaluate, int block, int &nb_blocks);
void check_for_controlled_exo_validity(FBEGINBLOCK_ *fb, vector<s_plan> sconstrained_extended_path);
bool MainLoop(string bin_basename, CodeLoad code, bool evaluate, int block, bool last_call, bool constrained, vector<s_plan> sconstrained_extended_path, vector_table_conditional_local_type vector_table_conditional_local);
bool extended_path(const string &file_name, const string &bin_basename, bool evaluate, int block, int &nb_blocks, int nb_periods, const vector<s_plan> &sextended_path, const vector<s_plan> &sconstrained_extended_path, const vector<string> &dates, const table_conditional_global_type &table_conditional_global);
bool compute_blocks(const string &file_name, const string &bin_basename, bool evaluate, int block, int &nb_blocks);
void check_for_controlled_exo_validity(FBEGINBLOCK_ *fb, const vector<s_plan> &sconstrained_extended_path);
bool MainLoop(const string &bin_basename, const CodeLoad &code, bool evaluate, int block, bool last_call, bool constrained, const vector<s_plan> &sconstrained_extended_path, const vector_table_conditional_local_type &vector_table_conditional_local);
void ReadCodeFile(string file_name, CodeLoad &code);
inline mxArray *
get_jacob(int block_num)
get_jacob(int block_num) const
{
return jacobian_block[block_num];
};
}
inline mxArray *
get_jacob_exo(int block_num)
get_jacob_exo(int block_num) const
{
return jacobian_exo_block[block_num];
};
}
inline mxArray *
get_jacob_exo_det(int block_num)
get_jacob_exo_det(int block_num) const
{
return jacobian_det_exo_block[block_num];
};
}
inline mxArray *
get_jacob_other_endo(int block_num)
get_jacob_other_endo(int block_num) const
{
return jacobian_other_endo_block[block_num];
};
}
inline vector<double>
get_residual()
get_residual() const
{
return residual;
};
}
inline mxArray *
get_Temporary_Terms()
get_Temporary_Terms() const
{
return GlobalTemporaryTerms;
};
}
};
#endif

22
mex/sources/bytecode/Mem_Mngr.cc
View File

@ -1,5 +1,5 @@
/*
* Copyright © 2007-2017 Dynare Team
* Copyright © 2007-2021 Dynare Team
*
* This file is part of Dynare.
*
@ -24,16 +24,6 @@ Mem_Mngr::Mem_Mngr()
swp_f = false;
swp_f_b = 0;
}
/*void
Mem_Mngr::Print_heap()
{
unsigned int i;
mexPrintf("i :");
for (i = 0; i < CHUNK_SIZE; i++)
mexPrintf("%3d ", i);
mexPrintf("\n");
}
*/
void
Mem_Mngr::init_Mem()
@ -50,7 +40,7 @@ Mem_Mngr::init_Mem()
void
Mem_Mngr::fixe_file_name(string filename_arg)
{
filename_mem = filename_arg;
filename_mem = move(filename_arg);
}
void
@ -67,12 +57,12 @@ Mem_Mngr::mxMalloc_NZE()
{
NonZeroElem *p1 = Chunk_Stack.back();
Chunk_Stack.pop_back();
return (p1);
return p1;
}
else if (CHUNK_heap_pos < CHUNK_SIZE) /*there is enough allocated memory space available we keep it at the top of the heap*/
{
i = CHUNK_heap_pos++;
return (NZE_Mem_add[i]);
return NZE_Mem_add[i];
}
else /*We have to allocate extra memory space*/
{
@ -99,7 +89,7 @@ Mem_Mngr::mxMalloc_NZE()
for (i = CHUNK_heap_pos; i < CHUNK_SIZE; i++)
NZE_Mem_add[i] = const_cast<NonZeroElem *>(NZE_Mem+(i-CHUNK_heap_pos));
i = CHUNK_heap_pos++;
return (NZE_Mem_add[i]);
return NZE_Mem_add[i];
}
}
@ -111,7 +101,7 @@ Mem_Mngr::mxFree_NZE(void *pos)
for (i = 0; i < Nb_CHUNK; i++)
{
gap = (reinterpret_cast<size_t>(pos)-reinterpret_cast<size_t>(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem);
if ((gap < CHUNK_BLCK_SIZE) && (gap >= 0))
if (gap < CHUNK_BLCK_SIZE && gap >= 0)
break;
}
Chunk_Stack.push_back(static_cast<NonZeroElem *>(pos));

2
mex/sources/bytecode/Mem_Mngr.hh
View File

@ -24,7 +24,6 @@
#include <vector>
#include <fstream>
#include <dynmex.h>
//using namespace std;
struct NonZeroElem
{
@ -38,7 +37,6 @@ using v_NonZeroElem = vector<NonZeroElem *>;
class Mem_Mngr
{
public:
//void Print_heap();
void init_Mem();
void mxFree_NZE(void *pos);
NonZeroElem *mxMalloc_NZE();

1405
mex/sources/bytecode/SparseMatrix.cc
View File

File diff suppressed because it is too large Load Diff

107
mex/sources/bytecode/SparseMatrix.hh
View File

@ -19,23 +19,19 @@
#ifndef SPARSEMATRIX_HH_INCLUDED
#define SPARSEMATRIX_HH_INCLUDED
#define PRINTF_ printf
#include <stack>
#include <cmath>
#include <map>
#include <ctime>
#include <tuple>
#include "dynblas.h"
#include "dynumfpack.h"
#include "Mem_Mngr.hh"
#include "ErrorHandling.hh"
//#include "Interpreter.hh"
#include "Evaluate.hh"
#define NEW_ALLOC
#define MARKOVITZ
using namespace std;
struct t_save_op_s
@ -44,44 +40,36 @@ struct t_save_op_s
int first, second;
};
const int IFLD = 0;
const int IFDIV = 1;
const int IFLESS = 2;
const int IFSUB = 3;
const int IFLDZ = 4;
const int IFMUL = 5;
const int IFSTP = 6;
const int IFADD = 7;
const double eps = 1e-15;
const double very_big = 1e24;
const int alt_symbolic_count_max = 1;
const double mem_increasing_factor = 1.1;
constexpr int IFLD = 0, IFDIV = 1, IFLESS = 2, IFSUB = 3, IFLDZ = 4, IFMUL = 5, IFSTP = 6, IFADD = 7;
constexpr double eps = 1e-15, very_big = 1e24;
constexpr int alt_symbolic_count_max = 1;
constexpr double mem_increasing_factor = 1.1;
class dynSparseMatrix : public Evaluate
{
public:
dynSparseMatrix();
dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg, const double slowc_arg);
dynSparseMatrix(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, int periods_arg, int minimal_solving_periods_arg, double slowc_arg);
void Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin, int y_kmax, int Size, int periods, bool cvg, int minimal_solving_periods, int stack_solve_algo, unsigned int endo_name_length, char *P_endo_names, vector_table_conditional_local_type vector_table_conditional_local);
void Simulate_Newton_One_Boundary(bool forward);
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool two_boundaries, int stack_solve_algo, int solve_algo);
void Read_SparseMatrix(const string &file_name, int Size, int periods, int y_kmin, int y_kmax, bool two_boundaries, int stack_solve_algo, int solve_algo);
void Close_SaveCode();
void Read_file(string file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax, int &nb_endo, int &u_count, int &u_count_init, double *u);
void Singular_display(int block, int Size);
void End_Solver();
double g0, gp0, glambda2;
int try_at_iteration;
int find_exo_num(vector<s_plan> sconstrained_extended_path, int value);
int find_int_date(vector<pair<int, double>> per_value, int value);
static int find_exo_num(const vector<s_plan> &sconstrained_extended_path, int value);
static int find_int_date(const vector<pair<int, double>> &per_value, int value);
private:
void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray *x0_m);
void Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, mxArray *x0_m, vector_table_conditional_local_type vector_table_conditional_local, int block_num);
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution, mxArray *x0_m);
void Init_UMFPACK_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, bool &zero_solution, mxArray *x0_m);
void Simple_Init(int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
void Init_GE(int periods, int y_kmin, int y_kmax, int Size, const map<tuple<int, int, int>, int> &IM);
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, const map<tuple<int, int, int>, int> &IM, mxArray *A_m, mxArray *b_m, const mxArray *x0_m) const;
void Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Size, const map<tuple<int, int, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, const mxArray *x0_m, const vector_table_conditional_local_type &vector_table_conditional_local, int block_num) const;
void Init_Matlab_Sparse_Simple(int Size, const map<tuple<int, int, int>, int> &IM, const mxArray *A_m, const mxArray *b_m, bool &zero_solution, const mxArray *x0_m) const;
void Init_UMFPACK_Sparse_Simple(int Size, const map<tuple<int, int, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, bool &zero_solution, const mxArray *x0_m) const;
void Simple_Init(int Size, const map<tuple<int, int, int>, int> &IM, bool &zero_solution);
void End_GE(int Size);
bool mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc);
bool golden(double ax, double bx, double cx, double tol, double solve_tolf, double *xmin);
@ -89,11 +77,11 @@ private:
bool Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, int it_);
void Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
void Print_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, int n);
void Printfull_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n);
void PrintM(int n, double *Ax, mwIndex *Ap, mwIndex *Ai);
static void Print_UMFPack(const SuiteSparse_long *Ap, const SuiteSparse_long *Ai, const double *Ax, int n);
static void Printfull_UMFPack(const SuiteSparse_long *Ap, const SuiteSparse_long *Ai, const double *Ax, const double *b, int n);
static void PrintM(int n, const double *Ax, const mwIndex *Ap, const mwIndex *Ai);
void Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, vector_table_conditional_local_type vector_table_conditional_local);
void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_, const vector_table_conditional_local_type &vector_table_conditional_local);
void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_);
void End_Matlab_LU_UMFPack();
@ -101,52 +89,43 @@ private:
void Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray *x0_m, int precond);
void Check_and_Correct_Previous_Iteration(int block_num, int y_size, int size, double crit_opt_old);
bool Simulate_One_Boundary(int blck, int y_size, int y_kmin, int y_kmax, int Size, bool cvg);
bool solve_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size, const int iter);
void solve_non_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size);
bool solve_linear(int block_num, int y_size, int y_kmin, int y_kmax, int size, int iter);
void solve_non_linear(int block_num, int y_size, int y_kmin, int y_kmax, int size);
string preconditioner_print_out(string s, int preconditioner, bool ss);
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
#ifdef PROFILER
, long int *ndiv, long int *nsub
#endif
);
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long nop4, int Size);
void Grad_f_product(int n, mxArray *b_m, double *vectr, mxArray *A_m, SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b);
void Insert(const int r, const int c, const int u_index, const int lag_index);
void Delete(const int r, const int c);
int At_Row(int r, NonZeroElem **first);
int At_Pos(int r, int c, NonZeroElem **first);
int At_Col(int c, NonZeroElem **first);
int At_Col(int c, int lag, NonZeroElem **first);
int NRow(int r);
int NCol(int c);
int Union_Row(int row1, int row2);
void Print(int Size, int *b);
void Insert(int r, int c, int u_index, int lag_index);
void Delete(int r, int c);
int At_Row(int r, NonZeroElem **first) const;
int At_Pos(int r, int c, NonZeroElem **first) const;
int At_Col(int c, NonZeroElem **first) const;
int At_Col(int c, int lag, NonZeroElem **first) const;
int NRow(int r) const;
int NCol(int c) const;
int Union_Row(int row1, int row2) const;
void Print(int Size, const int *b) const;
int Get_u();
void Delete_u(int pos);
void Clear_u();
void Print_u();
void Print_u() const;
void *Symbolic, *Numeric;
void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods);
void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b);
int complete(int beg_t, int Size, int periods, int *b);
void bksub(int tbreak, int last_period, int Size, double slowc_l
#ifdef PROFILER
, long int *nmul
#endif
);
void bksub(int tbreak, int last_period, int Size, double slowc_l);
void simple_bksub(int it_, int Size, double slowc_l);
mxArray *Sparse_transpose(mxArray *A_m);
mxArray *Sparse_mult_SAT_SB(mxArray *A_m, mxArray *B_m);
mxArray *Sparse_mult_SAT_B(mxArray *A_m, mxArray *B_m);
mxArray *mult_SAT_B(mxArray *A_m, mxArray *B_m);
mxArray *Sparse_substract_SA_SB(mxArray *A_m, mxArray *B_m);
mxArray *Sparse_substract_A_SB(mxArray *A_m, mxArray *B_m);
mxArray *substract_A_B(mxArray *A_m, mxArray *B_m);
static mxArray *Sparse_transpose(const mxArray *A_m);
static mxArray *Sparse_mult_SAT_SB(const mxArray *A_m, const mxArray *B_m);
static mxArray *Sparse_mult_SAT_B(const mxArray *A_m, const mxArray *B_m);
static mxArray *mult_SAT_B(const mxArray *A_m, const mxArray *B_m);
static mxArray *Sparse_substract_SA_SB(const mxArray *A_m, const mxArray *B_m);
static mxArray *Sparse_substract_A_SB(const mxArray *A_m, const mxArray *B_m);
static mxArray *substract_A_B(const mxArray *A_m, const mxArray *B_m);
protected:
stack<double> Stack;
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
int middle_count_loop;
//char type;
fstream SaveCode;
string filename;
int max_u, min_u;
@ -171,7 +150,7 @@ protected:
double markowitz_c_s;
double res1a;
long int nop_all, nop1, nop2;
map<pair<pair<int, int>, int>, int> IM_i;
map<tuple<int, int, int>, int> IM_i;
protected:
vector<double> residual;
int u_count_alloc, u_count_alloc_save;
@ -186,7 +165,7 @@ protected:
int restart;
double g_lambda1, g_lambda2, gp_0;
double lu_inc_tol;
//private:
SuiteSparse_long *Ap_save, *Ai_save;
double *Ax_save, *b_save;
mxArray *A_m_save, *b_m_save;

188
mex/sources/bytecode/bytecode.cc
View File

@ -95,7 +95,7 @@ Get_Arguments_and_global_variables(int nrhs,
steady_col_y = mxGetN(prhs[i]);
break;
case 4:
periods = int (mxGetScalar(prhs[i]));
periods = static_cast<int>(mxGetScalar(prhs[i]));
break;
case 5:
*block_structur = mxDuplicateArray(prhs[i]);
@ -165,11 +165,7 @@ Get_Arguments_and_global_variables(int nrhs,
*plan_struct_name = deblank(Get_Argument(prhs[i]).substr(pos, string::npos));
}
else
{
ostringstream tmp;
tmp << " in main, unknown argument : " << Get_Argument(prhs[i]) << "\n";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" in main, unknown argument : " + Get_Argument(prhs[i]) + "\n");
}
}
if (count_array_argument > 0 && count_array_argument < 5)
@ -177,34 +173,20 @@ Get_Arguments_and_global_variables(int nrhs,
if (count_array_argument == 3 && steady_state)
periods = 1;
else
{
ostringstream tmp;
tmp << " in main, missing arguments. All the following arguments have to be indicated y, x, params, it_, ys\n";
throw FatalExceptionHandling(tmp.str());
}
throw FatalExceptionHandling(" in main, missing arguments. All the following arguments have to be indicated y, x, params, it_, ys\n");
}
*M_ = mexGetVariable("global", "M_");
if (*M_ == nullptr)
{
ostringstream tmp;
tmp << " in main, global variable not found: M_\n";
throw FatalExceptionHandling(tmp.str());
}
if (!*M_)
throw FatalExceptionHandling(" in main, global variable not found: M_\n");
/* Gets variables and parameters from global workspace of Matlab */
*oo_ = mexGetVariable("global", "oo_");
if (*oo_ == nullptr)
{
ostringstream tmp;
tmp << " in main, global variable not found: oo_\n";
throw FatalExceptionHandling(tmp.str());
}
if (!*oo_)
throw FatalExceptionHandling(" in main, global variable not found: oo_\n");
*options_ = mexGetVariable("global", "options_");
if (*options_ == nullptr)
{
ostringstream tmp;
tmp << " in main, global variable not found: options_\n";
throw FatalExceptionHandling(tmp.str());
}
if (!*options_)
throw FatalExceptionHandling(" in main, global variable not found: options_\n");
}
/* The gateway routine */
@ -352,9 +334,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
controlled_varexo = mxGetField(options_cond_fcst_, 0, "controlled_varexo");
nb_controlled = mxGetM(controlled_varexo) * mxGetN(controlled_varexo);
if (nb_controlled != nb_constrained)
{
mexErrMsgTxt("The number of exogenized variables and the number of exogenous controlled variables should be equal.");
}
mexErrMsgTxt("The number of exogenized variables and the number of exogenous controlled variables should be equal.");
}
double *controlled_varexo_value = nullptr;
if (controlled_varexo != nullptr)
@ -392,34 +372,25 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int *constrained_int_date = static_cast<int *>(mxMalloc(nb_local_periods * sizeof(int)));
error_msg.test_mxMalloc(constrained_int_date, __LINE__, __FILE__, __func__, nb_local_periods * sizeof(int));
if (nb_periods < nb_local_periods)
{
ostringstream oss;
oss << nb_periods;
string tmp = oss.str();
tmp.insert(0, "The total number of simulation periods (");
tmp.append(") is lesser than the number of periods in the shock definitions (");
oss << nb_local_periods;
string tmp1 = oss.str();
tmp.append(tmp1);
tmp.append(")");
mexErrMsgTxt(tmp.c_str());
}
(sconditional_extended_path[i]).per_value.resize(nb_local_periods);
(sconditional_extended_path[i]).value.resize(nb_periods);
mexErrMsgTxt((string{"The total number of simulation periods ("} + to_string(nb_periods)
+ ") is lesser than the number of periods in the shock definitions ("
+ to_string(nb_local_periods)).c_str());
sconditional_extended_path[i].per_value.resize(nb_local_periods);
sconditional_extended_path[i].value.resize(nb_periods);
for (int j = 0; j < nb_periods; j++)
sconditional_extended_path[i].value[j] = 0;
for (int j = 0; j < nb_local_periods; j++)
{
constrained_int_date[j] = int (specific_constrained_int_date_[j]) - 1;
constrained_int_date[j] = static_cast<int>(specific_constrained_int_date_[j]) - 1;
conditional_local.is_cond = true;
conditional_local.var_exo = sconditional_extended_path[i].var_num - 1;
conditional_local.var_endo = sconditional_extended_path[i].exo_num - 1;
conditional_local.constrained_value = specific_constrained_paths_[j];
table_conditional_global[constrained_int_date[j]][sconditional_extended_path[i].exo_num - 1] = conditional_local;
sconditional_extended_path[i].per_value[j] = make_pair(constrained_int_date[j], specific_constrained_paths_[j]);
sconditional_extended_path[i].per_value[j] = { constrained_int_date[j], specific_constrained_paths_[j] };
sconditional_extended_path[i].value[constrained_int_date[j]] = specific_constrained_paths_[j];
if (max_periods < constrained_int_date[j] + 1)
max_periods = constrained_int_date[j] + 1;
max_periods = max(max_periods, constrained_int_date[j] + 1);
}
mxFree(constrained_int_date);
}
@ -435,28 +406,18 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
double *specific_shock_int_date_ = mxGetPr(mxGetCell(shock_int_date_, i));
int nb_local_periods = mxGetM(Array_shock_paths_) * mxGetN(Array_shock_paths_);
if (nb_periods < nb_local_periods)
{
ostringstream oss;
oss << nb_periods;
string tmp = oss.str();
tmp.insert(0, "The total number of simulation periods (");
tmp.append(") is lesser than the number of periods in the shock definitions (");
oss << nb_local_periods;
string tmp1 = oss.str();
tmp.append(tmp1);
tmp.append(")");
mexErrMsgTxt(tmp.c_str());
}
(sextended_path[i]).per_value.resize(nb_local_periods);
(sextended_path[i]).value.resize(nb_periods);
mexErrMsgTxt((string{"The total number of simulation periods ("} + to_string(nb_periods)
+ ") is lesser than the number of periods in the shock definitions ("
+ to_string(nb_local_periods)).c_str());
sextended_path[i].per_value.resize(nb_local_periods);
sextended_path[i].value.resize(nb_periods);
for (int j = 0; j < nb_periods; j++)
sextended_path[i].value[j] = 0;
for (int j = 0; j < nb_local_periods; j++)
{
sextended_path[i].per_value[j] = make_pair(int (specific_shock_int_date_[j]), specific_shock_paths_[j]);
sextended_path[i].value[int (specific_shock_int_date_[j]-1)] = specific_shock_paths_[j];
if (max_periods < int (specific_shock_int_date_[j]))
max_periods = int (specific_shock_int_date_[j]);
sextended_path[i].per_value[j] = { static_cast<int>(specific_shock_int_date_[j]), specific_shock_paths_[j] };
sextended_path[i].value[static_cast<int>(specific_shock_int_date_[j]-1)] = specific_shock_paths_[j];
max_periods = max(max_periods, static_cast<int>(specific_shock_int_date_[j]));
}
}
for (int i = 0; i < nb_periods; i++)
@ -528,10 +489,10 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
if (tmp)
{
size_t num_shocks = mxGetM(tmp);
(splan[i]).per_value.resize(num_shocks);
splan[i].per_value.resize(num_shocks);
double *per_value = mxGetPr(tmp);
for (int j = 0; j < static_cast<int>(num_shocks); j++)
(splan[i]).per_value[j] = make_pair(ceil(per_value[j]), per_value[j + num_shocks]);
splan[i].per_value[j] = { ceil(per_value[j]), per_value[j + num_shocks] };
}
}
int i = 0;
@ -543,10 +504,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexPrintf(" plan fliping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num, it.exo.c_str(), it.exo_num);
else
mexPrintf(" plan shocks on var=%s for the following periods and with the following values:\n", it.var.c_str());
for (auto & it1 : it.per_value)
{
mexPrintf(" %3d %10.5f\n", it1.first, it1.second);
}
for (auto &[period, value]: it.per_value)
mexPrintf(" %3d %10.5f\n", period, value);
i++;
}
}
@ -555,11 +514,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
pfplan_struct = mexGetVariable("base", pfplan.c_str());
if (!pfplan_struct)
{
string tmp = pfplan;
tmp.insert(0, "Can't find the pfplan: ");
mexErrMsgTxt(tmp.c_str());
}
mexErrMsgTxt(("Can't find the pfplan: " + pfplan).c_str());
size_t n_plan = mxGetN(pfplan_struct);
spfplan.resize(n_plan);
for (int i = 0; i < static_cast<int>(n_plan); i++)
@ -607,9 +562,9 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
size_t num_shocks = mxGetM(tmp);
double *per_value = mxGetPr(tmp);
(spfplan[i]).per_value.resize(num_shocks);
spfplan[i].per_value.resize(num_shocks);
for (int j = 0; j < static_cast<int>(num_shocks); j++)
spfplan[i].per_value[j] = make_pair(ceil(per_value[j]), per_value[j+ num_shocks]);
spfplan[i].per_value[j] = { ceil(per_value[j]), per_value[j+ num_shocks] };
}
}
int i = 0;
@ -621,10 +576,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexPrintf(" plan flipping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num, it.exo.c_str(), it.exo_num);
else
mexPrintf(" plan shocks on var=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num);
for (auto & it1 : it.per_value)
{
mexPrintf(" %3d %10.5f\n", it1.first, it1.second);
}
for (auto &[period, value] : it.per_value)
mexPrintf(" %3d %10.5f\n", period, value);
i++;
}
}
@ -660,17 +613,17 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
}
int field = mxGetFieldNumber(M_, "maximum_lag");
if (field >= 0)
y_kmin = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
y_kmin = static_cast<int>(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
else
mexErrMsgTxt("maximum_lag is not a field of M_");
field = mxGetFieldNumber(M_, "maximum_lead");
if (field >= 0)
y_kmax = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
y_kmax = static_cast<int>(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
else
mexErrMsgTxt("maximum_lead is not a field of M_");
field = mxGetFieldNumber(M_, "maximum_endo_lag");
if (field >= 0)
y_decal = max(0, y_kmin-int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field))))));
y_decal = max(0, y_kmin-static_cast<int>(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field))))));
else
mexErrMsgTxt("maximum_endo_lag is not a field of M_");
@ -678,7 +631,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
int field = mxGetFieldNumber(options_, "periods");
if (field >= 0)
periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, field)))));
periods = static_cast<int>(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, field)))));
else
mexErrMsgTxt("options_ is not a field of options_");
}
@ -711,7 +664,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int field = mxGetFieldNumber(options_, "verbosity");
int verbose = 0;
if (field >= 0)
verbose = int (*mxGetPr((mxGetFieldByNumber(options_, 0, field))));
verbose = static_cast<int>(*mxGetPr((mxGetFieldByNumber(options_, 0, field))));
else
mexErrMsgTxt("verbosity is not a field of options_");
if (verbose)
@ -738,23 +691,23 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
else
mexErrMsgTxt("maxit is not a field of options_.steady");
}
int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(temporaryfield, 0, field)))));
int maxit_ = static_cast<int>(floor(*mxGetPr(mxGetFieldByNumber(temporaryfield, 0, field))));
field = mxGetFieldNumber(options_, "slowc");
if (field < 0)
mexErrMsgTxt("slows is not a field of options_");
auto slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
auto slowc = static_cast<double>(*mxGetPr(mxGetFieldByNumber(options_, 0, field)));
field = mxGetFieldNumber(options_, "markowitz");
if (field < 0)
mexErrMsgTxt("markowitz is not a field of options_");
auto markowitz_c = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
auto markowitz_c = static_cast<double>(*mxGetPr(mxGetFieldByNumber(options_, 0, field)));
field = mxGetFieldNumber(options_, "minimal_solving_periods");
if (field < 0)
mexErrMsgTxt("minimal_solving_periods is not a field of options_");
int minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
int minimal_solving_periods = static_cast<int>(*mxGetPr(mxGetFieldByNumber(options_, 0, field)));
field = mxGetFieldNumber(options_, "stack_solve_algo");
if (field < 0)
mexErrMsgTxt("stack_solve_algo is not a field of options_");
int stack_solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
int stack_solve_algo = static_cast<int>(*mxGetPr(mxGetFieldByNumber(options_, 0, field)));
int solve_algo;
double solve_tolf;
@ -762,12 +715,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
int field = mxGetFieldNumber(options_, "solve_algo");
if (field >= 0)
solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
solve_algo = static_cast<int>(*mxGetPr(mxGetFieldByNumber(options_, 0, field)));
else
mexErrMsgTxt("solve_algo is not a field of options_");
field = mxGetFieldNumber(options_, "solve_tolf");
if (field >= 0)
solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, field)));
solve_tolf = *mxGetPr(mxGetFieldByNumber(options_, 0, field));
else
mexErrMsgTxt("solve_tolf is not a field of options_");
}
@ -782,7 +735,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexErrMsgTxt("dynatol is not a field of options_");
field = mxGetFieldNumber(dynatol, "f");
if (field >= 0)
solve_tolf = *mxGetPr((mxGetFieldByNumber(dynatol, 0, field)));
solve_tolf = *mxGetPr(mxGetFieldByNumber(dynatol, 0, field));
else
mexErrMsgTxt("f is not a field of options_.dynatol");
}
@ -793,9 +746,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
else
mexErrMsgTxt("fname is not a field of M_");
size_t buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
char *fname;
fname = static_cast<char *>(mxCalloc(buflen+1, sizeof(char)));
size_t status = mxGetString(mxa, fname, int (buflen));
char *fname = static_cast<char *>(mxCalloc(buflen+1, sizeof(char)));
size_t status = mxGetString(mxa, fname, static_cast<int>(buflen));
fname[buflen] = ' ';
if (status != 0)
mexWarnMsgTxt("Not enough space. Filename is truncated.");
@ -812,17 +764,14 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
direction = static_cast<double *>(mxMalloc(size_of_direction));
error_msg.test_mxMalloc(direction, __LINE__, __FILE__, __func__, size_of_direction);
memset(direction, 0, size_of_direction);
/*mexPrintf("col_x : %d, row_x : %d\n",col_x, row_x);*/
auto *x = static_cast<double *>(mxMalloc(col_x*row_x*sizeof(double)));
error_msg.test_mxMalloc(x, __LINE__, __FILE__, __func__, col_x*row_x*sizeof(double));
for (i = 0; i < row_x*col_x; i++)
{
x[i] = double (xd[i]);
}
x[i] = static_cast<double>(xd[i]);
for (i = 0; i < row_y*col_y; i++)
{
y[i] = double (yd[i]);
ya[i] = double (yd[i]);
y[i] = static_cast<double>(yd[i]);
ya[i] = static_cast<double>(yd[i]);
}
size_t y_size = row_y;
size_t nb_row_x = row_x;
@ -861,7 +810,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
clock_t t1 = clock();
if (!steady_state && !evaluate && print)
mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
mexPrintf("Simulation Time=%f milliseconds\n",
1000.0*(static_cast<double>(t1)-static_cast<double>(t0))/static_cast<double>(CLOCKS_PER_SEC));
bool dont_store_a_structure = false;
if (nlhs > 0)
{
@ -870,7 +820,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
if (evaluate)
{
vector<double> residual = interprete.get_residual();
plhs[0] = mxCreateDoubleMatrix(int (residual.size()/double (col_y)), int (col_y), mxREAL);
plhs[0] = mxCreateDoubleMatrix(static_cast<int>(residual.size()/static_cast<double>(col_y)),
static_cast<int>(col_y), mxREAL);
pind = mxGetPr(plhs[0]);
for (i = 0; i < residual.size(); i++)
pind[i] = residual[i];
@ -882,7 +833,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
out_periods = max_periods + y_kmin;
else
out_periods = row_y;
plhs[0] = mxCreateDoubleMatrix(out_periods, int (col_y), mxREAL);
plhs[0] = mxCreateDoubleMatrix(out_periods, static_cast<int>(col_y), mxREAL);
pind = mxGetPr(plhs[0]);
for (i = 0; i < out_periods*col_y; i++)
pind[i] = y[i];
@ -895,7 +846,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
out_periods = max_periods + y_kmin;
else
out_periods = col_y;
plhs[0] = mxCreateDoubleMatrix(int (row_y), out_periods, mxREAL);
plhs[0] = mxCreateDoubleMatrix(static_cast<int>(row_y), out_periods, mxREAL);
pind = mxGetPr(plhs[0]);
if (evaluate)
{
@ -919,13 +870,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
jacob_exo_field_number = 1;
jacob_exo_det_field_number = 2;
jacob_other_endo_field_number = 3;
mwSize dims[1] = {static_cast<mwSize>(nb_blocks) };
mwSize dims[1] = { static_cast<mwSize>(nb_blocks) };
plhs[1] = mxCreateStructArray(1, dims, 4, field_names);
}
else if (!mxIsStruct(block_structur))
{
plhs[1] = interprete.get_jacob(0);
//mexCallMATLAB(0,NULL, 1, &plhs[1], "disp");
dont_store_a_structure = true;
}
else
@ -960,17 +910,14 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
}
else
{
plhs[1] = mxCreateDoubleMatrix(int (row_x), int (col_x), mxREAL);
plhs[1] = mxCreateDoubleMatrix(static_cast<int>(row_x), static_cast<int>(col_x), mxREAL);
pind = mxGetPr(plhs[1]);
for (i = 0; i < row_x*col_x; i++)
{
pind[i] = x[i];
}
pind[i] = x[i];
}
if (nlhs > 2)
{
plhs[2] = mxCreateDoubleMatrix(int (row_y), int (col_y), mxREAL);
plhs[2] = mxCreateDoubleMatrix(static_cast<int>(row_y), static_cast<int>(col_y), mxREAL);
pind = mxGetPr(plhs[2]);
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i];
@ -978,7 +925,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mxArray *GlobalTemporaryTerms = interprete.get_Temporary_Terms();
size_t nb_temp_terms = mxGetM(GlobalTemporaryTerms);
plhs[3] = mxCreateDoubleMatrix(int (nb_temp_terms), 1, mxREAL);
plhs[3] = mxCreateDoubleMatrix(static_cast<int>(nb_temp_terms), 1, mxREAL);
pind = mxGetPr(plhs[3]);
double *tt = mxGetPr(GlobalTemporaryTerms);
for (i = 0; i < nb_temp_terms; i++)
@ -996,5 +943,4 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mxFree(ya);
if (direction)
mxFree(direction);
return;
}

2
preprocessor

@ -1 +1 @@
Subproject commit 07391e0f9fdc8e697131cd3e292d6c9f23afd63b
Subproject commit 702cf62ee386b6d240d9d4fcbb0ddeb33bb35c2f