From 4893f0e82cacae7ed0e9a7ca8d02c21c73cffcf5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 3 Feb 2021 18:10:01 +0100
Subject: [PATCH] Bytecode: various simplifications and modernizations
---
mex/sources/bytecode/ErrorHandling.hh | 124 +--
mex/sources/bytecode/Evaluate.cc | 228 ++--
mex/sources/bytecode/Evaluate.hh | 24 +-
mex/sources/bytecode/Interpreter.cc | 156 +--
mex/sources/bytecode/Interpreter.hh | 37 +-
mex/sources/bytecode/Mem_Mngr.cc | 22 +-
mex/sources/bytecode/Mem_Mngr.hh | 2 -
mex/sources/bytecode/SparseMatrix.cc | 1405 +++++++++----------------
mex/sources/bytecode/SparseMatrix.hh | 107 +-
mex/sources/bytecode/bytecode.cc | 188 ++--
preprocessor | 2 +-
11 files changed, 830 insertions(+), 1465 deletions(-)
diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh
index b376173d3e..6a019cff7b 100644
--- a/mex/sources/bytecode/ErrorHandling.hh
+++ b/mex/sources/bytecode/ErrorHandling.hh
@@ -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
diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc
index 4b93203970..bf3d803911 100644
--- a/mex/sources/bytecode/Evaluate.cc
+++ b/mex/sources/bytecode/Evaluate.cc
@@ -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;
diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh
index 40914e89cb..9526042e46 100644
--- a/mex/sources/bytecode/Evaluate.hh
+++ b/mex/sources/bytecode/Evaluate.hh
@@ -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);
};
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index c62d68c6f6..5645c3d745 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -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);
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index ece879401a..08c5944750 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -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
diff --git a/mex/sources/bytecode/Mem_Mngr.cc b/mex/sources/bytecode/Mem_Mngr.cc
index e706bd92cd..0446e069ba 100644
--- a/mex/sources/bytecode/Mem_Mngr.cc
+++ b/mex/sources/bytecode/Mem_Mngr.cc
@@ -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));
diff --git a/mex/sources/bytecode/Mem_Mngr.hh b/mex/sources/bytecode/Mem_Mngr.hh
index 3124f526b6..d02e967229 100644
--- a/mex/sources/bytecode/Mem_Mngr.hh
+++ b/mex/sources/bytecode/Mem_Mngr.hh
@@ -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();
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 7009f1d3d0..06b083ddcb 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -17,14 +17,11 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
-//define _GLIBCXX_USE_C99_FENV_TR1 1
-//include <cfenv>
-
#include <cstring>
#include <ctime>
#include <sstream>
-//#include <gsl/gsl_min.h>
-//#include <minimize.h>
+#include <algorithm>
+
#include "SparseMatrix.hh"
using namespace std;
@@ -50,8 +47,8 @@ dynSparseMatrix::dynSparseMatrix()
Numeric = nullptr;
}
-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::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) :
Evaluate(y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, periods_arg, minimal_solving_periods_arg, slowc_arg)
{
pivotva = nullptr;
@@ -74,26 +71,26 @@ dynSparseMatrix::dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, con
}
int
-dynSparseMatrix::NRow(int r)
+dynSparseMatrix::NRow(int r) const
{
return NbNZRow[r];
}
int
-dynSparseMatrix::NCol(int c)
+dynSparseMatrix::NCol(int c) const
{
return NbNZCol[c];
}
int
-dynSparseMatrix::At_Row(int r, NonZeroElem **first)
+dynSparseMatrix::At_Row(int r, NonZeroElem **first) const
{
- (*first) = FNZE_R[r];
+ *first = FNZE_R[r];
return NbNZRow[r];
}
int
-dynSparseMatrix::Union_Row(int row1, int row2)
+dynSparseMatrix::Union_Row(int row1, int row2) const
{
NonZeroElem *first1, *first2;
int n1 = At_Row(row1, &first1);
@@ -126,31 +123,31 @@ dynSparseMatrix::Union_Row(int row1, int row2)
}
int
-dynSparseMatrix::At_Pos(int r, int c, NonZeroElem **first)
+dynSparseMatrix::At_Pos(int r, int c, NonZeroElem **first) const
{
- (*first) = FNZE_R[r];
+ *first = FNZE_R[r];
while ((*first)->c_index != c)
- (*first) = (*first)->NZE_R_N;
+ *first = (*first)->NZE_R_N;
return NbNZRow[r];
}
int
-dynSparseMatrix::At_Col(int c, NonZeroElem **first)
+dynSparseMatrix::At_Col(int c, NonZeroElem **first) const
{
- (*first) = FNZE_C[c];
+ *first = FNZE_C[c];
return NbNZCol[c];
}
int
-dynSparseMatrix::At_Col(int c, int lag, NonZeroElem **first)
+dynSparseMatrix::At_Col(int c, int lag, NonZeroElem **first) const
{
- (*first) = FNZE_C[c];
+ *first = FNZE_C[c];
int i = 0;
- while ((*first)->lag_index != lag && (*first))
- (*first) = (*first)->NZE_C_N;
- if ((*first))
+ while ((*first)->lag_index != lag && *first)
+ *first = (*first)->NZE_C_N;
+ if (*first)
{
- NonZeroElem *firsta = (*first);
+ NonZeroElem *firsta = *first;
if (!firsta->NZE_C_N)
i++;
else
@@ -168,7 +165,7 @@ dynSparseMatrix::At_Col(int c, int lag, NonZeroElem **first)
}
void
-dynSparseMatrix::Delete(const int r, const int c)
+dynSparseMatrix::Delete(int r, int c)
{
NonZeroElem *first = FNZE_R[r], *firsta = nullptr;
@@ -202,31 +199,30 @@ dynSparseMatrix::Delete(const int r, const int c)
}
void
-dynSparseMatrix::Print(int Size, int *b)
+dynSparseMatrix::Print(int Size, const int *b) const
{
- int a, i, j, k, l;
mexPrintf(" ");
- for (k = 0; k < Size*periods; k++)
+ for (int k = 0; k < Size*periods; k++)
mexPrintf("%-2d ", k);
mexPrintf(" | ");
- for (k = 0; k < Size*periods; k++)
+ for (int k = 0; k < Size*periods; k++)
mexPrintf("%8d", k);
mexPrintf("\n");
- for (i = 0; i < Size*periods; i++)
+ for (int i = 0; i < Size*periods; i++)
{
NonZeroElem *first = FNZE_R[i];
- j = NbNZRow[i];
+ int j = NbNZRow[i];
mexPrintf("%-2d ", i);
- a = 0;
- for (k = 0; k < j; k++)
+ int a = 0;
+ for (int k = 0; k < j; k++)
{
- for (l = 0; l < (first->c_index-a); l++)
+ for (int l = 0; l < (first->c_index-a); l++)
mexPrintf(" ");
mexPrintf("%-2d ", first->u_index);
a = first->c_index+1;
first = first->NZE_R_N;
}
- for (k = a; k < Size*periods; k++)
+ for (int k = a; k < Size*periods; k++)
mexPrintf(" ");
mexPrintf("%-2d ", b[i]);
@@ -234,23 +230,23 @@ dynSparseMatrix::Print(int Size, int *b)
j = NbNZRow[i];
mexPrintf(" | %-2d ", i);
a = 0;
- for (k = 0; k < j; k++)
+ for (int k = 0; k < j; k++)
{
- for (l = 0; l < (first->c_index-a); l++)
+ for (int l = 0; l < (first->c_index-a); l++)
mexPrintf(" ");
- mexPrintf("%8.4f", double (u[first->u_index]));
+ mexPrintf("%8.4f", static_cast<double>(u[first->u_index]));
a = first->c_index+1;
first = first->NZE_R_N;
}
- for (k = a; k < Size*periods; k++)
+ for (int k = a; k < Size*periods; k++)
mexPrintf(" ");
- mexPrintf("%8.4f", double (u[b[i]]));
+ mexPrintf("%8.4f", static_cast<double>(u[b[i]]));
mexPrintf("\n");
}
}
void
-dynSparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_index)
+dynSparseMatrix::Insert(int r, int c, int u_index, int lag_index)
{
NonZeroElem *firstn, *first, *firsta, *a;
firstn = mem_mngr.mxMalloc_NZE();
@@ -310,13 +306,11 @@ dynSparseMatrix::Close_SaveCode()
}
void
-dynSparseMatrix::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)
+dynSparseMatrix::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)
{
unsigned int eq, var;
int lag;
mem_mngr.fixe_file_name(file_name);
- /*mexPrintf("steady_state=%d, size=%d, solve_algo=%d, stack_solve_algo=%d, two_boundaries=%d\n",steady_state, Size, solve_algo, stack_solve_algo, two_boundaries);
- mexEvalString("drawnow;");*/
if (!SaveCode.is_open())
{
if (steady_state)
@@ -345,10 +339,10 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&val), sizeof(val));
- IM_i[make_pair(make_pair(eq, var), lag)] = val;
+ IM_i[{ eq, var, lag }] = val;
}
for (int j = 0; j < Size; j++)
- IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
+ IM_i[{ j, Size*(periods+y_kmax), 0 }] = j;
}
else if (stack_solve_algo >= 0 && stack_solve_algo <= 4)
{
@@ -359,10 +353,10 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&val), sizeof(val));
- IM_i[make_pair(make_pair(var - lag*Size, -lag), eq)] = val;
+ IM_i[{ var - lag*Size, -lag, eq }] = val;
}
for (int j = 0; j < Size; j++)
- IM_i[make_pair(make_pair(Size*(periods+y_kmax), 0), j)] = j;
+ IM_i[{ Size*(periods+y_kmax), 0, j }] = j;
}
else if (stack_solve_algo == 7)
{
@@ -373,10 +367,10 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&val), sizeof(val));
- IM_i[make_pair(make_pair(eq, lag), var - lag * Size)] = val;
+ IM_i[{ eq, lag, var - lag * Size }] = val;
}
for (int j = 0; j < Size; j++)
- IM_i[make_pair(make_pair(Size*(periods+y_kmax), 0), j)] = j;
+ IM_i[{ Size*(periods+y_kmax), 0, j }] = j;
}
}
else
@@ -390,10 +384,11 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&val), sizeof(val));
- IM_i[make_pair(make_pair(eq, var), lag)] = val;
+ IM_i[{ eq, var, lag }] = val;
}
}
- else if (((stack_solve_algo >= 0 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 6 || solve_algo <= 8) && steady_state))
+ else if (((stack_solve_algo >= 0 || stack_solve_algo <= 4) && !steady_state)
+ || ((solve_algo >= 6 || solve_algo <= 8) && steady_state))
{
for (int i = 0; i < u_count_init; i++)
{
@@ -402,7 +397,7 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&val), sizeof(val));
- IM_i[make_pair(make_pair(var - lag*Size, -lag), eq)] = val;
+ IM_i[{ var - lag*Size, -lag, eq }] = val;
}
}
}
@@ -412,10 +407,8 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods
SaveCode.read(reinterpret_cast<char *>(&index_vara[j]), sizeof(*index_vara));
if (periods+y_kmin+y_kmax > 1)
for (int i = 1; i < periods+y_kmin+y_kmax; i++)
- {
- for (int j = 0; j < Size; j++)
- index_vara[j+Size*i] = index_vara[j+Size*(i-1)] + y_size;
- }
+ for (int j = 0; j < Size; j++)
+ index_vara[j+Size*i] = index_vara[j+Size*(i-1)] + y_size;
index_equa = static_cast<int *>(mxMalloc(Size*sizeof(int)));
test_mxMalloc(index_equa, __LINE__, __FILE__, __func__, Size*sizeof(int));
for (int j = 0; j < Size; j++)
@@ -423,11 +416,8 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods
}
void
-dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM, bool &zero_solution)
+dynSparseMatrix::Simple_Init(int Size, const map<tuple<int, int, int>, int> &IM, bool &zero_solution)
{
- int i, eq, var, lag;
- map<pair<pair<int, int>, int>, int>::iterator it4;
- NonZeroElem *first;
pivot = static_cast<int *>(mxMalloc(Size*sizeof(int)));
test_mxMalloc(pivot, __LINE__, __FILE__, __func__, Size*sizeof(int));
pivot_save = static_cast<int *>(mxMalloc(Size*sizeof(int)));
@@ -446,7 +436,7 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
mem_mngr.init_CHUNK_BLCK_SIZE(u_count);
g_save_op = nullptr;
g_nop_all = 0;
- i = Size*sizeof(NonZeroElem *);
+ int i = Size*sizeof(NonZeroElem *);
FNZE_R = static_cast<NonZeroElem **>(mxMalloc(i));
test_mxMalloc(FNZE_R, __LINE__, __FILE__, __func__, i);
FNZE_C = static_cast<NonZeroElem **>(mxMalloc(i));
@@ -460,8 +450,6 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
test_mxMalloc(NbNZRow, __LINE__, __FILE__, __func__, i);
NbNZCol = static_cast<int *>(mxMalloc(i));
test_mxMalloc(NbNZCol, __LINE__, __FILE__, __func__, i);
- it4 = IM.begin();
- eq = -1;
for (i = 0; i < Size; i++)
{
line_done[i] = false;
@@ -473,16 +461,14 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
NbNZCol[i] = 0;
}
int u_count1 = Size;
- while (it4 != IM.end())
+ for (auto &[key, value] : IM)
{
- var = it4->first.first.second;
- eq = it4->first.first.first;
- lag = it4->first.second;
+ auto &[eq, var, lag] = key;
if (lag == 0) /*Build the index for sparse matrix containing the jacobian : u*/
{
NbNZRow[eq]++;
NbNZCol[var]++;
- first = mem_mngr.mxMalloc_NZE();
+ NonZeroElem *first = mem_mngr.mxMalloc_NZE();
first->NZE_C_N = nullptr;
first->NZE_R_N = nullptr;
first->u_index = u_count1;
@@ -501,7 +487,6 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
temp_NZE_C[var] = first;
u_count1++;
}
- it4++;
}
double cum_abs_sum = 0;
for (int i = 0; i < Size; i++)
@@ -520,45 +505,24 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
}
void
-dynSparseMatrix::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)
+dynSparseMatrix::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
{
- int eq, var;
double *b = mxGetPr(b_m);
if (!b)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, can't retrieve b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, can't retrieve b vector\n");
double *x0 = mxGetPr(x0_m);
if (!x0)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, can't retrieve x0 vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, can't retrieve x0 vector\n");
mwIndex *Ai = mxGetIr(A_m);
if (!Ai)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, can't allocate Ai index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, can't allocate Ai index vector\n");
mwIndex *Aj = mxGetJc(A_m);
if (!Aj)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, can't allocate Aj index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, can't allocate Aj index vector\n");
double *A = mxGetPr(A_m);
if (!A)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, can't retrieve A matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
- map<pair<pair<int, int>, int>, int>::iterator it4;
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, can't retrieve A matrix\n");
+
for (int i = 0; i < y_size*(periods+y_kmin); i++)
ya[i] = y[i];
#ifdef DEBUG
@@ -580,104 +544,66 @@ dynSparseMatrix::Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, in
Aj[0] = 0;
last_var = 0;
- it4 = IM.begin();
- while (it4 != IM.end())
+ for (auto &[key, index] : IM)
{
- var = it4->first.first.first;
+ auto &[var, ignore, eq] = key;
if (var != last_var)
{
Aj[1+last_var] = NZE;
last_var = var;
}
- eq = it4->first.second;
- int index = it4->second;
#ifdef DEBUG
if (index < 0 || index >= u_count_alloc || index > Size + Size*Size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, index (" << index << ") out of range for u vector max = " << Size+Size*Size << " allocated = " << u_count_alloc << "\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, index (" + to_string(index)
+ + ") out of range for u vector max = "
+ + to_string(Size+Size*Size)
+ + " allocated = " + to_string(u_count_alloc) + "\n");
if (NZE >= max_nze)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, exceeds the capacity of A_m sparse matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, exceeds the capacity of A_m sparse matrix\n");
#endif
A[NZE] = u[index];
Ai[NZE] = eq;
NZE++;
#ifdef DEBUG
if (eq < 0 || eq >= Size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, index (" << eq << ") out of range for b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, index (" + to_string(eq)
+ + ") out of range for b vector\n");
if (var < 0 || var >= Size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, index (" << var << ") out of range for index_vara vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, index (" + to_string(var)
+ + ") out of range for index_vara vector\n");
if (index_vara[var] < 0 || index_vara[var] >= y_size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, index (" << index_vara[var] << ") out of range for y vector max=" << y_size << " (0)\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, index ("
+ + to_string(index_vara[var])
+ + ") out of range for y vector max=" + to_string(y_size)
+ +" (0)\n");
#endif
- it4++;
}
Aj[Size] = NZE;
}
void
-dynSparseMatrix::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)
+dynSparseMatrix::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
{
- int eq, var;
*b = static_cast<double *>(mxMalloc(Size * sizeof(double)));
test_mxMalloc(*b, __LINE__, __FILE__, __func__, Size * sizeof(double));
if (!(*b))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, can't retrieve b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, can't retrieve b vector\n");
double *x0 = mxGetPr(x0_m);
if (!x0)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse_Simple, can't retrieve x0 vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse_Simple, can't retrieve x0 vector\n");
*Ap = static_cast<SuiteSparse_long *>(mxMalloc((Size+1) * sizeof(SuiteSparse_long)));
test_mxMalloc(*Ap, __LINE__, __FILE__, __func__, (Size+1) * sizeof(SuiteSparse_long));
if (!(*Ap))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, can't allocate Ap index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, can't allocate Ap index vector\n");
size_t prior_nz = IM.size();
*Ai = static_cast<SuiteSparse_long *>(mxMalloc(prior_nz * sizeof(SuiteSparse_long)));
test_mxMalloc(*Ai, __LINE__, __FILE__, __func__, prior_nz * sizeof(SuiteSparse_long));
if (!(*Ai))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, can't allocate Ai index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, can't allocate Ai index vector\n");
*Ax = static_cast<double *>(mxMalloc(prior_nz * sizeof(double)));
test_mxMalloc(*Ax, __LINE__, __FILE__, __func__, prior_nz * sizeof(double));
if (!(*Ax))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, can't retrieve Ax matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
- map<pair<pair<int, int>, int>, int>::iterator it4;
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, can't retrieve Ax matrix\n");
for (int i = 0; i < Size; i++)
{
int eq = index_vara[i];
@@ -702,137 +628,87 @@ dynSparseMatrix::Init_UMFPACK_Sparse_Simple(int Size, map<pair<pair<int, int>, i
(*Ap)[0] = 0;
last_var = 0;
- it4 = IM.begin();
- while (it4 != IM.end())
+ for (auto &[key, index] : IM)
{
- var = it4->first.first.first;
+ auto &[var, ignore, eq] = key;
if (var != last_var)
{
(*Ap)[1+last_var] = NZE;
last_var = var;
}
- eq = it4->first.second;
- int index = it4->second;
#ifdef DEBUG
if (index < 0 || index >= u_count_alloc || index > Size + Size*Size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, index (" << index << ") out of range for u vector max = " << Size+Size*Size << " allocated = " << u_count_alloc << "\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, index (" + to_string(index)
+ + ") out of range for u vector max = "
+ + to_string(Size+Size*Size)
+ + " allocated = " + to_string(u_count_alloc) + "\n");
if (NZE >= max_nze)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, exceeds the capacity of A_m sparse matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, exceeds the capacity of A_m sparse matrix\n");
#endif
(*Ax)[NZE] = u[index];
(*Ai)[NZE] = eq;
NZE++;
#ifdef DEBUG
if (eq < 0 || eq >= Size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, index (" << eq << ") out of range for b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, index (" + to_string(eq)
+ + ") out of range for b vector\n");
if (var < 0 || var >= Size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, index (" << var << ") out of range for index_vara vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, index (" + to_string(var)
+ + ") out of range for index_vara vector\n");
if (index_vara[var] < 0 || index_vara[var] >= y_size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, index (" << index_vara[var] << ") out of range for y vector max=" << y_size << " (0)\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, index ("
+ + to_string(index_vara[var])
+ + ") out of range for y vector max=" + to_string(y_size)
+ + " (0)\n");
#endif
- it4++;
}
(*Ap)[Size] = NZE;
}
int
-dynSparseMatrix::find_exo_num(vector<s_plan> sconstrained_extended_path, int value)
+dynSparseMatrix::find_exo_num(const vector<s_plan> &sconstrained_extended_path, int value)
{
- int res = -1;
- int i = 0;
- for (auto it = sconstrained_extended_path.begin(); it != sconstrained_extended_path.end(); it++, i++)
- if (it->exo_num == value)
- {
- res = i;
- break;
- }
- return res;
+ auto it = find_if(sconstrained_extended_path.begin(), sconstrained_extended_path.end(),
+ [=](auto v) { return v.exo_num == value; });
+ if (it != sconstrained_extended_path.end())
+ return it - sconstrained_extended_path.begin();
+ else
+ return -1;
}
int
-dynSparseMatrix::find_int_date(vector<pair<int, double>> per_value, int value)
+dynSparseMatrix::find_int_date(const vector<pair<int, double>> &per_value, int value)
{
- int res = -1;
- int i = 0;
- for (auto it = per_value.begin(); it != per_value.end(); it++, i++)
- if (it->first == value)
- {
- res = i;
- break;
- }
- return res;
+ auto it = find_if(per_value.begin(), per_value.end(), [=](auto v) { return v.first == value; });
+ if (it != per_value.end())
+ return it - per_value.begin();
+ else
+ return -1;
}
void
-dynSparseMatrix::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)
+dynSparseMatrix::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
{
- int t, eq, var, lag, ti_y_kmin, ti_y_kmax;
- double *jacob_exo;
- int row_x = 0;
-#ifdef DEBUG
- int col_x;
-#endif
int n = periods * Size;
*b = static_cast<double *>(mxMalloc(n * sizeof(double)));
if (!(*b))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, can't retrieve b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, can't retrieve b vector\n");
double *x0 = mxGetPr(x0_m);
if (!x0)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse_Simple, can't retrieve x0 vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse_Simple, can't retrieve x0 vector\n");
*Ap = static_cast<SuiteSparse_long *>(mxMalloc((n+1) * sizeof(SuiteSparse_long)));
test_mxMalloc(*Ap, __LINE__, __FILE__, __func__, (n+1) * sizeof(SuiteSparse_long));
if (!(*Ap))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, can't allocate Ap index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, can't allocate Ap index vector\n");
size_t prior_nz = IM.size() * periods;
*Ai = static_cast<SuiteSparse_long *>(mxMalloc(prior_nz * sizeof(SuiteSparse_long)));
test_mxMalloc(*Ai, __LINE__, __FILE__, __func__, prior_nz * sizeof(SuiteSparse_long));
if (!(*Ai))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, can't allocate Ai index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, can't allocate Ai index vector\n");
*Ax = static_cast<double *>(mxMalloc(prior_nz * sizeof(double)));
test_mxMalloc(*Ax, __LINE__, __FILE__, __func__, prior_nz * sizeof(double));
if (!(*Ax))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, can't retrieve Ax matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
- map<pair<pair<int, int>, int>, int>::iterator it4, it5;
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, can't retrieve Ax matrix\n");
for (int i = 0; i < y_size*(periods+y_kmin); i++)
ya[i] = y[i];
#ifdef DEBUG
@@ -845,6 +721,11 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
(*b)[i] = 0;
x0[i] = y[index_vara[Size*y_kmin+i]];
}
+ double *jacob_exo;
+ int row_x = 0;
+#ifdef DEBUG
+ int col_x;
+#endif
if (vector_table_conditional_local.size())
{
jacob_exo = mxGetPr(jacobian_exo_block[block_num]);
@@ -854,9 +735,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
#endif
}
else
- {
- jacob_exo = nullptr;
- }
+ jacob_exo = nullptr;
#ifdef DEBUG
int local_index;
#endif
@@ -865,33 +744,16 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
bool fliped_exogenous_derivatives_updated = false;
int flip_exo;
(*Ap)[0] = 0;
- for (t = 0; t < periods; t++)
+ for (int t = 0; t < periods; t++)
{
last_var = -1;
- it4 = IM.begin();
- var = 0;
- while (it4 != IM.end())
+ int var = 0;
+ for (auto &[key, value] : IM)
{
- var = it4->first.first.first;
-// #ifdef DEBUG
-// if (var < 0 || var >= Size)
-// {
-// ostringstream tmp;
-// tmp << " in Init_UMFPACK_Sparse, var (" << var << ") out of range\n";
-// throw FatalExceptionHandling(tmp.str());
-// }
-// #endif
- eq = it4->first.second+Size*t;
-// #ifdef DEBUG
-// if (eq < 0 || eq >= Size)
-// {
-// ostringstream tmp;
-// tmp << " in Init_UMFPACK_Sparse, eq (" << eq << ") out of range\n";
-// throw FatalExceptionHandling(tmp.str());
-// }
-// #endif
- lag = -it4->first.first.second;
- int index = it4->second+ (t-lag) * u_count_init;
+ var = get<0>(key);
+ int eq = get<2>(key)+Size*t;
+ int lag = -get<1>(key);
+ int index = value + (t-lag) * u_count_init;
if (var != last_var)
{
(*Ap)[1+last_var + t * Size] = NZE;
@@ -911,7 +773,8 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
}
if (fliped)
{
- if ((t == 0) && (var < (periods+y_kmax)*Size) && (lag == 0) && (vector_table_conditional_local.size()))
+ if (t == 0 && var < (periods+y_kmax)*Size
+ && lag == 0 && vector_table_conditional_local.size())
{
flip_exo = vector_table_conditional_local[var].var_exo;
#ifdef DEBUG
@@ -930,23 +793,20 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
#ifdef DEBUG
if (local_index < 0 || local_index >= Size * periods)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << local_index << ") out of range for b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(local_index)
+ + ") out of range for b vector\n");
if (k + row_x*flip_exo < 0 || k + row_x*flip_exo >= row_x * col_x)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << var+Size*(y_kmin+t+lag) << ") out of range for jacob_exo vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
- if (t+y_kmin+flip_exo*nb_row_x < 0 || t+y_kmin+flip_exo*nb_row_x >= nb_row_x * this->col_x)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << index_vara[var+Size*(y_kmin+t+lag)] << ") out of range for x vector max=" << nb_row_x * this->col_x << "\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(var+Size*(y_kmin+t+lag))
+ + ") out of range for jacob_exo vector\n");
+ if (t+y_kmin+flip_exo*nb_row_x < 0
+ || t+y_kmin+flip_exo*nb_row_x >= nb_row_x * this->col_x)
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(index_vara[var+Size*(y_kmin+t+lag)])
+ + ") out of range for x vector max="
+ + to_string(nb_row_x * this->col_x)
+ + "\n");
#endif
u[k] -= jacob_exo[k + row_x*flip_exo] * x[t+y_kmin+flip_exo*nb_row_x];
}
@@ -954,37 +814,20 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
}
}
}
- /*if (t==0)
- {
- if (min_lag > lag)
- min_lag = lag;
- if (max_lag < lag)
- max_lag = lag;
- }*/
if (var < (periods+y_kmax)*Size)
{
- ti_y_kmin = -min(t, y_kmin);
- ti_y_kmax = min(periods-(t +1), y_kmax);
+ int ti_y_kmin = -min(t, y_kmin);
+ int ti_y_kmax = min(periods-(t+1), y_kmax);
int ti_new_y_kmax = min(t, y_kmax);
int ti_new_y_kmin = -min(periods-(t+1), y_kmin);
if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/
{
#ifdef DEBUG
- // if (index < 0 || index >= u_count_alloc || index > Size + Size*Size)
- // {
- // ostringstream tmp;
- // tmp << " in Init_UMFPACK_Sparse, index (" << index << ") out of range for u vector max = " << Size+Size*Size << " allocated = " << u_count_alloc << "\n";
- // throw FatalExceptionHandling(tmp.str());
- // }
if (NZE >= max_nze)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, exceeds the capacity of A_m sparse matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, exceeds the capacity of A_m sparse matrix\n");
#endif
- if ((!fliped /*|| lag != 0*/) /*&& (!(vector_table_conditional_local[eq-lag*Size].is_cond && (t-lag == 0)))*/)
+ if (!fliped)
{
(*Ax)[NZE] = u[index];
(*Ai)[NZE] = eq - lag * Size;
@@ -994,25 +837,23 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
{
#ifdef DEBUG
if (eq - lag * Size < 0 || eq - lag * Size >= Size * periods)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << eq - lag * Size << ") out of range for b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
- if (var+Size*(y_kmin+t) < 0 || var+Size*(y_kmin+t) >= Size*(periods+y_kmin+y_kmax))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << var+Size*(y_kmin+t) << ") out of range for index_vara vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
- if (index_vara[var+Size*(y_kmin+t /*+lag*/)] < 0 || index_vara[var+Size*(y_kmin+t /*+lag*/)] >= y_size*(periods+y_kmin+y_kmax))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << index_vara[var+Size*(y_kmin+t /*+lag*/)] << ") out of range for y vector max=" << y_size*(periods+y_kmin+y_kmax) << "\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(eq - lag * Size)
+ + ") out of range for b vector\n");
+ if (var+Size*(y_kmin+t) < 0
+ || var+Size*(y_kmin+t) >= Size*(periods+y_kmin+y_kmax))
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(var+Size*(y_kmin+t))
+ + ") out of range for index_vara vector\n");
+ if (index_vara[var+Size*(y_kmin+t)] < 0
+ || index_vara[var+Size*(y_kmin+t)] >= y_size*(periods+y_kmin+y_kmax))
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(index_vara[var+Size*(y_kmin+t)])
+ + ") out of range for y vector max="
+ + to_string(y_size*(periods+y_kmin+y_kmax))
+ + "\n");
#endif
- (*b)[eq - lag * Size] += u[index] * y[index_vara[var+Size*(y_kmin+t /*+lag*/)]];
+ (*b)[eq - lag * Size] += u[index] * y[index_vara[var+Size*(y_kmin+t)]];
}
}
@@ -1020,23 +861,20 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
{
#ifdef DEBUG
if (eq < 0 || eq >= Size * periods)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << eq << ") out of range for b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
- if (var+Size*(y_kmin+t+lag) < 0 || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << var+Size*(y_kmin+t+lag) << ") out of range for index_vara vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
- if (index_vara[var+Size*(y_kmin+t+lag)] < 0 || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << index_vara[var+Size*(y_kmin+t+lag)] << ") out of range for y vector max=" << y_size*(periods+y_kmin+y_kmax) << "\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(eq)
+ + ") out of range for b vector\n");
+ if (var+Size*(y_kmin+t+lag) < 0
+ || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(var+Size*(y_kmin+t+lag))
+ + ") out of range for index_vara vector\n");
+ if (index_vara[var+Size*(y_kmin+t+lag)] < 0
+ || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index ("
+ + to_string(index_vara[var+Size*(y_kmin+t+lag)])
+ + ") out of range for y vector max="
+ + to_string(y_size*(periods+y_kmin+y_kmax)) + "\n");
#endif
(*b)[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]];
}
@@ -1045,21 +883,14 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
{
#ifdef DEBUG
if (index < 0 || index >= u_count_alloc)
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << index << ") out of range for u vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index (" + to_string(index)
+ + ") out of range for u vector\n");
if (eq < 0 || eq >= (Size*periods))
- {
- ostringstream tmp;
- tmp << " in Init_UMFPACK_Sparse, index (" << eq << ") out of range for b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_UMFPACK_Sparse, index (" + to_string(eq)
+ + ") out of range for b vector\n");
#endif
(*b)[eq] += u[index];
}
- it4++;
}
}
(*Ap)[Size*periods] = NZE;
@@ -1082,26 +913,23 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
}
void
-dynSparseMatrix::PrintM(int n, double *Ax, mwIndex *Ap, mwIndex *Ai)
+dynSparseMatrix::PrintM(int n, const double *Ax, const mwIndex *Ap, const mwIndex *Ai)
{
int nnz = Ap[n];
auto *A = static_cast<double *>(mxMalloc(n * n * sizeof(double)));
test_mxMalloc(A, __LINE__, __FILE__, __func__, n * n * sizeof(double));
- memset(A, 0, n * n * sizeof(double));
+ fill_n(A, n * n, 0);
int k = 0;
for (int i = 0; i < n; i++)
- {
- for (int j = Ap[i]; j < static_cast<int>(Ap[i + 1]); j++)
- {
- int row = Ai[j];
- A[row *n + i] = Ax[j];
- k++;
- }
- }
+ for (int j = Ap[i]; j < static_cast<int>(Ap[i + 1]); j++)
+ {
+ int row = Ai[j];
+ A[row * n + i] = Ax[j];
+ k++;
+ }
if (nnz != k)
mexPrintf("Problem nnz(%d) != number of elements(%d)\n", nnz, k);
mexPrintf("----------------------\n");
- //mexEvalString("drawnow;");
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
@@ -1112,47 +940,25 @@ dynSparseMatrix::PrintM(int n, double *Ax, mwIndex *Ap, mwIndex *Ai)
}
void
-dynSparseMatrix::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)
+dynSparseMatrix::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
{
- int t, eq, var, lag, ti_y_kmin, ti_y_kmax;
double *b = mxGetPr(b_m);
if (!b)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, can't retrieve b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, can't retrieve b vector\n");
double *x0 = mxGetPr(x0_m);
if (!x0)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse_Simple, can't retrieve x0 vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse_Simple, can't retrieve x0 vector\n");
mwIndex *Aj = mxGetJc(A_m);
if (!Aj)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, can't allocate Aj index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, can't allocate Aj index vector\n");
mwIndex *Ai = mxGetIr(A_m);
if (!Ai)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, can't allocate Ai index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, can't allocate Ai index vector\n");
double *A = mxGetPr(A_m);
if (!A)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, can't retrieve A matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, can't retrieve A matrix\n");
- map<pair<pair<int, int>, int>, int>::iterator it4;
for (int i = 0; i < y_size*(periods+y_kmin); i++)
ya[i] = y[i];
#ifdef DEBUG
@@ -1166,42 +972,36 @@ dynSparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Siz
x0[i] = y[index_vara[Size*y_kmin+i]];
}
Aj[0] = 0;
- for (t = 0; t < periods; t++)
+ for (int t = 0; t < periods; t++)
{
last_var = 0;
- it4 = IM.begin();
- while (it4 != IM.end())
+ for (auto &[key, value] : IM)
{
- var = it4->first.first.first;
+ int var = get<0>(key);
if (var != last_var)
{
Aj[1+last_var + t * Size] = NZE;
last_var = var;
}
- eq = it4->first.second+Size*t;
- lag = -it4->first.first.second;
- int index = it4->second+ (t-lag) * u_count_init;
+ int eq = get<2>(key)+Size*t;
+ int lag = -get<1>(key);
+ int index = value + (t-lag)*u_count_init;
if (var < (periods+y_kmax)*Size)
{
- ti_y_kmin = -min(t, y_kmin);
- ti_y_kmax = min(periods-(t +1), y_kmax);
+ int ti_y_kmin = -min(t, y_kmin);
+ int ti_y_kmax = min(periods-(t +1), y_kmax);
int ti_new_y_kmax = min(t, y_kmax);
int ti_new_y_kmin = -min(periods-(t+1), y_kmin);
if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/
{
#ifdef DEBUG
if (index < 0 || index >= u_count_alloc || index > Size + Size*Size)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, index (" << index << ") out of range for u vector max = " << Size+Size*Size << " allocated = " << u_count_alloc << "\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, index (" + to_string(index)
+ + ") out of range for u vector max = "
+ + to_string(Size+Size*Size) + " allocated = "
+ + to_string(u_count_alloc) + "\n");
if (NZE >= max_nze)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, exceeds the capacity of A_m sparse matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, exceeds the capacity of A_m sparse matrix\n");
#endif
A[NZE] = u[index];
Ai[NZE] = eq - lag * Size;
@@ -1211,23 +1011,19 @@ dynSparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Siz
{
#ifdef DEBUG
if (eq < 0 || eq >= Size * periods)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, index (" << eq << ") out of range for b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
- if (var+Size*(y_kmin+t+lag) < 0 || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, index (" << var+Size*(y_kmin+t+lag) << ") out of range for index_vara vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
- if (index_vara[var+Size*(y_kmin+t+lag)] < 0 || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, index (" << index_vara[var+Size*(y_kmin+t+lag)] << ") out of range for y vector max=" << y_size*(periods+y_kmin+y_kmax) << "\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, index (" + to_string(eq)
+ + ") out of range for b vector\n");
+ if (var+Size*(y_kmin+t+lag) < 0
+ || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, index ("
+ + to_string(var+Size*(y_kmin+t+lag))
+ + ") out of range for index_vara vector\n");
+ if (index_vara[var+Size*(y_kmin+t+lag)] < 0
+ || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, index ("
+ + to_string(index_vara[var+Size*(y_kmin+t+lag)])
+ + ") out of range for y vector max="
+ + to_string(y_size*(periods+y_kmin+y_kmax)) + "\n");
#endif
b[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]];
}
@@ -1236,33 +1032,23 @@ dynSparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Siz
{
#ifdef DEBUG
if (index < 0 || index >= u_count_alloc)
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, index (" << index << ") out of range for u vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, index (" + to_string(index)
+ + ") out of range for u vector\n");
if (eq < 0 || eq >= (Size*periods))
- {
- ostringstream tmp;
- tmp << " in Init_Matlab_Sparse, index (" << eq << ") out of range for b vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Init_Matlab_Sparse, index (" + to_string(eq)
+ + ") out of range for b vector\n");
#endif
b[eq] += u[index];
}
- it4++;
}
}
Aj[Size*periods] = NZE;
}
void
-dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM)
+dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, const map<tuple<int, int, int>, int> &IM)
{
- int t, i, eq, var, lag, ti_y_kmin, ti_y_kmax;
double tmp_b = 0.0;
- map<pair<pair<int, int>, int>, int>::iterator it4;
- NonZeroElem *first;
pivot = static_cast<int *>(mxMalloc(Size*periods*sizeof(int)));
test_mxMalloc(pivot, __LINE__, __FILE__, __func__, Size*periods*sizeof(int));
pivot_save = static_cast<int *>(mxMalloc(Size*periods*sizeof(int)));
@@ -1280,7 +1066,7 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair
mem_mngr.init_CHUNK_BLCK_SIZE(u_count);
g_save_op = nullptr;
g_nop_all = 0;
- i = (periods+y_kmax+1)*Size*sizeof(NonZeroElem *);
+ int i = (periods+y_kmax+1)*Size*sizeof(NonZeroElem *);
FNZE_R = static_cast<NonZeroElem **>(mxMalloc(i));
test_mxMalloc(FNZE_R, __LINE__, __FILE__, __func__, i);
FNZE_C = static_cast<NonZeroElem **>(mxMalloc(i));
@@ -1311,33 +1097,32 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair
}
int nnz = 0;
//pragma omp parallel for ordered private(it4, ti_y_kmin, ti_y_kmax, eq, var, lag) schedule(dynamic)
- for (t = 0; t < periods; t++)
+ for (int t = 0; t < periods; t++)
{
- ti_y_kmin = -min(t, y_kmin);
- ti_y_kmax = min(periods-(t+1), y_kmax);
- it4 = IM.begin();
- eq = -1;
+ int ti_y_kmin = -min(t, y_kmin);
+ int ti_y_kmax = min(periods-(t+1), y_kmax);
+ int eq = -1;
//pragma omp ordered
- while (it4 != IM.end())
+ for (auto &[key, value] : IM)
{
- var = it4->first.first.second;
- if (eq != it4->first.first.first+Size*t)
+ int var = get<1>(key);
+ if (eq != get<0>(key)+Size*t)
tmp_b = 0;
- eq = it4->first.first.first+Size*t;
- lag = it4->first.second;
+ eq = get<0>(key)+Size*t;
+ int lag = get<2>(key);
if (var < (periods+y_kmax)*Size)
{
- lag = it4->first.second;
+ lag = get<2>(key);
if (lag <= ti_y_kmax && lag >= ti_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/
{
nnz++;
var += Size*t;
NbNZRow[eq]++;
NbNZCol[var]++;
- first = mem_mngr.mxMalloc_NZE();
+ NonZeroElem *first = mem_mngr.mxMalloc_NZE();
first->NZE_C_N = nullptr;
first->NZE_R_N = nullptr;
- first->u_index = it4->second+u_count_init*t;
+ first->u_index = value+u_count_init*t;
first->r_index = eq;
first->c_index = var;
first->lag_index = lag;
@@ -1355,23 +1140,17 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair
else /*Build the additive terms ooutside the simulation periods related to the first lags and the last leads...*/
{
if (lag < ti_y_kmin)
- {
- tmp_b += u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
- }
+ tmp_b += u[value+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
else
- {
- tmp_b += u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
-
- }
+ tmp_b += u[value+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
}
}
else /* ...and store it in the u vector*/
{
- b[eq] = it4->second+u_count_init*t;
+ b[eq] = value+u_count_init*t;
u[b[eq]] += tmp_b;
tmp_b = 0;
}
- it4++;
}
}
mxFree(temp_NZE_R);
@@ -1400,11 +1179,8 @@ dynSparseMatrix::Get_u()
u_count_alloc += 5*u_count_alloc_save;
u = static_cast<double *>(mxRealloc(u, u_count_alloc*sizeof(double)));
if (!u)
- {
- ostringstream tmp;
- tmp << " in Get_u, memory exhausted (realloc(" << u_count_alloc*sizeof(double) << "))\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Get_u, memory exhausted (realloc("
+ + to_string(u_count_alloc*sizeof(double)) + "))\n");
int i = u_count;
u_count++;
return i;
@@ -1425,7 +1201,7 @@ dynSparseMatrix::Clear_u()
}
void
-dynSparseMatrix::Print_u()
+dynSparseMatrix::Print_u() const
{
for (int i : u_liste)
mexPrintf("%d ", i);
@@ -1449,29 +1225,24 @@ dynSparseMatrix::End_GE(int Size)
}
bool
-dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size)
+dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long nop4, int Size)
{
- long int i, j, nop = nop4/2;
+ long nop = nop4/2;
double r = 0.0;
bool OK = true;
- t_save_op_s *save_op_s, *save_opa_s, *save_opaa_s;
- int *diff1, *diff2;
- diff1 = static_cast<int *>(mxMalloc(nop*sizeof(int)));
+ int *diff1 = static_cast<int *>(mxMalloc(nop*sizeof(int)));
test_mxMalloc(diff1, __LINE__, __FILE__, __func__, nop*sizeof(int));
- diff2 = static_cast<int *>(mxMalloc(nop*sizeof(int)));
+ int *diff2 = static_cast<int *>(mxMalloc(nop*sizeof(int)));
test_mxMalloc(diff2, __LINE__, __FILE__, __func__, nop*sizeof(int));
int max_save_ops_first = -1;
- j = i = 0;
+ long j = 0, i = 0;
while (i < nop4 && OK)
{
- save_op_s = reinterpret_cast<t_save_op_s *>(&(save_op[i]));
- save_opa_s = reinterpret_cast<t_save_op_s *>(&(save_opa[i]));
- save_opaa_s = reinterpret_cast<t_save_op_s *>(&(save_opaa[i]));
+ t_save_op_s *save_op_s = reinterpret_cast<t_save_op_s *>(&(save_op[i]));
+ t_save_op_s *save_opa_s = reinterpret_cast<t_save_op_s *>(&(save_opa[i]));
+ t_save_op_s *save_opaa_s = reinterpret_cast<t_save_op_s *>(&(save_opaa[i]));
diff1[j] = save_op_s->first-save_opa_s->first;
- if (max_save_ops_first < save_op_s->first+diff1[j]*(periods-beg_t))
- {
- max_save_ops_first = save_op_s->first+diff1[j]*(periods-beg_t);
- }
+ max_save_ops_first = max(max_save_ops_first, save_op_s->first+diff1[j]*(periods-beg_t));
switch (save_op_s->operat)
{
case IFLD:
@@ -1489,9 +1260,8 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t,
i += 3;
break;
default:
- ostringstream tmp;
- tmp << " in compare, unknown operator = " << save_op_s->operat << "\n";
- throw FatalExceptionHandling(tmp.str());
+ throw FatalExceptionHandling(" in compare, unknown operator = "
+ + to_string(save_op_s->operat) + "\n");
}
j++;
}
@@ -1499,29 +1269,23 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t,
if (OK)
{
for (int i = beg_t; i < periods; i++)
- {
- for (int j = 0; j < Size; j++)
- pivot[i*Size+j] = pivot[(i-1)*Size+j]+Size;
- }
+ for (int j = 0; j < Size; j++)
+ pivot[i*Size+j] = pivot[(i-1)*Size+j]+Size;
if (max_save_ops_first >= u_count_alloc)
{
u_count_alloc += max_save_ops_first;
u = static_cast<double *>(mxRealloc(u, u_count_alloc*sizeof(double)));
if (!u)
- {
- ostringstream tmp;
- tmp << " in compare, memory exhausted (realloc(" << u_count_alloc*sizeof(double) << "))\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in compare, memory exhausted (realloc("
+ + to_string(u_count_alloc*sizeof(double)) + "))\n");
}
for (int t = 1; t < periods-beg_t-y_kmax; t++)
{
int i = j = 0;
- double *up;
while (i < nop4)
{
auto *save_op_s = reinterpret_cast<t_save_op_s *>(&(save_op[i]));
- up = &u[save_op_s->first+t*diff1[j]];
+ double *up = &u[save_op_s->first+t*diff1[j]];
switch (save_op_s->operat)
{
case IFLD:
@@ -1552,8 +1316,8 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t,
int gap = periods_beg_t-t;
while (i < nop4)
{
- auto *save_op_s = reinterpret_cast<t_save_op_s *>(&(save_op[i]));
- if (save_op_s->lag < gap)
+ if (auto *save_op_s = reinterpret_cast<t_save_op_s *>(&(save_op[i]));
+ save_op_s->lag < gap)
{
double *up = &u[save_op_s->first+t*diff1[j]];
switch (save_op_s->operat)
@@ -1577,19 +1341,17 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t,
}
}
else
- {
- switch (save_op_s->operat)
- {
- case IFLD:
- case IFDIV:
- i += 2;
- break;
- case IFSUB:
- case IFLESS:
- i += 3;
- break;
- }
- }
+ switch (save_op_s->operat)
+ {
+ case IFLD:
+ case IFDIV:
+ i += 2;
+ break;
+ case IFSUB:
+ case IFLESS:
+ i += 3;
+ break;
+ }
j++;
}
}
@@ -1602,26 +1364,23 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t,
int
dynSparseMatrix::complete(int beg_t, int Size, int periods, int *b)
{
- long int i, j, k, nop, nopa, nop1, cal_y, nb_var, pos, max_var, min_var;
- NonZeroElem *first;
- int *save_code;
- int *diff;
- double yy = 0.0, err;
+ double yy = 0.0;
int size_of_save_code = (1+y_kmax)*Size*(Size+1+4)/2*4;
- save_code = static_cast<int *>(mxMalloc(size_of_save_code*sizeof(int)));
+ int *save_code = static_cast<int *>(mxMalloc(size_of_save_code*sizeof(int)));
test_mxMalloc(save_code, __LINE__, __FILE__, __func__, size_of_save_code*sizeof(int));
int size_of_diff = (1+y_kmax)*Size*(Size+1+4);
- diff = static_cast<int *>(mxMalloc(size_of_diff*sizeof(int)));
+ int *diff = static_cast<int *>(mxMalloc(size_of_diff*sizeof(int)));
test_mxMalloc(diff, __LINE__, __FILE__, __func__, size_of_diff*sizeof(int));
- cal_y = y_size*y_kmin;
+ long cal_y = y_size*y_kmin;
- i = (beg_t+1)*Size-1;
- nop = 0;
- for (j = i; j > i-Size; j--)
+ long i = (beg_t+1)*Size-1;
+ long nop = 0;
+ for (long j = i; j > i-Size; j--)
{
- pos = pivot[j];
- nb_var = At_Row(pos, &first);
+ long pos = pivot[j];
+ NonZeroElem *first;
+ long nb_var = At_Row(pos, &first);
first = first->NZE_R_N;
nb_var--;
save_code[nop] = IFLDZ;
@@ -1633,7 +1392,7 @@ dynSparseMatrix::complete(int beg_t, int Size, int periods, int *b)
mexPrintf("out of save_code[%d] (bound=%d)\n", nop+2, size_of_save_code);
#endif
nop += 4;
- for (k = 0; k < nb_var; k++)
+ for (long k = 0; k < nb_var; k++)
{
save_code[nop] = IFMUL;
save_code[nop+1] = index_vara[first->c_index]+cal_y;
@@ -1666,18 +1425,19 @@ dynSparseMatrix::complete(int beg_t, int Size, int periods, int *b)
nop += 4;
}
i = beg_t*Size-1;
- nop1 = nopa = 0;
- for (j = i; j > i-Size; j--)
+ long nop1 = 0, nopa = 0;
+ for (long j = i; j > i-Size; j--)
{
- pos = pivot[j];
- nb_var = At_Row(pos, &first);
+ long pos = pivot[j];
+ NonZeroElem *first;
+ long nb_var = At_Row(pos, &first);
first = first->NZE_R_N;
nb_var--;
diff[nopa] = 0;
diff[nopa+1] = 0;
nopa += 2;
nop1 += 4;
- for (k = 0; k < nb_var; k++)
+ for (long k = 0; k < nb_var; k++)
{
diff[nopa] = save_code[nop1+1]-(index_vara[first->c_index]+cal_y);
diff[nopa+1] = save_code[nop1+2]-(first->u_index);
@@ -1712,8 +1472,8 @@ dynSparseMatrix::complete(int beg_t, int Size, int periods, int *b)
nopa += 2;
nop1 += 4;
}
- max_var = (periods+y_kmin)*y_size;
- min_var = y_kmin*y_size;
+ long max_var = (periods+y_kmin)*y_size;
+ long min_var = y_kmin*y_size;
for (int t = periods+y_kmin-1; t >= beg_t+y_kmin; t--)
{
int j = 0, k;
@@ -1728,16 +1488,14 @@ dynSparseMatrix::complete(int beg_t, int Size, int periods, int *b)
case IFMUL:
k = save_code[i+1]+ti*diff[j];
if (k < max_var && k > min_var)
- {
- yy += y[k]*u[save_code[i+2]+ti*diff[j+1]];
- }
+ yy += y[k]*u[save_code[i+2]+ti*diff[j+1]];
break;
case IFADD:
yy = -(yy+u[save_code[i+1]+ti*diff[j]]);
break;
case IFSTP:
k = save_code[i+1]+ti*diff[j];
- err = yy - y[k];
+ double err = yy - y[k];
y[k] += slowc*(err);
break;
}
@@ -1752,9 +1510,6 @@ dynSparseMatrix::complete(int beg_t, int Size, int periods, int *b)
void
dynSparseMatrix::bksub(int tbreak, int last_period, int Size, double slowc_l)
{
- NonZeroElem *first;
- int i, j, k;
- double yy;
for (int i = 0; i < y_size*(periods+y_kmin); i++)
y[i] = ya[i];
if (symbolic && tbreak)
@@ -1766,16 +1521,17 @@ dynSparseMatrix::bksub(int tbreak, int last_period, int Size, double slowc_l)
int ti = (t-y_kmin)*Size;
int cal = y_kmin*Size;
int cal_y = y_size*y_kmin;
- for (i = ti-1; i >= ti-Size; i--)
+ for (int i = ti-1; i >= ti-Size; i--)
{
- j = i+cal;
+ int j = i+cal;
int pos = pivot[i+Size];
+ NonZeroElem *first;
int nb_var = At_Row(pos, &first);
first = first->NZE_R_N;
nb_var--;
int eq = index_vara[j]+y_size;
- yy = 0;
- for (k = 0; k < nb_var; k++)
+ double yy = 0;
+ for (int k = 0; k < nb_var; k++)
{
yy += y[index_vara[first->c_index]+cal_y]*u[first->u_index];
first = first->NZE_R_N;
@@ -1790,20 +1546,18 @@ dynSparseMatrix::bksub(int tbreak, int last_period, int Size, double slowc_l)
void
dynSparseMatrix::simple_bksub(int it_, int Size, double slowc_l)
{
- int i, k;
- double yy;
- NonZeroElem *first;
for (int i = 0; i < y_size; i++)
y[i+it_*y_size] = ya[i+it_*y_size];
- for (i = Size-1; i >= 0; i--)
+ for (int i = Size-1; i >= 0; i--)
{
int pos = pivot[i];
+ NonZeroElem *first;
int nb_var = At_Row(pos, &first);
first = first->NZE_R_N;
nb_var--;
int eq = index_vara[i];
- yy = 0;
- for (k = 0; k < nb_var; k++)
+ double yy = 0;
+ for (int k = 0; k < nb_var; k++)
{
yy += y[index_vara[first->c_index]+it_*y_size]*u[first->u_index];
first = first->NZE_R_N;
@@ -1817,29 +1571,24 @@ dynSparseMatrix::simple_bksub(int it_, int Size, double slowc_l)
void
dynSparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods)
{
- const double epsilon = 1e-7;
+ constexpr double epsilon = 1e-7;
fstream SaveResult;
ostringstream out;
out << "Result" << iter;
- SaveResult.open(out.str().c_str(), ios::in);
+ SaveResult.open(out.str(), ios::in);
if (!SaveResult.is_open())
- {
- ostringstream tmp;
- tmp << " in CheckIt, Result file cannot be opened\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in CheckIt, Result file cannot be opened\n");
mexPrintf("Reading Result...");
int row, col;
SaveResult >> row;
mexPrintf("row=%d\n", row);
SaveResult >> col;
mexPrintf("col=%d\n", col);
- double G1a;
mexPrintf("Allocated\n");
- NonZeroElem *first;
for (int j = 0; j < col; j++)
{
mexPrintf("j=%d ", j);
+ NonZeroElem *first;
int nb_equ = At_Col(j, &first);
mexPrintf("nb_equ=%d\n", nb_equ);
int line;
@@ -1849,6 +1598,7 @@ dynSparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int perio
line = -9999999;
for (int i = 0; i < row; i++)
{
+ double G1a;
SaveResult >> G1a;
if (line == i)
{
@@ -1869,8 +1619,7 @@ dynSparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int perio
}
SaveResult >> row;
mexPrintf("row(2)=%d\n", row);
- double *B;
- B = static_cast<double *>(mxMalloc(row*sizeof(double)));
+ double *B = static_cast<double *>(mxMalloc(row*sizeof(double)));
test_mxMalloc(B, __LINE__, __FILE__, __func__, row*sizeof(double));
for (int i = 0; i < row; i++)
SaveResult >> B[i];
@@ -1878,19 +1627,16 @@ dynSparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int perio
mexPrintf("done\n");
mexPrintf("Comparing...");
for (int i = 0; i < row; i++)
- {
- if (abs(u[b[i]]+B[i]) > epsilon)
- mexPrintf("Problem at i=%d u[b[i]]=%f B[i]=%f\n", i, u[b[i]], B[i]);
- }
+ if (abs(u[b[i]]+B[i]) > epsilon)
+ mexPrintf("Problem at i=%d u[b[i]]=%f B[i]=%f\n", i, u[b[i]], B[i]);
mxFree(B);
}
void
dynSparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b)
{
- const double epsilon = 1e-10;
+ constexpr double epsilon = 1e-10;
Init_GE(periods, y_kmin, y_kmax, Size, IM_i);
- NonZeroElem *first;
int cal_y = y_kmin*Size;
mexPrintf(" ");
for (int i = 0; i < Size; i++)
@@ -1908,6 +1654,7 @@ dynSparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Siz
double res = 0;
int pos = pivot[i];
mexPrintf("pos[%d]=%d", i, pos);
+ NonZeroElem *first;
int nb_var = At_Row(pos, &first);
mexPrintf(" nb_var=%d\n", nb_var);
for (int j = 0; j < nb_var; j++)
@@ -1924,7 +1671,7 @@ dynSparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Siz
}
mxArray *
-dynSparseMatrix::substract_A_B(mxArray *A_m, mxArray *B_m)
+dynSparseMatrix::substract_A_B(const mxArray *A_m, const mxArray *B_m)
{
size_t n_A = mxGetN(A_m);
size_t m_A = mxGetM(A_m);
@@ -1943,7 +1690,7 @@ dynSparseMatrix::substract_A_B(mxArray *A_m, mxArray *B_m)
}
mxArray *
-dynSparseMatrix::Sparse_substract_A_SB(mxArray *A_m, mxArray *B_m)
+dynSparseMatrix::Sparse_substract_A_SB(const mxArray *A_m, const mxArray *B_m)
{
size_t n_B = mxGetN(B_m);
size_t m_B = mxGetM(B_m);
@@ -1966,7 +1713,7 @@ dynSparseMatrix::Sparse_substract_A_SB(mxArray *A_m, mxArray *B_m)
}
mxArray *
-dynSparseMatrix::Sparse_substract_SA_SB(mxArray *A_m, mxArray *B_m)
+dynSparseMatrix::Sparse_substract_SA_SB(const mxArray *A_m, const mxArray *B_m)
{
size_t n_A = mxGetN(A_m);
size_t m_A = mxGetM(A_m);
@@ -2050,7 +1797,7 @@ dynSparseMatrix::Sparse_substract_SA_SB(mxArray *A_m, mxArray *B_m)
}
mxArray *
-dynSparseMatrix::mult_SAT_B(mxArray *A_m, mxArray *B_m)
+dynSparseMatrix::mult_SAT_B(const mxArray *A_m, const mxArray *B_m)
{
size_t n_A = mxGetN(A_m);
size_t m_A = mxGetM(A_m);
@@ -2062,24 +1809,22 @@ dynSparseMatrix::mult_SAT_B(mxArray *A_m, mxArray *B_m)
mxArray *C_m = mxCreateDoubleMatrix(m_A, n_B, mxREAL);
double *C_d = mxGetPr(C_m);
for (int j = 0; j < static_cast<int>(n_B); j++)
- {
- for (unsigned int i = 0; i < n_A; i++)
- {
- double sum = 0;
- size_t nze_A = A_j[i];
- while (nze_A < static_cast<unsigned int>(A_j[i+1]))
- {
- size_t i_A = A_i[nze_A];
- sum += A_d[nze_A++] * B_d[i_A];
- }
- C_d[j*n_A+i] = sum;
- }
- }
+ for (unsigned int i = 0; i < n_A; i++)
+ {
+ double sum = 0;
+ size_t nze_A = A_j[i];
+ while (nze_A < static_cast<unsigned int>(A_j[i+1]))
+ {
+ size_t i_A = A_i[nze_A];
+ sum += A_d[nze_A++] * B_d[i_A];
+ }
+ C_d[j*n_A+i] = sum;
+ }
return C_m;
}
mxArray *
-dynSparseMatrix::Sparse_mult_SAT_B(mxArray *A_m, mxArray *B_m)
+dynSparseMatrix::Sparse_mult_SAT_B(const mxArray *A_m, const mxArray *B_m)
{
size_t n_A = mxGetN(A_m);
size_t m_A = mxGetM(A_m);
@@ -2099,26 +1844,24 @@ dynSparseMatrix::Sparse_mult_SAT_B(mxArray *A_m, mxArray *B_m)
C_j[C_col] = 0;
//#pragma omp parallel for
for (unsigned int j = 0; j < n_B; j++)
- {
- for (unsigned int i = 0; i < n_A; i++)
- {
- double sum = 0;
- size_t nze_A = A_j[i];
- while (nze_A < static_cast<unsigned int>(A_j[i+1]))
- {
- size_t i_A = A_i[nze_A];
- sum += A_d[nze_A++] * B_d[i_A];
- }
- if (fabs(sum) > 1e-10)
- {
- C_d[nze_C] = sum;
- C_i[nze_C] = i;
- while (C_col < j)
- C_j[++C_col] = nze_C;
- nze_C++;
- }
- }
- }
+ for (unsigned int i = 0; i < n_A; i++)
+ {
+ double sum = 0;
+ size_t nze_A = A_j[i];
+ while (nze_A < static_cast<unsigned int>(A_j[i+1]))
+ {
+ size_t i_A = A_i[nze_A];
+ sum += A_d[nze_A++] * B_d[i_A];
+ }
+ if (fabs(sum) > 1e-10)
+ {
+ C_d[nze_C] = sum;
+ C_i[nze_C] = i;
+ while (C_col < j)
+ C_j[++C_col] = nze_C;
+ nze_C++;
+ }
+ }
while (C_col < m_B)
C_j[++C_col] = nze_C;
mxSetNzmax(C_m, nze_C);
@@ -2126,7 +1869,7 @@ dynSparseMatrix::Sparse_mult_SAT_B(mxArray *A_m, mxArray *B_m)
}
mxArray *
-dynSparseMatrix::Sparse_mult_SAT_SB(mxArray *A_m, mxArray *B_m)
+dynSparseMatrix::Sparse_mult_SAT_SB(const mxArray *A_m, const mxArray *B_m)
{
size_t n_A = mxGetN(A_m);
size_t m_A = mxGetM(A_m);
@@ -2145,33 +1888,31 @@ dynSparseMatrix::Sparse_mult_SAT_SB(mxArray *A_m, mxArray *B_m)
unsigned int C_col = 0;
C_j[C_col] = 0;
for (unsigned int j = 0; j < n_B; j++)
- {
- for (unsigned int i = 0; i < n_A; i++)
- {
- double sum = 0;
- nze_B = B_j[j];
- nze_A = A_j[i];
- while (nze_A < static_cast<unsigned int>(A_j[i+1]) && nze_B < static_cast<unsigned int>(B_j[j+1]))
- {
- size_t i_A = A_i[nze_A];
- size_t i_B = B_i[nze_B];
- if (i_A == i_B)
- sum += A_d[nze_A++] * B_d[nze_B++];
- else if (i_A < i_B)
- nze_A++;
- else
- nze_B++;
- }
- if (fabs(sum) > 1e-10)
- {
- C_d[nze_C] = sum;
- C_i[nze_C] = i;
- while (C_col < j)
- C_j[++C_col] = nze_C;
- nze_C++;
- }
- }
- }
+ for (unsigned int i = 0; i < n_A; i++)
+ {
+ double sum = 0;
+ nze_B = B_j[j];
+ nze_A = A_j[i];
+ while (nze_A < static_cast<unsigned int>(A_j[i+1]) && nze_B < static_cast<unsigned int>(B_j[j+1]))
+ {
+ size_t i_A = A_i[nze_A];
+ size_t i_B = B_i[nze_B];
+ if (i_A == i_B)
+ sum += A_d[nze_A++] * B_d[nze_B++];
+ else if (i_A < i_B)
+ nze_A++;
+ else
+ nze_B++;
+ }
+ if (fabs(sum) > 1e-10)
+ {
+ C_d[nze_C] = sum;
+ C_i[nze_C] = i;
+ while (C_col < j)
+ C_j[++C_col] = nze_C;
+ nze_C++;
+ }
+ }
while (C_col < n_B)
C_j[++C_col] = nze_C;
mxSetNzmax(C_m, nze_C);
@@ -2179,7 +1920,7 @@ dynSparseMatrix::Sparse_mult_SAT_SB(mxArray *A_m, mxArray *B_m)
}
mxArray *
-dynSparseMatrix::Sparse_transpose(mxArray *A_m)
+dynSparseMatrix::Sparse_transpose(const mxArray *A_m)
{
size_t n_A = mxGetN(A_m);
size_t m_A = mxGetM(A_m);
@@ -2195,60 +1936,52 @@ dynSparseMatrix::Sparse_transpose(mxArray *A_m)
memset(C_j, 0, m_A);
map<pair<mwIndex, unsigned int>, double> B2;
for (unsigned int i = 0; i < n_A; i++)
- {
- while (nze_A < static_cast<unsigned int>(A_j[i+1]))
- {
- C_j[A_i[nze_A]+1]++;
- B2[make_pair(A_i[nze_A], i)] = A_d[nze_A];
- nze_A++;
- }
- }
+ while (nze_A < static_cast<unsigned int>(A_j[i+1]))
+ {
+ C_j[A_i[nze_A]+1]++;
+ B2[{ A_i[nze_A], i }] = A_d[nze_A];
+ nze_A++;
+ }
for (unsigned int i = 0; i < m_A; i++)
C_j[i+1] += C_j[i];
- for (map<pair<mwIndex, unsigned int>, double>::const_iterator it = B2.begin(); it != B2.end(); it++)
+ for (auto &[key, val] : B2)
{
- C_d[nze_C] = it->second;
- C_i[nze_C++] = it->first.second;
+ C_d[nze_C] = val;
+ C_i[nze_C++] = key.second;
}
return C_m;
}
-#define sign(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
bool
dynSparseMatrix::mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc)
{
- const double GOLD = 1.618034;
- const double GLIMIT = 100.0;
- const double TINY = 1.0e-20;
+ constexpr double GOLD = 1.618034;
+ constexpr double GLIMIT = 100.0;
+ constexpr double TINY = 1.0e-20;
+
+ auto sign = [](double a, double b) { return b >= 0.0 ? fabs(a) : -fabs(a); };
- double tmp;
mexPrintf("bracketing *ax=%f, *bx=%f\n", *ax, *bx);
- //mexEvalString("drawnow;");
- double ulim, u, r, q, fu;
if (!compute_complete(*ax, fa))
return false;
if (!compute_complete(*bx, fb))
return false;
if (*fb > *fa)
{
- tmp = *ax;
- *ax = *bx;
- *bx = tmp;
-
- tmp = *fa;
- *fa = *fb;
- *fb = tmp;
+ swap(*ax, *bx);
+ swap(*fa, *fb);
}
*cx = (*bx)+GOLD*(*bx-*ax);
if (!compute_complete(*cx, fc))
return false;
while (*fb > *fc)
{
- r = (*bx-*ax)*(*fb-*fc);
- q = (*bx-*cx)*(*fb-*fa);
- u = (*bx)-((*bx-*cx)*q-(*bx-*ax)*r)
+ double r = (*bx-*ax)*(*fb-*fc);
+ double q = (*bx-*cx)*(*fb-*fa);
+ double u = (*bx)-((*bx-*cx)*q-(*bx-*ax)*r)
/(2.0*sign(fmax(fabs(q-r), TINY), q-r));
- ulim = (*bx)+GLIMIT*(*cx-*bx);
+ double ulim = (*bx)+GLIMIT*(*cx-*bx);
+ double fu;
if ((*bx-u)*(u-*cx) > 0.0)
{
if (!compute_complete(u, &fu))
@@ -2314,11 +2047,10 @@ dynSparseMatrix::golden(double ax, double bx, double cx, double tol, double solv
const double R = 0.61803399;
const double C = (1.0-R);
mexPrintf("golden\n");
- //mexEvalString("drawnow;");
- double f1, f2, x0, x1, x2, x3;
int iter = 0, max_iter = 100;
- x0 = ax;
- x3 = cx;
+ double f1, f2, x1, x2;
+ double x0 = ax;
+ double x3 = cx;
if (fabs(cx-bx) > fabs(bx-ax))
{
x1 = bx;
@@ -2333,7 +2065,8 @@ dynSparseMatrix::golden(double ax, double bx, double cx, double tol, double solv
return false;
if (!compute_complete(x2, &f2))
return false;
- while ((fabs(x3-x0) > tol*(fabs(x1)+fabs(x2)) && (f1 > solve_tolf && f2 > solve_tolf)) && (iter < max_iter) && (abs(x1 - x2) > 1e-4))
+ while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2)) && f1 > solve_tolf && f2 > solve_tolf
+ && iter < max_iter && abs(x1 - x2) > 1e-4)
{
if (f2 < f1)
{
@@ -2370,81 +2103,59 @@ dynSparseMatrix::golden(double ax, double bx, double cx, double tol, double solv
void
dynSparseMatrix::Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_)
{
- mxArray *B1, *C1, *A2, *B2, *A3, *b1, *b2;
double *b_m_d = mxGetPr(b_m);
if (!b_m_d)
- {
- ostringstream tmp;
- tmp << " in Solve_Matlab_Relaxation, can't retrieve b_m vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Solve_Matlab_Relaxation, can't retrieve b_m vector\n");
mwIndex *A_m_i = mxGetIr(A_m);
if (!A_m_i)
- {
- ostringstream tmp;
- tmp << " in Solve_Matlab_Relaxation, can't allocate A_m_i index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Solve_Matlab_Relaxation, can't allocate A_m_i index vector\n");
mwIndex *A_m_j = mxGetJc(A_m);
if (!A_m_j)
- {
- ostringstream tmp;
- tmp << " in Solve_Matlab_Relaxation, can't allocate A_m_j index vector\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Solve_Matlab_Relaxation, can't allocate A_m_j index vector\n");
double *A_m_d = mxGetPr(A_m);
if (!A_m_d)
- {
- ostringstream tmp;
- tmp << " in Solve_Matlab_Relaxation, can't retrieve A matrix\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" in Solve_Matlab_Relaxation, can't retrieve A matrix\n");
size_t max_nze = A_m_j[Size*periods];
- unsigned int nze = 0;
+ unsigned nze = 0;
size_t var = A_m_j[nze];
- B1 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
+ mxArray *B1 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
mwIndex *B1_i = mxGetIr(B1);
mwIndex *B1_j = mxGetJc(B1);
double *B1_d = mxGetPr(B1);
- unsigned int B1_nze = 0;
- unsigned int B1_var = 0;
+ unsigned B1_nze = 0, B1_var = 0;
B1_i[B1_nze] = 0;
B1_j[B1_var] = 0;
- C1 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
+ mxArray *C1 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
mwIndex *C1_i = mxGetIr(C1);
mwIndex *C1_j = mxGetJc(C1);
double *C1_d = mxGetPr(C1);
- unsigned int C1_nze = 0;
- unsigned int C1_var = 0;
+ unsigned C1_nze = 0, C1_var = 0;
C1_i[C1_nze] = 0;
C1_j[C1_var] = 0;
- A2 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
+ mxArray *A2 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
mwIndex *A2_i = mxGetIr(A2);
mwIndex *A2_j = mxGetJc(A2);
double *A2_d = mxGetPr(A2);
- unsigned int A2_nze = 0;
- unsigned int A2_var = 0;
+ unsigned A2_nze = 0, A2_var = 0;
A2_i[A2_nze] = 0;
A2_j[A2_var] = 0;
- B2 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
+ mxArray *B2 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
mwIndex *B2_i = mxGetIr(B2);
mwIndex *B2_j = mxGetJc(B2);
double *B2_d = mxGetPr(B2);
- unsigned int B2_nze = 0;
- unsigned int B2_var = 0;
+ unsigned B2_nze = 0, B2_var = 0;
B2_i[B2_nze] = 0;
B2_j[B2_var] = 0;
- A3 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
+ mxArray *A3 = mxCreateSparse(Size, Size, Size*Size, mxREAL);
mwIndex *A3_i = mxGetIr(A3);
mwIndex *A3_j = mxGetJc(A3);
double *A3_d = mxGetPr(A3);
- unsigned int A3_nze = 0;
- unsigned int A3_var = 0;
+ unsigned A3_nze = 0, A3_var = 0;
A3_i[A3_nze] = 0;
A3_j[A3_var] = 0;
- b1 = mxCreateDoubleMatrix(Size, 1, mxREAL);
+ mxArray *b1 = mxCreateDoubleMatrix(Size, 1, mxREAL);
double *b1_d = mxGetPr(b1);
- b2 = mxCreateDoubleMatrix(Size, 1, mxREAL);
+ mxArray *b2 = mxCreateDoubleMatrix(Size, 1, mxREAL);
double *b2_d = mxGetPr(b2);
size_t eq = 0;
/*B1 C1
@@ -2646,9 +2357,7 @@ dynSparseMatrix::Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, d
{
size_t n = mxGetM(A_m);
mxArray *z;
- mxArray *rhs[2];
- rhs[0] = A_m;
- rhs[1] = b_m;
+ mxArray *rhs[2] = { A_m, b_m };
mexCallMATLAB(1, &z, 2, rhs, "mldivide");
double *res = mxGetPr(z);
if (is_two_boundaries)
@@ -2684,12 +2393,13 @@ dynSparseMatrix::End_Matlab_LU_UMFPack()
void
dynSparseMatrix::End_Solver()
{
- if (((stack_solve_algo == 0 || stack_solve_algo == 4) && !steady_state) || (solve_algo == 6 && steady_state))
+ if (((stack_solve_algo == 0 || stack_solve_algo == 4) && !steady_state)
+ || (solve_algo == 6 && steady_state))
End_Matlab_LU_UMFPack();
}
void
-dynSparseMatrix::Printfull_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n)
+dynSparseMatrix::Printfull_UMFPack(const SuiteSparse_long *Ap, const SuiteSparse_long *Ai, const double *Ax, const double *b, int n)
{
double A[n*n];
for (int i = 0; i < n*n; i++)
@@ -2707,7 +2417,7 @@ dynSparseMatrix::Printfull_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, d
}
void
-dynSparseMatrix::Print_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, int n)
+dynSparseMatrix::Print_UMFPack(const SuiteSparse_long *Ap, const SuiteSparse_long *Ai, const double *Ax, int n)
{
int k = 0;
for (int i = 0; i < n; i++)
@@ -2716,14 +2426,14 @@ dynSparseMatrix::Print_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, doubl
}
void
-dynSparseMatrix::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)
+dynSparseMatrix::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)
{
- SuiteSparse_long status, sys = 0;
+ SuiteSparse_long sys = 0;
double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO], res[n];
umfpack_dl_defaults(Control);
Control[UMFPACK_PRL] = 5;
- status = 0;
+ SuiteSparse_long status = 0;
if (iter == 0)
{
status = umfpack_dl_symbolic(n, n, Ap, Ai, Ax, &Symbolic, Control, Info);
@@ -2731,9 +2441,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do
{
umfpack_dl_report_info(Control, Info);
umfpack_dl_report_status(Control, status);
- ostringstream Error;
- Error << " umfpack_dl_symbolic failed\n";
- throw FatalExceptionHandling(Error.str());
+ throw FatalExceptionHandling(" umfpack_dl_symbolic failed\n");
}
}
if (iter > 0)
@@ -2743,18 +2451,14 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do
{
umfpack_dl_report_info(Control, Info);
umfpack_dl_report_status(Control, status);
- ostringstream Error;
- Error << " umfpack_dl_numeric failed\n";
- throw FatalExceptionHandling(Error.str());
+ throw FatalExceptionHandling(" umfpack_dl_numeric failed\n");
}
status = umfpack_dl_solve(sys, Ap, Ai, Ax, res, b, Numeric, Control, Info);
if (status != UMFPACK_OK)
{
umfpack_dl_report_info(Control, Info);
umfpack_dl_report_status(Control, status);
- ostringstream Error;
- Error << " umfpack_dl_solve failed\n";
- throw FatalExceptionHandling(Error.str());
+ throw FatalExceptionHandling(" umfpack_dl_solve failed\n");
}
if (vector_table_conditional_local.size())
@@ -2831,12 +2535,12 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do
void
dynSparseMatrix::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_)
{
- SuiteSparse_long status, sys = 0;
+ SuiteSparse_long sys = 0;
double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO], res[n];
umfpack_dl_defaults(Control);
Control[UMFPACK_PRL] = 5;
- status = 0;
+ SuiteSparse_long status = 0;
if (iter == 0)
{
status = umfpack_dl_symbolic(n, n, Ap, Ai, Ax, &Symbolic, Control, Info);
@@ -2844,9 +2548,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do
{
umfpack_dl_report_info(Control, Info);
umfpack_dl_report_status(Control, status);
- ostringstream Error;
- Error << " umfpack_dl_symbolic failed\n";
- throw FatalExceptionHandling(Error.str());
+ throw FatalExceptionHandling(" umfpack_dl_symbolic failed\n");
}
}
if (iter > 0)
@@ -2856,18 +2558,14 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do
{
umfpack_dl_report_info(Control, Info);
umfpack_dl_report_status(Control, status);
- ostringstream Error;
- Error << " umfpack_dl_numeric failed\n";
- throw FatalExceptionHandling(Error.str());
+ throw FatalExceptionHandling(" umfpack_dl_numeric failed\n");
}
status = umfpack_dl_solve(sys, Ap, Ai, Ax, res, b, Numeric, Control, Info);
if (status != UMFPACK_OK)
{
umfpack_dl_report_info(Control, Info);
umfpack_dl_report_status(Control, status);
- ostringstream Error;
- Error << " umfpack_dl_solve failed\n";
- throw FatalExceptionHandling(Error.str());
+ throw FatalExceptionHandling(" umfpack_dl_solve failed\n");
}
if (is_two_boundaries)
@@ -2918,7 +2616,7 @@ dynSparseMatrix::Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double s
status = umfpack_dl_solve(sys, Ap, Ai, Ax, res, B, Numeric, Control, Info);
if (status != UMFPACK_OK)
umfpack_dl_report_info(nullptr, Info);
- //double *res = mxGetPr(z);
+
if (is_two_boundaries)
for (int i = 0; i < n; i++)
{
@@ -2939,78 +2637,6 @@ dynSparseMatrix::Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double s
mxDestroyArray(b_m);
}
-void
-Solve(double *Ax, int *Ap, int *Ai, double *b, int n, bool Lower, double *x)
-{
- if (Lower)
- {
- for (int i = 0; i < n; i++)
- {
- double sum = 0;
- for (int j = Ap[i]; j < Ap[i+1]; j++)
- {
- int k = Ai[j];
- if (k < i)
- sum += x[k] * Ax[j];
- }
- x[i] = b[i] - sum;
- }
- }
- else
- {
- for (int i = n-1; i >= 0; i--)
- {
- double sum = 0, mul = 1;
- for (int j = Ap[i]; j < Ap[i+1]; j++)
- {
- int k = Ai[j];
- if (k > i)
- sum += x[k] * Ax[j];
- else if (k == i)
- mul = Ax[j];
- }
- x[i] = (b[i] - sum) / mul;
- }
- }
-}
-
-void
-Check(int n, double *Ax, int *Ap, int *Ai, double *b, double *x, bool Lower)
-{
- if (Lower)
- {
- for (int i = 0; i < n; i++)
- {
- double sum = 0;
- for (int j = Ap[i]; j < Ap[i+1]; j++)
- {
- int k = Ai[j];
- if (k < i)
- sum += x[k] * Ax[j];
- }
- double err = b[i] - sum - x[i];
- if (abs(err) > 1e-10)
- mexPrintf("error at i=%d\n", i);
- }
- }
- else
- {
- for (int i = n-1; i >= 0; i--)
- {
- double sum = 0;
- for (int j = Ap[i]; j < Ap[i+1]; j++)
- {
- int k = Ai[j];
- if (k >= i)
- sum += x[k] * Ax[j];
- }
- double err = b[i] - sum;
- if (abs(err) > 1e-10)
- mexPrintf("error at i=%d\n", i);
- }
- }
-}
-
void
dynSparseMatrix::Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray *x0_m)
{
@@ -3021,23 +2647,14 @@ dynSparseMatrix::Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double
mxSetFieldByNumber(Setup, 0, 0, mxCreateDoubleScalar(lu_inc_tol));
mxSetFieldByNumber(Setup, 0, 1, mxCreateString("ilutp"));
mxArray *lhs0[2];
- mxArray *rhs0[2];
- rhs0[0] = A_m;
- rhs0[1] = Setup;
+ mxArray *rhs0[2] = { A_m, Setup };
if (mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"))
throw FatalExceptionHandling("In GMRES, the incomplet LU decomposition (ilu) ahs failed.");
mxArray *L1 = lhs0[0];
mxArray *U1 = lhs0[1];
/*[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);*/
- mxArray *rhs[8];
- rhs[0] = A_m;
- rhs[1] = b_m;
- rhs[2] = mxCreateDoubleScalar(Size);
- rhs[3] = mxCreateDoubleScalar(1e-6);
- rhs[4] = mxCreateDoubleScalar(static_cast<double>(n));
- rhs[5] = L1;
- rhs[6] = U1;
- rhs[7] = x0_m;
+ mxArray *rhs[8] = { A_m, b_m, mxCreateDoubleScalar(Size), mxCreateDoubleScalar(1e-6),
+ mxCreateDoubleScalar(static_cast<double>(n)), L1, U1, x0_m };
mxArray *lhs[2];
mexCallMATLAB(2, lhs, 8, rhs, "gmres");
mxArray *z = lhs[0];
@@ -3051,22 +2668,15 @@ dynSparseMatrix::Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double
mxDestroyArray(rhs[6]);
if (*flag1 > 0)
{
- ostringstream tmp;
if (*flag1 == 1)
- {
- tmp << "Error in bytecode: No convergence inside GMRES, in block " << block+1;
- mexWarnMsgTxt(tmp.str().c_str());
- }
+ mexWarnMsgTxt(("Error in bytecode: No convergence inside GMRES, in block "
+ + to_string(block+1)).c_str());
else if (*flag1 == 2)
- {
- tmp << "Error in bytecode: Preconditioner is ill-conditioned, in block " << block+1;
- mexWarnMsgTxt(tmp.str().c_str());
- }
+ mexWarnMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
+ + to_string(block+1)).c_str());
else if (*flag1 == 3)
- {
- tmp << "Error in bytecode: GMRES stagnated (Two consecutive iterates were the same.), in block " << block+1;
- mexWarnMsgTxt(tmp.str().c_str());
- }
+ mexWarnMsgTxt(("Error in bytecode: GMRES stagnated (Two consecutive iterates were the same.), in block "
+ + to_string(block+1)).c_str());
lu_inc_tol /= 10;
}
else
@@ -3101,17 +2711,12 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
/* precond = 0 => Jacobi
precond = 1 => Incomplet LU decomposition*/
size_t n = mxGetM(A_m);
- mxArray *L1, *U1, *Diag;
- L1 = nullptr;
- U1 = nullptr;
- Diag = nullptr;
+ mxArray *L1 = nullptr, *U1 = nullptr, *Diag = nullptr;
- mxArray *rhs0[4];
if (preconditioner == 0)
{
mxArray *lhs0[1];
- rhs0[0] = A_m;
- rhs0[1] = mxCreateDoubleScalar(0);
+ mxArray *rhs0[2] = { A_m, mxCreateDoubleScalar(0) };
mexCallMATLAB(1, lhs0, 2, rhs0, "spdiags");
mxArray *tmp = lhs0[0];
double *tmp_val = mxGetPr(tmp);
@@ -3131,12 +2736,8 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
{
/*[L1, U1] = ilu(g1a=;*/
const char *field_names[] = {"type", "droptol", "milu", "udiag", "thresh"};
- const int type = 0;
- const int droptol = 1;
- const int milu = 2;
- const int udiag = 3;
- const int thresh = 4;
- mwSize dims[1] = {static_cast<mwSize>(1) };
+ const int type = 0, droptol = 1, milu = 2, udiag = 3, thresh = 4;
+ mwSize dims[1] = { static_cast<mwSize>(1) };
mxArray *Setup = mxCreateStructArray(1, dims, 5, field_names);
mxSetFieldByNumber(Setup, 0, type, mxCreateString("ilutp"));
mxSetFieldByNumber(Setup, 0, droptol, mxCreateDoubleScalar(lu_inc_tol));
@@ -3144,22 +2745,15 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
mxSetFieldByNumber(Setup, 0, udiag, mxCreateDoubleScalar(0));
mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(1));
mxArray *lhs0[2];
- mxArray *rhs0[2];
- rhs0[0] = A_m;
- rhs0[1] = Setup;
+ mxArray *rhs0[2] = { A_m, Setup };
if (mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"))
- {
- ostringstream tmp;
- tmp << " In BiCGStab, the incomplet LU decomposition (ilu) ahs failed.\n";
- throw FatalExceptionHandling(tmp.str());
- }
+ throw FatalExceptionHandling(" In BiCGStab, the incomplet LU decomposition (ilu) ahs failed.\n");
L1 = lhs0[0];
U1 = lhs0[1];
mxDestroyArray(Setup);
}
double flags = 2;
- mxArray *z;
- z = nullptr;
+ mxArray *z = nullptr;
if (steady_state) /*Octave BicStab algorihtm involves a 0 division in case of a preconditionner equal to the LU decomposition of A matrix*/
{
mxArray *res = mult_SAT_B(Sparse_transpose(A_m), x0_m);
@@ -3167,14 +2761,11 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
double *b = mxGetPr(b_m);
for (int i = 0; i < static_cast<int>(n); i++)
resid[i] = b[i] - resid[i];
- mxArray *rhs[2];
mxArray *lhs[1];
- rhs[0] = L1;
- rhs[1] = res;
- mexCallMATLAB(1, lhs, 2, rhs, "mldivide");
- rhs[0] = U1;
- rhs[1] = lhs[0];
+ mxArray *rhs[2] = { L1, res };
mexCallMATLAB(1, lhs, 2, rhs, "mldivide");
+ mxArray *rhs2[2] = { U1, lhs[0] };
+ mexCallMATLAB(1, lhs, 2, rhs2, "mldivide");
z = lhs[0];
double *phat = mxGetPr(z);
double *x0 = mxGetPr(x0_m);
@@ -3196,20 +2787,14 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
flags = 0;
mxDestroyArray(res);
}
- //else
if (flags == 2)
{
if (preconditioner == 0)
{
/*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
- mxArray *rhs[5];
- rhs[0] = A_m;
- rhs[1] = b_m;
- rhs[2] = mxCreateDoubleScalar(1e-6);
- rhs[3] = mxCreateDoubleScalar(static_cast<double>(n));
- rhs[4] = Diag;
- //rhs[5] = x0_m;
+ mxArray *rhs[5] = { A_m, b_m, mxCreateDoubleScalar(1e-6),
+ mxCreateDoubleScalar(static_cast<double>(n)), Diag };
mxArray *lhs[2];
mexCallMATLAB(2, lhs, 5, rhs, "bicgstab");
z = lhs[0];
@@ -3224,14 +2809,8 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
else if (preconditioner == 1)
{
/*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
- mxArray *rhs[7];
- rhs[0] = A_m;
- rhs[1] = b_m;
- rhs[2] = mxCreateDoubleScalar(1e-6);
- rhs[3] = mxCreateDoubleScalar(static_cast<double>(n));
- rhs[4] = L1;
- rhs[5] = U1;
- rhs[6] = x0_m;
+ mxArray *rhs[7] = { A_m, b_m, mxCreateDoubleScalar(1e-6),
+ mxCreateDoubleScalar(static_cast<double>(n)), L1, U1, x0_m };
mxArray *lhs[2];
mexCallMATLAB(2, lhs, 7, rhs, "bicgstab");
z = lhs[0];
@@ -3248,22 +2827,15 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
if (flags > 0)
{
- ostringstream tmp;
if (flags == 1)
- {
- tmp << "Error in bytecode: No convergence inside BiCGStab, in block " << block+1;
- mexWarnMsgTxt(tmp.str().c_str());
- }
+ mexWarnMsgTxt(("Error in bytecode: No convergence inside BiCGStab, in block "
+ + to_string(block+1)).c_str());
else if (flags == 2)
- {
- tmp << "Error in bytecode: Preconditioner is ill-conditioned, in block " << block+1;
- mexWarnMsgTxt(tmp.str().c_str());
- }
+ mexWarnMsgTxt(("Error in bytecode: Preconditioner is ill-conditioned, in block "
+ + to_string(block+1)).c_str());
else if (flags == 3)
- {
- tmp << "Error in bytecode: BiCGStab stagnated (Two consecutive iterates were the same.), in block " << block+1;
- mexWarnMsgTxt(tmp.str().c_str());
- }
+ mexWarnMsgTxt(("Error in bytecode: BiCGStab stagnated (Two consecutive iterates were the same.), in block "
+ + to_string(block+1)).c_str());
lu_inc_tol /= 10;
}
else
@@ -3297,8 +2869,7 @@ dynSparseMatrix::Singular_display(int block, int Size)
bool zero_solution;
Simple_Init(Size, IM_i, zero_solution);
NonZeroElem *first;
- mxArray *rhs[1];
- rhs[0] = mxCreateDoubleMatrix(Size, Size, mxREAL);
+ mxArray *rhs[1] = { mxCreateDoubleMatrix(Size, Size, mxREAL) };
double *pind;
pind = mxGetPr(rhs[0]);
for (int j = 0; j < Size * Size; j++)
@@ -3483,7 +3054,7 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i
if (fabs(piv_v[j]) > 0)
{
if (markowitz_c > 0)
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(static_cast<double>(NR[j])/static_cast<double>(N_max)));
else
markovitz = fabs(piv_v[j])/piv_abs;
}
@@ -3510,7 +3081,7 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i
if (fabs(piv_v[j]) > 0)
{
if (markowitz_c > 0)
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(static_cast<double>(NR[j])/static_cast<double>(N_max)));
else
markovitz = fabs(piv_v[j])/piv_abs;
}
@@ -3757,7 +3328,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
if (fabs(piv_v[j]) > 0)
{
if (markowitz_c > 0)
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(static_cast<double>(NR[j])/static_cast<double>(N_max)));
else
markovitz = fabs(piv_v[j])/piv_abs;
}
@@ -3785,7 +3356,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
if (fabs(piv_v[j]) > 0)
{
if (markowitz_c > 0)
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(static_cast<double>(NR[j])/static_cast<double>(N_max)));
else
markovitz = fabs(piv_v[j])/piv_abs;
}
@@ -4149,7 +3720,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
}
if (symbolic)
{
- if (t > int (periods*0.35))
+ if (t > static_cast<int>(periods*0.35))
{
symbolic = false;
mxFree(save_opaa);
@@ -4158,7 +3729,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
}
else if (record && (nop == nop1))
{
- if (t > int (periods*0.35))
+ if (t > static_cast<int>(periods*0.35))
{
symbolic = false;
if (save_opaa)
@@ -4183,7 +3754,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
{
tbreak = t;
tbreak_g = tbreak;
- //mexPrintf("time=%f\n",(1000.0*(double (clock())-double (time11)))/double (CLOCKS_PER_SEC));
+ //mexPrintf("time=%f\n",(1000.0*(static_cast<double>(clock())-static_cast<double>(time11)))/static_cast<double>(CLOCKS_PER_SEC));
break;
}
}
@@ -4220,14 +3791,14 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
nop2 = nop1;
nop1 = nop;
}
- //mexPrintf("time=%f\n",(1000.0*(double (clock())-double (time11)))/double (CLOCKS_PER_SEC));
+ //mexPrintf("time=%f\n",(1000.0*(static_cast<double>(clock())-static_cast<double>(time11)))/static_cast<double>(CLOCKS_PER_SEC));
}
mxFree(bc);
mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
mxFree(NR);
- /*mexPrintf("tbreak=%d, periods=%d time required=%f\n",tbreak,periods, (1000.0*(double (clock())-double (time00)))/double (CLOCKS_PER_SEC));
+ /*mexPrintf("tbreak=%d, periods=%d time required=%f\n",tbreak,periods, (1000.0*(static_cast<double>(clock())-static_cast<double>(time00)))/static_cast<double>(CLOCKS_PER_SEC));
mexEvalString("drawnow;");
time00 = clock();*/
nop_all += nop;
@@ -4247,7 +3818,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
ya[i] = y[i];
slowc_save = slowc;
bksub(tbreak, last_period, Size, slowc_lbx);
- /*mexPrintf("remaining operations and bksub time required=%f\n",tbreak,periods, (1000.0*(double (clock())-double (time00)))/double (CLOCKS_PER_SEC));
+ /*mexPrintf("remaining operations and bksub time required=%f\n",tbreak,periods, (1000.0*(static_cast<double>(clock())-static_cast<double>(time00)))/static_cast<double>(CLOCKS_PER_SEC));
mexEvalString("drawnow;");*/
End_GE(Size);
}
@@ -4539,9 +4110,9 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in
mexPrintf("-----------------------------------\n");
mexPrintf(" Simulate iteration no %d \n", iter+1);
- mexPrintf(" max. error=%.10e \n", double (max_res));
- mexPrintf(" sqr. error=%.10e \n", double (res2));
- mexPrintf(" abs. error=%.10e \n", double (res1));
+ mexPrintf(" max. error=%.10e \n", static_cast<double>(max_res));
+ mexPrintf(" sqr. error=%.10e \n", static_cast<double>(res2));
+ mexPrintf(" abs. error=%.10e \n", static_cast<double>(res1));
mexPrintf("-----------------------------------\n");
}
bool zero_solution;
@@ -4557,7 +4128,7 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in
tmp << " in Simulate_One_Boundary, can't allocate b_m vector\n";
throw FatalExceptionHandling(tmp.str());
}
- A_m = mxCreateSparse(size, size, min(int (IM_i.size()*2), size * size), mxREAL);
+ A_m = mxCreateSparse(size, size, min(static_cast<int>(IM_i.size()*2), size * size), mxREAL);
if (!A_m)
{
ostringstream tmp;
@@ -4620,7 +4191,7 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in
}
bool
-dynSparseMatrix::solve_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size, const int iter)
+dynSparseMatrix::solve_linear(int block_num, int y_size, int y_kmin, int y_kmax, int size, int iter)
{
bool cvg = false;
double crit_opt_old = res2/2;
@@ -4638,7 +4209,7 @@ dynSparseMatrix::solve_linear(const int block_num, const int y_size, const int y
}
void
-dynSparseMatrix::solve_non_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size)
+dynSparseMatrix::solve_non_linear(int block_num, int y_size, int y_kmin, int y_kmax, int size)
{
max_res_idx = 0;
bool cvg = false;
@@ -4776,7 +4347,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
{
if (print_it)
{
- mexPrintf("Sim : %f ms\n", (1000.0*(double (clock())-double (time00)))/double (CLOCKS_PER_SEC));
+ mexPrintf("Sim : %f ms\n", (1000.0*(static_cast<double>(clock())-static_cast<double>(time00)))/static_cast<double>(CLOCKS_PER_SEC));
mexEvalString("drawnow;");
}
time00 = clock();
@@ -4925,9 +4496,9 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
}
mexPrintf("-----------------------------------\n");
mexPrintf(" Simulate iteration no %d \n", iter+1);
- mexPrintf(" max. error=%.10e \n", double (max_res));
- mexPrintf(" sqr. error=%.10e \n", double (res2));
- mexPrintf(" abs. error=%.10e \n", double (res1));
+ mexPrintf(" max. error=%.10e \n", static_cast<double>(max_res));
+ mexPrintf(" sqr. error=%.10e \n", static_cast<double>(res2));
+ mexPrintf(" abs. error=%.10e \n", static_cast<double>(res1));
mexPrintf("-----------------------------------\n");
mexEvalString("drawnow;");
}
@@ -4985,7 +4556,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
if (print_it)
{
clock_t t2 = clock();
- mexPrintf("(** %f milliseconds **)\n", 1000.0*(double (t2) - double (t1))/double (CLOCKS_PER_SEC));
+ mexPrintf("(** %f milliseconds **)\n", 1000.0*(static_cast<double>(t2) - static_cast<double>(t1))/static_cast<double>(CLOCKS_PER_SEC));
mexEvalString("drawnow;");
}
if ((!steady_state && (stack_solve_algo == 4 /*|| stack_solve_algo == 0*/)) /* || steady_state*/)
@@ -5000,7 +4571,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
return;
slowc = xmin;
clock_t t3 = clock();
- mexPrintf("(** %f milliseconds **)\n", 1000.0*(double (t3) - double (t2))/double (CLOCKS_PER_SEC));
+ mexPrintf("(** %f milliseconds **)\n", 1000.0*(static_cast<double>(t3) - static_cast<double>(t2))/static_cast<double>(CLOCKS_PER_SEC));
mexEvalString("drawnow;");
}
time00 = clock();
diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh
index e662c0a218..2721652b63 100644
--- a/mex/sources/bytecode/SparseMatrix.hh
+++ b/mex/sources/bytecode/SparseMatrix.hh
@@ -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;
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index b93e357d8d..53764e59c4 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -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;
}
diff --git a/preprocessor b/preprocessor
index 07391e0f9f..702cf62ee3 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 07391e0f9fdc8e697131cd3e292d6c9f23afd63b
+Subproject commit 702cf62ee386b6d240d9d4fcbb0ddeb33bb35c2f
--
GitLab