diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh index 6a019cff7..24c6bac7b 100644 --- a/mex/sources/bytecode/ErrorHandling.hh +++ b/mex/sources/bytecode/ErrorHandling.hh @@ -20,17 +20,22 @@ #ifndef ERROR_HANDLING #define ERROR_HANDLING -#include -#include -#include +#include +#include +#include #include +#include +#include +#include +#include #include -#define BYTE_CODE -#include "CodeInterpreter.hh" - #define _USE_MATH_DEFINES #include -#include + +#include "dynmex.h" + +#define BYTE_CODE +#include "CodeInterpreter.hh" #ifdef OCTAVE_MEX_FILE # define CHAR_LENGTH 1 @@ -187,9 +192,9 @@ public: ExpressionType EQN_type; it_code_type it_code_expr; - /*unsigned int*/ size_t nb_endo, nb_exo, nb_param; + size_t nb_endo, nb_exo, nb_param; char *P_endo_names, *P_exo_names, *P_param_names; - size_t /*unsigned int*/ endo_name_length, exo_name_length, param_name_length; + size_t 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> Variable_list; @@ -243,10 +248,7 @@ public: string temp; int pos1 = -1, pos2 = -1; string tmp_n(str.length(), ' '); - string dollar, pound, tilde; - dollar = "$"; - pound = "£"; - tilde = "~"; + string dollar{"$"}, pound{"£"}, tilde{"~"}; for (const char & i : str) { if (dollar.compare(&i) != 0 && pound.compare(&i) != 0) @@ -359,7 +361,7 @@ public: inline string error_location(bool evaluate, bool steady_state, int size, int block_num, int it_, int Per_u_) { - stringstream Error_loc; + ostringstream Error_loc; if (!steady_state) switch (EQN_type) { @@ -406,7 +408,7 @@ public: Error_loc << "first order derivative of equation " << EQN_equation+1 << " with respect to parameter " << get_variable(SymbolType::endogenous, EQN_dvar1) << " at time " << it_; break; default: - return ("???"); + return "???"; } else switch (EQN_type) diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc index bf3d80391..ce65868cd 100644 --- a/mex/sources/bytecode/Evaluate.cc +++ b/mex/sources/bytecode/Evaluate.cc @@ -17,7 +17,6 @@ * along with Dynare. If not, see . */ -#include #include #include #include diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh index 9526042e4..eb1e8307e 100644 --- a/mex/sources/bytecode/Evaluate.hh +++ b/mex/sources/bytecode/Evaluate.hh @@ -20,13 +20,13 @@ #ifndef EVALUATE_HH_INCLUDED #define EVALUATE_HH_INCLUDED -#include #include #include -#include + +#include "dynmex.h" + #define BYTE_CODE #include "CodeInterpreter.hh" -#include #include "ErrorHandling.hh" class Evaluate : public ErrorMsg diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index 88a8cb0a2..c02fe3a4b 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -17,9 +17,10 @@ * along with Dynare. If not, see . */ -#include #include #include +#include + #include "Interpreter.hh" constexpr double BIG = 1.0e+8, SMALL = 1.0e-5; @@ -440,7 +441,7 @@ Interpreter::simulate_a_block(const vector_table_conditional_local_type &vector_ res1 = 0; max_res = 0; max_res_idx = 0; - memcpy(y_save, y, y_size*sizeof(double)*(periods+y_kmax+y_kmin)); + copy_n(y, y_size*(periods+y_kmax+y_kmin), y_save); if (vector_table_conditional_local.size()) for (auto & it1 : vector_table_conditional_local) if (it1.is_cond) @@ -450,7 +451,7 @@ Interpreter::simulate_a_block(const vector_table_conditional_local_type &vector_ if (!(isnan(res1) || isinf(res1))) cvg = (max_res < solve_tolf); if (isnan(res1) || isinf(res1) || (stack_solve_algo == 4 && iter > 0)) - memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin)); + copy_n(y_save, y_size*(periods+y_kmax+y_kmin), y); u_count = u_count_saved; int prev_iter = iter; Simulate_Newton_Two_Boundaries(block_num, symbol_table_endo_nbr, y_kmin, y_kmax, size, periods, cvg, minimal_solving_periods, stack_solve_algo, endo_name_length, P_endo_names, vector_table_conditional_local); @@ -736,9 +737,9 @@ Interpreter::MainLoop(const string &bin_basename, const CodeLoad &code, bool eva mxFree(T); if (global_temporary_terms) { - if (GlobalTemporaryTerms == nullptr) + if (!GlobalTemporaryTerms) { - mexPrintf("GlobalTemporaryTerms is NULL\n"); + mexPrintf("GlobalTemporaryTerms is nullptr\n"); mexEvalString("drawnow;"); } if (var != static_cast(mxGetNumberOfElements(GlobalTemporaryTerms))) diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh index 08c594475..bd29bd67d 100644 --- a/mex/sources/bytecode/Interpreter.hh +++ b/mex/sources/bytecode/Interpreter.hh @@ -20,15 +20,16 @@ #ifndef INTERPRETER_HH_INCLUDED #define INTERPRETER_HH_INCLUDED -#include #include #include -#include +#include + +#include "dynmex.h" + +#include "ErrorHandling.hh" +#include "SparseMatrix.hh" #define BYTE_CODE #include "CodeInterpreter.hh" -#include "SparseMatrix.hh" -#include "Evaluate.hh" -#include using namespace std; diff --git a/mex/sources/bytecode/Mem_Mngr.cc b/mex/sources/bytecode/Mem_Mngr.cc index 0446e069b..de6c92ded 100644 --- a/mex/sources/bytecode/Mem_Mngr.cc +++ b/mex/sources/bytecode/Mem_Mngr.cc @@ -19,6 +19,8 @@ #include "Mem_Mngr.hh" +#include "dynmex.h" + Mem_Mngr::Mem_Mngr() { swp_f = false; diff --git a/mex/sources/bytecode/Mem_Mngr.hh b/mex/sources/bytecode/Mem_Mngr.hh index d02e96722..58420a9a6 100644 --- a/mex/sources/bytecode/Mem_Mngr.hh +++ b/mex/sources/bytecode/Mem_Mngr.hh @@ -20,10 +20,13 @@ #ifndef MEM_MNGR_HH_INCLUDED #define MEM_MNGR_HH_INCLUDED -#include "ErrorHandling.hh" +#include #include #include -#include + +#include "ErrorHandling.hh" + +using namespace std; struct NonZeroElem { diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 06b083ddc..bb4e2eaf4 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -17,15 +17,11 @@ * along with Dynare. If not, see . */ -#include -#include #include #include #include "SparseMatrix.hh" -using namespace std; - dynSparseMatrix::dynSparseMatrix() { pivotva = nullptr; @@ -174,7 +170,7 @@ dynSparseMatrix::Delete(int r, int c) firsta = first; first = first->NZE_R_N; } - if (firsta != nullptr) + if (firsta) firsta->NZE_R_N = first->NZE_R_N; if (first == FNZE_R[r]) FNZE_R[r] = first->NZE_R_N; @@ -188,7 +184,7 @@ dynSparseMatrix::Delete(int r, int c) first = first->NZE_C_N; } - if (firsta != nullptr) + if (firsta) firsta->NZE_C_N = first->NZE_C_N; if (first == FNZE_C[c]) FNZE_C[c] = first->NZE_C_N; @@ -265,7 +261,7 @@ dynSparseMatrix::Insert(int r, int c, int u_index, int lag_index) { if (first == FNZE_R[r]) FNZE_R[r] = firstn; - if (firsta != nullptr) + if (firsta) firsta->NZE_R_N = firstn; firstn->NZE_R_N = first; } @@ -286,7 +282,7 @@ dynSparseMatrix::Insert(int r, int c, int u_index, int lag_index) { if (first == FNZE_C[c]) FNZE_C[c] = firstn; - if (firsta != nullptr) + if (firsta) firsta->NZE_C_N = firstn; firstn->NZE_C_N = first; } @@ -475,13 +471,13 @@ dynSparseMatrix::Simple_Init(int Size, const map, int> &IM, first->r_index = eq; first->c_index = var; first->lag_index = lag; - if (FNZE_R[eq] == nullptr) + if (!FNZE_R[eq]) FNZE_R[eq] = first; - if (FNZE_C[var] == nullptr) + if (!FNZE_C[var]) FNZE_C[var] = first; - if (temp_NZE_R[eq] != nullptr) + if (temp_NZE_R[eq]) temp_NZE_R[eq]->NZE_R_N = first; - if (temp_NZE_C[var] != nullptr) + if (temp_NZE_C[var]) temp_NZE_C[var]->NZE_C_N = first; temp_NZE_R[eq] = first; temp_NZE_C[var] = first; @@ -1238,9 +1234,9 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, long j = 0, i = 0; while (i < nop4 && OK) { - t_save_op_s *save_op_s = reinterpret_cast(&(save_op[i])); - t_save_op_s *save_opa_s = reinterpret_cast(&(save_opa[i])); - t_save_op_s *save_opaa_s = reinterpret_cast(&(save_opaa[i])); + t_save_op_s *save_op_s = reinterpret_cast(&save_op[i]); + t_save_op_s *save_opa_s = reinterpret_cast(&save_opa[i]); + t_save_op_s *save_opaa_s = reinterpret_cast(&save_opaa[i]); diff1[j] = save_op_s->first-save_opa_s->first; max_save_ops_first = max(max_save_ops_first, save_op_s->first+diff1[j]*(periods-beg_t)); switch (save_op_s->operat) @@ -1284,7 +1280,7 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int i = j = 0; while (i < nop4) { - auto *save_op_s = reinterpret_cast(&(save_op[i])); + auto *save_op_s = reinterpret_cast(&save_op[i]); double *up = &u[save_op_s->first+t*diff1[j]]; switch (save_op_s->operat) { @@ -1316,7 +1312,7 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int gap = periods_beg_t-t; while (i < nop4) { - if (auto *save_op_s = reinterpret_cast(&(save_op[i])); + if (auto *save_op_s = reinterpret_cast(&save_op[i]); save_op_s->lag < gap) { double *up = &u[save_op_s->first+t*diff1[j]]; @@ -2868,14 +2864,13 @@ dynSparseMatrix::Singular_display(int block, int Size) { bool zero_solution; Simple_Init(Size, IM_i, zero_solution); - NonZeroElem *first; mxArray *rhs[1] = { mxCreateDoubleMatrix(Size, Size, mxREAL) }; - double *pind; - pind = mxGetPr(rhs[0]); + double *pind = mxGetPr(rhs[0]); for (int j = 0; j < Size * Size; j++) pind[j] = 0.0; for (int ii = 0; ii < Size; ii++) { + NonZeroElem *first; int nb_eq = At_Col(ii, &first); for (int j = 0; j < nb_eq; j++) { @@ -2889,96 +2884,78 @@ dynSparseMatrix::Singular_display(int block, int Size) mexCallMATLAB(3, lhs, 1, rhs, "svd"); mxArray *SVD_u = lhs[0]; mxArray *SVD_s = lhs[1]; - //mxArray* SVD_v = lhs[2]; double *SVD_ps = mxGetPr(SVD_s); double *SVD_pu = mxGetPr(SVD_u); for (int i = 0; i < Size; i++) - { - if (abs(SVD_ps[i * (1 + Size)]) < 1e-12) - { - mexPrintf(" The following equations form a linear combination:\n "); - double max_u = 0; - for (int j = 0; j < Size; j++) - if (abs(SVD_pu[j + i * Size]) > abs(max_u)) - max_u = SVD_pu[j + i * Size]; - vector equ_list; - for (int j = 0; j < Size; j++) - { - double rr = SVD_pu[j + i * Size] / max_u; - if (rr < -1e-10) - { - equ_list.push_back(j); - if (rr != -1) - mexPrintf(" - %3.2f*Dequ_%d_dy", abs(rr), j+1); + if (abs(SVD_ps[i * (1 + Size)]) < 1e-12) + { + mexPrintf(" The following equations form a linear combination:\n "); + double max_u = 0; + for (int j = 0; j < Size; j++) + if (abs(SVD_pu[j + i * Size]) > abs(max_u)) + max_u = SVD_pu[j + i * Size]; + vector equ_list; + for (int j = 0; j < Size; j++) + { + double rr = SVD_pu[j + i * Size] / max_u; + if (rr < -1e-10) + { + equ_list.push_back(j); + if (rr != -1) + mexPrintf(" - %3.2f*Dequ_%d_dy", abs(rr), j+1); + else + mexPrintf(" - Dequ_%d_dy", j+1); + } + else if (rr > 1e-10) + { + equ_list.push_back(j); + if (j > 0) + if (rr != 1) + mexPrintf(" + %3.2f*Dequ_%d_dy", rr, j+1); else - mexPrintf(" - Dequ_%d_dy", j+1); - } - else if (rr > 1e-10) - { - equ_list.push_back(j); - if (j > 0) - if (rr != 1) - mexPrintf(" + %3.2f*Dequ_%d_dy", rr, j+1); - else - mexPrintf(" + Dequ_%d_dy", j+1); - else if (rr != 1) - mexPrintf(" %3.2f*Dequ_%d_dy", rr, j+1); - else - mexPrintf(" Dequ_%d_dy", j+1); - } - } - mexPrintf(" = 0\n"); - /*mexPrintf(" with:\n"); - it_code = get_begin_block(block); - for (int j=0; j < Size; j++) - { - if (find(equ_list.begin(), equ_list.end(), j) != equ_list.end()) - mexPrintf(" equ_%d: %s\n",j, print_expression(it_code_expr, false, Size, block, steady_state, 0, 0, it_code, true).c_str()); - }*/ - } - } + mexPrintf(" + Dequ_%d_dy", j+1); + else if (rr != 1) + mexPrintf(" %3.2f*Dequ_%d_dy", rr, j+1); + else + mexPrintf(" Dequ_%d_dy", j+1); + } + } + mexPrintf(" = 0\n"); + } mxDestroyArray(lhs[0]); mxDestroyArray(lhs[1]); mxDestroyArray(lhs[2]); - ostringstream tmp; if (block > 1) - tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " << block+1 << "\n"; + throw FatalExceptionHandling(" in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " + + to_string(block+1) + "\n"); else - tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system\n"; - throw FatalExceptionHandling(tmp.str()); + throw FatalExceptionHandling(" in Solve_ByteCode_Sparse_GaussianElimination, singular system\n"); } bool dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, int it_) { - bool one; int pivj = 0, pivk = 0; - double *piv_v; - int *pivj_v, *pivk_v, *NR; - int l, N_max; - NonZeroElem *first, *firsta, *first_suba; - double piv_abs; - NonZeroElem **bc; - bc = static_cast(mxMalloc(Size*sizeof(*bc))); + NonZeroElem **bc = static_cast(mxMalloc(Size*sizeof(*bc))); test_mxMalloc(bc, __LINE__, __FILE__, __func__, Size*sizeof(*bc)); - piv_v = static_cast(mxMalloc(Size*sizeof(double))); + double *piv_v = static_cast(mxMalloc(Size*sizeof(double))); test_mxMalloc(piv_v, __LINE__, __FILE__, __func__, Size*sizeof(double)); - pivj_v = static_cast(mxMalloc(Size*sizeof(int))); + int *pivj_v = static_cast(mxMalloc(Size*sizeof(int))); test_mxMalloc(pivj_v, __LINE__, __FILE__, __func__, Size*sizeof(int)); - pivk_v = static_cast(mxMalloc(Size*sizeof(int))); + int *pivk_v = static_cast(mxMalloc(Size*sizeof(int))); test_mxMalloc(pivk_v, __LINE__, __FILE__, __func__, Size*sizeof(int)); - NR = static_cast(mxMalloc(Size*sizeof(int))); + int *NR = static_cast(mxMalloc(Size*sizeof(int))); test_mxMalloc(NR, __LINE__, __FILE__, __func__, Size*sizeof(int)); for (int i = 0; i < Size; i++) { /*finding the max-pivot*/ - double piv = piv_abs = 0; + double piv = 0, piv_abs = 0; + NonZeroElem *first; int nb_eq = At_Col(i, &first); - l = 0; - N_max = 0; - one = false; - piv_abs = 0; + int l = 0; + int N_max = 0; + bool one = false; for (int j = 0; j < nb_eq; j++) { if (!line_done[first->r_index]) @@ -3005,15 +2982,12 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i if (NRow_jj > N_max) N_max = NRow_jj; } - else + else if (NRow_jj == 1) { - if (NRow_jj == 1) - { - if (piv_fabs > piv_abs) - piv_abs = piv_fabs; - if (NRow_jj > N_max) - N_max = NRow_jj; - } + if (piv_fabs > piv_abs) + piv_abs = piv_fabs; + if (NRow_jj > N_max) + N_max = NRow_jj; } l++; } @@ -3036,27 +3010,27 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i } else { - ostringstream tmp; if (blck > 1) - tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " << blck+1 << "\n"; + throw FatalExceptionHandling(" in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " + + to_string(blck+1) + "\n"); else - tmp << " in Solve_ByteCode_Sparse_GaussianElimination, singular system\n"; - throw FatalExceptionHandling(tmp.str()); + throw FatalExceptionHandling(" in Solve_ByteCode_Sparse_GaussianElimination, singular system\n"); } } double markovitz = 0, markovitz_max = -9e70; if (!one) - { - for (int j = 0; j < l; j++) - { - if (N_max > 0 && NR[j] > 0) - { - if (fabs(piv_v[j]) > 0) - { - if (markowitz_c > 0) - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(static_cast(NR[j])/static_cast(N_max))); - else - markovitz = fabs(piv_v[j])/piv_abs; + for (int j = 0; j < l; j++) + { + if (N_max > 0 && NR[j] > 0) + { + if (fabs(piv_v[j]) > 0) + { + if (markowitz_c > 0) + markovitz = exp(log(fabs(piv_v[j])/piv_abs) + -markowitz_c*log(static_cast(NR[j]) + /static_cast(N_max))); + else + markovitz = fabs(piv_v[j])/piv_abs; } else markovitz = 0; @@ -3070,35 +3044,34 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i pivk = pivk_v[j]; //positi markovitz_max = markovitz; } - } - } + } else - { - for (int j = 0; j < l; j++) - { - if (N_max > 0 && NR[j] > 0) - { - if (fabs(piv_v[j]) > 0) - { - if (markowitz_c > 0) - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(static_cast(NR[j])/static_cast(N_max))); - else - markovitz = fabs(piv_v[j])/piv_abs; - } - else - markovitz = 0; - } - else - markovitz = fabs(piv_v[j])/piv_abs; - if (NR[j] == 1) - { - piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi - markovitz_max = markovitz; - } - } - } + for (int j = 0; j < l; j++) + { + if (N_max > 0 && NR[j] > 0) + { + if (fabs(piv_v[j]) > 0) + { + if (markowitz_c > 0) + markovitz = exp(log(fabs(piv_v[j])/piv_abs) + -markowitz_c*log(static_cast(NR[j]) + /static_cast(N_max))); + else + markovitz = fabs(piv_v[j])/piv_abs; + } + else + markovitz = 0; + } + else + markovitz = fabs(piv_v[j])/piv_abs; + if (NR[j] == 1) + { + piv = piv_v[j]; + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi + markovitz_max = markovitz; + } + } pivot[i] = pivj; pivotk[i] = pivk; pivotv[i] = piv; @@ -3123,7 +3096,6 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i bc[nb_eq_todo++] = first; first = first->NZE_C_N; } - //pragma omp parallel for for (int j = 0; j < nb_eq_todo; j++) { first = bc[j]; @@ -3137,67 +3109,65 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i int l_sub = 0, l_piv = 0; int sub_c_index = first_sub->c_index, piv_c_index = first_piv->c_index; while (l_sub < nb_var_sub || l_piv < nb_var_piv) - { - if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv)) - { - first_sub = first_sub->NZE_R_N; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size; - l_sub++; - } - else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub) - { - int tmp_u_count = Get_u(); - Insert(row, first_piv->c_index, tmp_u_count, 0); - u[tmp_u_count] = -u[first_piv->u_index]*first_elem; - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size; - l_piv++; - } - else - { - if (i == sub_c_index) - { - firsta = first; - first_suba = first_sub->NZE_R_N; - Delete(first_sub->r_index, first_sub->c_index); - first = firsta->NZE_C_N; - first_sub = first_suba; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size; - l_sub++; - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size; - l_piv++; - } - else - { - u[first_sub->u_index] -= u[first_piv->u_index]*first_elem; - first_sub = first_sub->NZE_R_N; - if (first_sub) - sub_c_index = first_sub->c_index; - else - sub_c_index = Size; - l_sub++; - first_piv = first_piv->NZE_R_N; - if (first_piv) - piv_c_index = first_piv->c_index; - else - piv_c_index = Size; - l_piv++; - } - } - } + if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv)) + { + first_sub = first_sub->NZE_R_N; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size; + l_sub++; + } + else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub) + { + int tmp_u_count = Get_u(); + Insert(row, first_piv->c_index, tmp_u_count, 0); + u[tmp_u_count] = -u[first_piv->u_index]*first_elem; + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size; + l_piv++; + } + else + { + if (i == sub_c_index) + { + NonZeroElem *firsta = first; + NonZeroElem *first_suba = first_sub->NZE_R_N; + Delete(first_sub->r_index, first_sub->c_index); + first = firsta->NZE_C_N; + first_sub = first_suba; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size; + l_sub++; + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size; + l_piv++; + } + else + { + u[first_sub->u_index] -= u[first_piv->u_index]*first_elem; + first_sub = first_sub->NZE_R_N; + if (first_sub) + sub_c_index = first_sub->c_index; + else + sub_c_index = Size; + l_sub++; + first_piv = first_piv->NZE_R_N; + if (first_piv) + piv_c_index = first_piv->c_index; + else + piv_c_index = Size; + l_piv++; + } + } u[b[row]] -= u[b[pivj]]*first_elem; } } @@ -3220,35 +3190,26 @@ void dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number) { /*Triangularisation at each period of a block using a simple gaussian Elimination*/ - t_save_op_s *save_op_s; int *save_op = nullptr, *save_opa = nullptr, *save_opaa = nullptr; long int nop = 0, nopa = 0; bool record = false; - double *piv_v; - double piv_abs; - int *pivj_v, *pivk_v, *NR; + int pivj = 0, pivk = 0; - NonZeroElem *first; - int tmp_u_count, lag; int tbreak = 0, last_period = periods; - piv_v = static_cast(mxMalloc(Size*sizeof(double))); + double *piv_v = static_cast(mxMalloc(Size*sizeof(double))); test_mxMalloc(piv_v, __LINE__, __FILE__, __func__, Size*sizeof(double)); - pivj_v = static_cast(mxMalloc(Size*sizeof(int))); + int *pivj_v = static_cast(mxMalloc(Size*sizeof(int))); test_mxMalloc(pivj_v, __LINE__, __FILE__, __func__, Size*sizeof(int)); - pivk_v = static_cast(mxMalloc(Size*sizeof(int))); + int *pivk_v = static_cast(mxMalloc(Size*sizeof(int))); test_mxMalloc(pivk_v, __LINE__, __FILE__, __func__, Size*sizeof(int)); - NR = static_cast(mxMalloc(Size*sizeof(int))); + int *NR = static_cast(mxMalloc(Size*sizeof(int))); test_mxMalloc(NR, __LINE__, __FILE__, __func__, Size*sizeof(int)); - //clock_t time00 = clock(); - NonZeroElem **bc; - bc = static_cast(mxMalloc(Size*sizeof(first))); - test_mxMalloc(bc, __LINE__, __FILE__, __func__, Size*sizeof(first)); + NonZeroElem **bc = static_cast(mxMalloc(Size*sizeof(NonZeroElem *))); + test_mxMalloc(bc, __LINE__, __FILE__, __func__, Size*sizeof(NonZeroElem *)); for (int t = 0; t < periods; t++) { - /*clock_t time11 = clock(); - mexPrintf("t=%d, record = %d\n",t, record);*/ #ifdef MATLAB_MEX_FILE if (utIsInterruptPending()) throw UserExceptionHandling(); @@ -3256,11 +3217,6 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo if (record && symbolic) { - /*if (save_op) - { - mxFree(save_op); - save_op = NULL; - }*/ save_op = static_cast(mxMalloc(nop*sizeof(int))); test_mxMalloc(save_op, __LINE__, __FILE__, __func__, nop*sizeof(int)); nopa = nop; @@ -3271,7 +3227,8 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo for (int i = ti; i < Size+ti; i++) { /*finding the max-pivot*/ - double piv = piv_abs = 0; + double piv = 0, piv_abs = 0; + NonZeroElem *first; int nb_eq = At_Col(i, 0, &first); if ((symbolic && t <= start_compare) || !symbolic) { @@ -3303,15 +3260,12 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo if (NRow_jj > N_max) N_max = NRow_jj; } - else + else if (NRow_jj == 1) { - if (NRow_jj == 1) - { - if (piv_fabs > piv_abs) - piv_abs = piv_fabs; - if (NRow_jj > N_max) - N_max = NRow_jj; - } + if (piv_fabs > piv_abs) + piv_abs = piv_fabs; + if (NRow_jj > N_max) + N_max = NRow_jj; } l++; } @@ -3320,61 +3274,61 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo double markovitz = 0, markovitz_max = -9e70; int NR_max = 0; if (!one) - { - for (int j = 0; j < l; j++) - { - if (N_max > 0 && NR[j] > 0) - { - if (fabs(piv_v[j]) > 0) - { - if (markowitz_c > 0) - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(static_cast(NR[j])/static_cast(N_max))); - else - markovitz = fabs(piv_v[j])/piv_abs; - } - else - markovitz = 0; - } - else - markovitz = fabs(piv_v[j])/piv_abs; - if (markovitz > markovitz_max) - { - piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi - markovitz_max = markovitz; - NR_max = NR[j]; - } - } - } + for (int j = 0; j < l; j++) + { + if (N_max > 0 && NR[j] > 0) + { + if (fabs(piv_v[j]) > 0) + { + if (markowitz_c > 0) + markovitz = exp(log(fabs(piv_v[j])/piv_abs) + -markowitz_c*log(static_cast(NR[j]) + /static_cast(N_max))); + else + markovitz = fabs(piv_v[j])/piv_abs; + } + else + markovitz = 0; + } + else + markovitz = fabs(piv_v[j])/piv_abs; + if (markovitz > markovitz_max) + { + piv = piv_v[j]; + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi + markovitz_max = markovitz; + NR_max = NR[j]; + } + } else - { - for (int j = 0; j < l; j++) - { - if (N_max > 0 && NR[j] > 0) - { - if (fabs(piv_v[j]) > 0) - { - if (markowitz_c > 0) - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(static_cast(NR[j])/static_cast(N_max))); - else - markovitz = fabs(piv_v[j])/piv_abs; - } - else - markovitz = 0; - } - else - markovitz = fabs(piv_v[j])/piv_abs; - if (NR[j] == 1) - { - piv = piv_v[j]; - pivj = pivj_v[j]; //Line number - pivk = pivk_v[j]; //positi - markovitz_max = markovitz; - NR_max = NR[j]; - } - } - } + for (int j = 0; j < l; j++) + { + if (N_max > 0 && NR[j] > 0) + { + if (fabs(piv_v[j]) > 0) + { + if (markowitz_c > 0) + markovitz = exp(log(fabs(piv_v[j])/piv_abs) + -markowitz_c*log(static_cast(NR[j]) + /static_cast(N_max))); + else + markovitz = fabs(piv_v[j])/piv_abs; + } + else + markovitz = 0; + } + else + markovitz = fabs(piv_v[j])/piv_abs; + if (NR[j] == 1) + { + piv = piv_v[j]; + pivj = pivj_v[j]; //Line number + pivk = pivk_v[j]; //positi + markovitz_max = markovitz; + NR_max = NR[j]; + } + } if (fabs(piv) < eps) mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n", NR_max, N_max, piv, piv_abs, markovitz_max); if (NR_max == 0) @@ -3399,22 +3353,21 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo { if (nop+1 >= nopa) { - nopa = long (mem_increasing_factor*static_cast(nopa)); + nopa = static_cast(mem_increasing_factor*static_cast(nopa)); save_op = static_cast(mxRealloc(save_op, nopa*sizeof(int))); } - save_op_s = reinterpret_cast(&(save_op[nop])); + t_save_op_s *save_op_s = reinterpret_cast(&save_op[nop]); save_op_s->operat = IFLD; save_op_s->first = pivk; save_op_s->lag = 0; nop += 2; if (piv_abs < eps) { - ostringstream tmp; if (Block_number > 1) - tmp << " in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system in block " << Block_number+1 << "\n"; + throw FatalExceptionHandling(" in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system in block " + + to_string(Block_number+1) + "\n"); else - tmp << " in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system\n"; - throw FatalExceptionHandling(tmp.str()); + throw FatalExceptionHandling(" in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system\n"); } /*divide all the non zeros elements of the line pivj by the max_pivot*/ int nb_var = At_Row(pivj, &first); @@ -3426,7 +3379,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo nopa = static_cast(mem_increasing_factor*static_cast(nopa)); save_op = static_cast(mxRealloc(save_op, nopa*sizeof(int))); } - save_op_s = reinterpret_cast(&(save_op[nop+j*2])); + save_op_s = reinterpret_cast(&save_op[nop+j*2]); save_op_s->operat = IFDIV; save_op_s->first = first->u_index; save_op_s->lag = first->lag_index; @@ -3439,12 +3392,12 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo nopa = static_cast(mem_increasing_factor*static_cast(nopa)); save_op = static_cast(mxRealloc(save_op, nopa*sizeof(int))); } - save_op_s = reinterpret_cast(&(save_op[nop])); + save_op_s = reinterpret_cast(&save_op[nop]); save_op_s->operat = IFDIV; save_op_s->first = b[pivj]; save_op_s->lag = 0; nop += 2; - /*substract the elements on the non treated lines*/ + // Substract the elements on the non treated lines nb_eq = At_Col(i, &first); NonZeroElem *first_piva; int nb_var_piva = At_Row(pivj, &first_piva); @@ -3456,7 +3409,6 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo bc[nb_eq_todo++] = first; first = first->NZE_C_N; } - //#pragma omp parallel for shared(nb_var_piva, first_piva, nopa, save_op) reduction(+:nop) for (int j = 0; j < nb_eq_todo; j++) { t_save_op_s *save_op_s_l; @@ -3468,7 +3420,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo nopa = static_cast(mem_increasing_factor*static_cast(nopa)); save_op = static_cast(mxRealloc(save_op, nopa*sizeof(int))); } - save_op_s_l = reinterpret_cast(&(save_op[nop])); + save_op_s_l = reinterpret_cast(&save_op[nop]); save_op_s_l->operat = IFLD; save_op_s_l->first = first->u_index; save_op_s_l->lag = abs(first->lag_index); @@ -3483,11 +3435,13 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo int sub_c_index = first_sub->c_index; int piv_c_index = first_piv->c_index; int tmp_lag = first_sub->lag_index; - while (l_sub < (nb_var_sub /*=NRow(row)*/) || l_piv < nb_var_piv) + while (l_sub < nb_var_sub /*=NRow(row)*/ || l_piv < nb_var_piv) { if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv)) { - //There is no nonzero element at row pivot for this column=> Nothing to do for the current element got to next column + /* There is no nonzero element at row pivot for this + column ⇒ Nothing to do for the current element got to + next column */ first_sub = first_sub->NZE_R_N; if (first_sub) sub_c_index = first_sub->c_index; @@ -3498,19 +3452,16 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub) { // There is an nonzero element at row pivot but not at the current row=> insert a negative element in the current row - tmp_u_count = Get_u(); - lag = first_piv->c_index/Size-row/Size; - //#pragma omp critical - { - Insert(row, first_piv->c_index, tmp_u_count, lag); - } + int tmp_u_count = Get_u(); + int lag = first_piv->c_index/Size-row/Size; + Insert(row, first_piv->c_index, tmp_u_count, lag); u[tmp_u_count] = -u[first_piv->u_index]*first_elem; if (nop+2 >= nopa) { nopa = static_cast(mem_increasing_factor*static_cast(nopa)); save_op = static_cast(mxRealloc(save_op, nopa*sizeof(int))); } - save_op_s_l = reinterpret_cast(&(save_op[nop])); + save_op_s_l = reinterpret_cast(&save_op[nop]); save_op_s_l->operat = IFLESS; save_op_s_l->first = tmp_u_count; save_op_s_l->second = first_piv->u_index; @@ -3529,10 +3480,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo { NonZeroElem *firsta = first; NonZeroElem *first_suba = first_sub->NZE_R_N; - //#pragma omp critical - { - Delete(first_sub->r_index, first_sub->c_index); - } + Delete(first_sub->r_index, first_sub->c_index); first = firsta->NZE_C_N; first_sub = first_suba; if (first_sub) @@ -3555,7 +3503,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo nopa = static_cast(mem_increasing_factor*static_cast(nopa)); save_op = static_cast(mxRealloc(save_op, nopa*sizeof(int))); } - save_op_s_l = reinterpret_cast(&(save_op[nop])); + save_op_s_l = reinterpret_cast(&save_op[nop]); save_op_s_l->operat = IFSUB; save_op_s_l->first = first_sub->u_index; save_op_s_l->second = first_piv->u_index; @@ -3583,7 +3531,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo nopa = static_cast(mem_increasing_factor*static_cast(nopa)); save_op = static_cast(mxRealloc(save_op, nopa*sizeof(int))); } - save_op_s_l = reinterpret_cast(&(save_op[nop])); + save_op_s_l = reinterpret_cast(&save_op[nop]); save_op_s_l->operat = IFSUB; save_op_s_l->first = b[row]; save_op_s_l->second = b[pivj]; @@ -3596,14 +3544,13 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo nop += 2; if (piv_abs < eps) { - ostringstream tmp; if (Block_number > 1) - tmp << " in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system in block " << Block_number+1 << "\n"; + throw FatalExceptionHandling(" in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system in block " + + to_string(Block_number+1) + "\n"); else - tmp << " in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system\n"; - throw FatalExceptionHandling(tmp.str()); + throw FatalExceptionHandling(" in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system\n"); } - /*divide all the non zeros elements of the line pivj by the max_pivot*/ + // Divide all the non zeros elements of the line pivj by the max_pivot int nb_var = At_Row(pivj, &first); for (int j = 0; j < nb_var; j++) { @@ -3613,7 +3560,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo nop += nb_var*2; u[b[pivj]] /= piv; nop += 2; - /*substract the elements on the non treated lines*/ + // Substract the elements on the non treated lines nb_eq = At_Col(i, &first); NonZeroElem *first_piva; int nb_var_piva = At_Row(pivj, &first_piva); @@ -3625,7 +3572,6 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo bc[nb_eq_todo++] = first; first = first->NZE_C_N; } - //#pragma omp parallel for shared(nb_var_piva, first_piva, nopa, save_op) reduction(+:nop) for (int j = 0; j < nb_eq_todo; j++) { NonZeroElem *first = bc[j]; @@ -3640,11 +3586,13 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo int l_piv = 0; int sub_c_index = first_sub->c_index; int piv_c_index = first_piv->c_index; - while (l_sub < (nb_var_sub /*= NRow(row)*/) || l_piv < nb_var_piv) + while (l_sub < nb_var_sub /*= NRow(row)*/ || l_piv < nb_var_piv) { if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv)) { - //There is no nonzero element at row pivot for this column=> Nothing to do for the current element got to next column + /* There is no nonzero element at row pivot for this + column ⇒ Nothing to do for the current element got to + next column */ first_sub = first_sub->NZE_R_N; if (first_sub) sub_c_index = first_sub->c_index; @@ -3654,13 +3602,12 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo } else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub) { - // There is an nonzero element at row pivot but not at the current row=> insert a negative element in the current row - tmp_u_count = Get_u(); - lag = first_piv->c_index/Size-row/Size; - //#pragma omp critical - { - Insert(row, first_piv->c_index, tmp_u_count, lag); - } + /* There is an nonzero element at row pivot but not + at the current row ⇒ insert a negative element in the + current row */ + int tmp_u_count = Get_u(); + int lag = first_piv->c_index/Size-row/Size; + Insert(row, first_piv->c_index, tmp_u_count, lag); u[tmp_u_count] = -u[first_piv->u_index]*first_elem; nop += 3; first_piv = first_piv->NZE_R_N; @@ -3676,10 +3623,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo { NonZeroElem *firsta = first; NonZeroElem *first_suba = first_sub->NZE_R_N; - //#pragma omp critical - { - Delete(first_sub->r_index, first_sub->c_index); - } + Delete(first_sub->r_index, first_sub->c_index); first = firsta->NZE_C_N; first_sub = first_suba; if (first_sub) @@ -3727,7 +3671,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo mxFree(save_opa); mxFree(save_op); } - else if (record && (nop == nop1)) + else if (record && nop == nop1) { if (t > static_cast(periods*0.35)) { @@ -3754,7 +3698,6 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo { tbreak = t; tbreak_g = tbreak; - //mexPrintf("time=%f\n",(1000.0*(static_cast(clock())-static_cast(time11)))/static_cast(CLOCKS_PER_SEC)); break; } } @@ -3791,16 +3734,12 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo nop2 = nop1; nop1 = nop; } - //mexPrintf("time=%f\n",(1000.0*(static_cast(clock())-static_cast(time11)))/static_cast(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*(static_cast(clock())-static_cast(time00)))/static_cast(CLOCKS_PER_SEC)); - mexEvalString("drawnow;"); - time00 = clock();*/ nop_all += nop; if (symbolic) { @@ -3812,84 +3751,21 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo mxFree(save_opaa); } - /*The backward substitution*/ + // The backward substitution double slowc_lbx = slowc; for (int i = 0; i < y_size*(periods+y_kmin); i++) 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*(static_cast(clock())-static_cast(time00)))/static_cast(CLOCKS_PER_SEC)); - mexEvalString("drawnow;");*/ End_GE(Size); } -void -dynSparseMatrix::Grad_f_product(int n, mxArray *b_m, double *vectr, mxArray *A_m, SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b_) -{ - if ((solve_algo == 5 && steady_state) || (stack_solve_algo == 5 && !steady_state)) - { - NonZeroElem *first; - for (int i = 0; i < n; i++) - { - double sum = 0; - first = FNZE_R[i]; - if (first) - for (int k = 0; k < NbNZRow[i]; k++) - { - sum += u[first->u_index] * u[b[first->c_index]]; - first = first->NZE_R_N; - } - vectr[i] = sum; - } - } - else - { - if (!((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state))) - { - 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()); - } - 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()); - } - 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()); - } - b_ = mxGetPr(b_m); - if (!b_) - { - ostringstream tmp; - tmp << " in Init_Matlab_Sparse_Simple, can't retrieve b matrix\n"; - throw FatalExceptionHandling(tmp.str()); - } - } - memset(vectr, 0, n * sizeof(double)); - for (int i = 0; i < n; i++) - for (SuiteSparse_long j = Ap[i]; j < Ap[i+1]; j++) - vectr[Ai[j]] += Ax[j] * b_[i]; - } -} - void dynSparseMatrix::Check_and_Correct_Previous_Iteration(int block_num, int y_size, int size, double crit_opt_old) { - double top = 1.0; - double bottom = 0.1; if (isnan(res1) || isinf(res1) || (res2 > g0 && iter > 0)) { - while ((isnan(res1) || isinf(res1))) + while (isnan(res1) || isinf(res1)) { prev_slowc_save = slowc_save; slowc_save /= 1.1; @@ -3898,9 +3774,6 @@ dynSparseMatrix::Check_and_Correct_Previous_Iteration(int block_num, int y_size, int eq = index_vara[i]; y[eq+it_*y_size] = ya[eq+it_*y_size] + slowc_save * direction[eq+it_*y_size]; } - /*mexPrintf("reducing solwc_save = %e, it_=%d, y_size=%d, size=%d, y[%d]=%e, ya[%d]=%e,\n y[%d]=%e, ya[%d]=%e\n",slowc_save, it_, y_size, size-1, index_vara[0]+it_*y_size, y[index_vara[0]+it_*y_size], index_vara[0]+it_*y_size, ya[index_vara[0]+it_*y_size] - , index_vara[size-1]+it_*y_size, y[index_vara[size-1]+it_*y_size], index_vara[size-1]+it_*y_size, ya[index_vara[size-1]+it_*y_size]);*/ - //mexPrintf("->slowc_save=%f\n",slowc_save); compute_complete(true, res1, res2, max_res, max_res_idx); } @@ -3913,92 +3786,8 @@ dynSparseMatrix::Check_and_Correct_Previous_Iteration(int block_num, int y_size, int eq = index_vara[i]; y[eq+it_*y_size] = ya[eq+it_*y_size] + slowc_save * direction[eq+it_*y_size]; } - /*mexPrintf("reducing solwc_save = %e, it_=%d, y_size=%d, size=%d, y[%d]=%e, ya[%d]=%e,\n y[%d]=%e, ya[%d]=%e\n",slowc_save, it_, y_size, size-1, index_vara[0]+it_*y_size, y[index_vara[0]+it_*y_size], index_vara[0]+it_*y_size, ya[index_vara[0]+it_*y_size] , index_vara[size-1]+it_*y_size, y[index_vara[size-1]+it_*y_size], index_vara[size-1]+it_*y_size, ya[index_vara[size-1]+it_*y_size]);*/ - //mexPrintf("->slowc_save=%f\n",slowc_save); compute_complete(true, res1, res2, max_res, max_res_idx); } - double ax = slowc_save-0.001, bx = slowc_save+0.001, cx = slowc_save, fa, fb, fc, xmin; - if (false /*slowc_save > 2e-1*/) - if (mnbrak(&ax, &bx, &cx, &fa, &fb, &fc)) - if (golden(ax, bx, cx, 1e-1, solve_tolf, &xmin)) - slowc_save = xmin; - //mexPrintf("cx=%f\n", cx); - //mexPrintf("ax= %f, bx=%f, cx=%f, fa=%f, fb=%f, fc=%d\n", ax, bx, cx, fa, fb, fc); - - //if (!(isnan(res1) || isinf(res1))/* && !(isnan(g0) || isinf(g0))*//*|| (res2 > g0 && iter > 1)*/) - if (false) - { - - auto *p = static_cast(mxMalloc(size * sizeof(double))); - test_mxMalloc(p, __LINE__, __FILE__, __func__, size * sizeof(double)); - Grad_f_product(size, b_m_save, p, A_m_save, Ap_save, Ai_save, Ax_save, b_save); - double slope = 0.0; - for (int i = 1; i < size; i++) - slope += -direction[i] * p[i]; - /*if (slope > 0) - mexPrintf("Roundoff in lnsearch\n"); - else*/ - { - prev_slowc_save = 1; - double crit_opt = res2/2; - double max_try_iteration = 100; - double small_ = 1.0e-4; - bool try_at_cvg = false; - while ((try_at_iteration < max_try_iteration) && (!try_at_cvg) && (abs(prev_slowc_save - slowc_save) > 1e-10)) - { - crit_opt = res2 / 2; - if (slowc_save < 1e-7) - { - try_at_cvg = true; - continue; - } - else if ((crit_opt <= crit_opt_old + small_ * slowc_save * slope) && !(isnan(res1) || isinf(res1))) - { - try_at_cvg = true; - continue; - } - else if (try_at_iteration == 0) - { - prev_slowc_save = slowc_save; - //slowc_save = max(- top * slope / ( (crit_opt - crit_opt_old - slope)), bottom); - slowc_save /= 1.2; - } - else - { - double t1 = crit_opt - slope * slowc_save - crit_opt_old; - double t2 = glambda2 - slope * prev_slowc_save - crit_opt_old; - double a = (1/(slowc_save * slowc_save) * t1 - 1/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save); - double b = (-prev_slowc_save/(slowc_save * slowc_save) * t1 + slowc_save/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save); - if (a == 0) - slowc_save = max(min(-slope/(2 * b), top * slowc_save), bottom * slowc_save); - else - { - double delta = b*b - 3 * a * slope; - if (delta <= 0) - slowc_save = top * slowc_save; - else if (b <= 0) - slowc_save = max(min(-b + sqrt(delta) / (3 * a), top * slowc_save), bottom * slowc_save); - else - slowc_save = max(min(-slope / (b + sqrt(delta)), top * slowc_save), bottom * slowc_save); - } - } - if (abs(prev_slowc_save - slowc_save) < 1e-10) - slowc_save /= 1.1; - //mexPrintf("=>slowc_save=%f, prev_slowc_save=%f\n",slowc_save, prev_slowc_save); - prev_slowc_save = slowc_save; - glambda2 = crit_opt; - try_at_iteration++; - for (int i = 0; i < size; i++) - { - int eq = index_vara[i]; - y[eq+it_*y_size] = ya[eq+it_*y_size] + slowc_save * direction[eq+it_*y_size]; - } - compute_complete(true, res1, res2, max_res, max_res_idx); - } - } - mxFree(p); - } - //if (print_it) mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save); for (int i = 0; i < size; i++) { @@ -4008,21 +3797,17 @@ dynSparseMatrix::Check_and_Correct_Previous_Iteration(int block_num, int y_size, compute_complete(false, res1, res2, max_res, max_res_idx); } else - { - //mexPrintf("slowc_save=%f res1=%f\n",slowc_save, res1); - for (int i = 0; i < size; i++) - { - int eq = index_vara[i]; - y[eq+it_*y_size] = ya[eq+it_*y_size] + slowc_save * 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] + slowc_save * direction[eq+it_*y_size]; + } slowc_save = slowc; } bool dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, int y_kmax, int size, bool cvg) { - //int i; mxArray *b_m = nullptr, *A_m = nullptr, *x0_m = nullptr; SuiteSparse_long *Ap = nullptr, *Ai = nullptr; double *Ax = nullptr, *b = nullptr; @@ -4058,16 +3843,16 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in else mexPrintf(" dynare cannot improve the simulation in block %d at time %d (variable %d)\n", block_num+1, it_+1, index_vara[max_res_idx]+1); mexEvalString("drawnow;"); - //return singular_system; } else { - ostringstream tmp; if (iter == 0) - tmp << " in Simulate_One_Boundary, The initial values of endogenous variables are too far from the solution.\nChange them!\n"; + throw FatalExceptionHandling(" in Simulate_One_Boundary, The initial values of endogenous variables are too far from the solution.\nChange them!\n"); else - tmp << " in Simulate_One_Boundary, Dynare cannot improve the simulation in block " << block_num+1 << " at time " << it_+1 << " (variable " << index_vara[max_res_idx]+1 << "%d)\n"; - throw FatalExceptionHandling(tmp.str()); + throw FatalExceptionHandling(" in Simulate_One_Boundary, Dynare cannot improve the simulation in block " + + to_string(block_num+1) + " at time " + to_string(it_+1) + + " (variable " + to_string(index_vara[max_res_idx]+1) + + "%d)\n"); } } if (print_it) @@ -4097,11 +3882,9 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in break; case 7: mexPrintf(preconditioner_print_out("MODEL STEADY STATE: (method=GMRES)\n", preconditioner, true).c_str()); - //mexPrintf("MODEL STEADY STATE: (method=GMRES)\n"); break; case 8: mexPrintf(preconditioner_print_out("MODEL STEADY STATE: (method=BiCGStab)\n", preconditioner, true).c_str()); - //mexPrintf("MODEL STEADY STATE: (method=BiCGStab)\n"); break; default: mexPrintf("MODEL STEADY STATE: (method=Unknown - %d - )\n", stack_solve_algo); @@ -4123,26 +3906,15 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in { b_m = mxCreateDoubleMatrix(size, 1, mxREAL); if (!b_m) - { - ostringstream tmp; - tmp << " in Simulate_One_Boundary, can't allocate b_m vector\n"; - throw FatalExceptionHandling(tmp.str()); - } + throw FatalExceptionHandling(" in Simulate_One_Boundary, can't allocate b_m vector\n"); A_m = mxCreateSparse(size, size, min(static_cast(IM_i.size()*2), size * size), mxREAL); if (!A_m) - { - ostringstream tmp; - tmp << " in Simulate_One_Boundary, can't allocate A_m matrix\n"; - throw FatalExceptionHandling(tmp.str()); - } + throw FatalExceptionHandling(" in Simulate_One_Boundary, can't allocate A_m matrix\n"); x0_m = mxCreateDoubleMatrix(size, 1, mxREAL); if (!x0_m) - { - ostringstream tmp; - tmp << " in Simulate_One_Boundary, can't allocate x0_m vector\n"; - throw FatalExceptionHandling(tmp.str()); - } - if (!((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 4) && !steady_state))) + throw FatalExceptionHandling(" in Simulate_One_Boundary, can't allocate x0_m vector\n"); + if (!((solve_algo == 6 && steady_state) + || ((stack_solve_algo == 0 || stack_solve_algo == 4) && !steady_state))) { Init_Matlab_Sparse_Simple(size, IM_i, A_m, b_m, zero_solution, x0_m); A_m_save = mxDuplicateArray(A_m); @@ -4160,22 +3932,20 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in Ax_save = static_cast(mxMalloc(Ap[size] * sizeof(double))); test_mxMalloc(Ax_save, __LINE__, __FILE__, __func__, Ap[size] * sizeof(double)); } - memcpy(Ap_save, Ap, (size + 1) * sizeof(SuiteSparse_long)); - memcpy(Ai_save, Ai, Ap[size] * sizeof(SuiteSparse_long)); - memcpy(Ax_save, Ax, Ap[size] * sizeof(double)); - memcpy(b_save, b, size * sizeof(double)); + copy_n(Ap, size + 1, Ap_save); + copy_n(Ai, Ap[size], Ai_save); + copy_n(Ax, Ap[size], Ax_save); + copy_n(b, size, b_save); } } if (zero_solution) - { - for (int i = 0; i < size; i++) - { - int eq = index_vara[i]; - double yy = -(y[eq+it_*y_size]); - direction[eq] = yy; - y[eq+it_*y_size] += slowc * yy; - } - } + for (int i = 0; i < size; i++) + { + int eq = index_vara[i]; + double yy = -(y[eq+it_*y_size]); + direction[eq] = yy; + y[eq+it_*y_size] += slowc * yy; + } else { if ((solve_algo == 5 && steady_state) || (stack_solve_algo == 5 && !steady_state)) @@ -4215,8 +3985,7 @@ dynSparseMatrix::solve_non_linear(int block_num, int y_size, int y_kmin, int y_k bool cvg = false; iter = 0; glambda2 = g0 = very_big; - //try_at_iteration = 0; - while ((!cvg) && (iter < maxit_)) + while (!cvg && iter < maxit_) { cvg = solve_linear(block_num, y_size, y_kmin, y_kmax, size, iter); g0 = res2; @@ -4224,24 +3993,27 @@ dynSparseMatrix::solve_non_linear(int block_num, int y_size, int y_kmin, int y_k } if (!cvg) { - ostringstream tmp; if (steady_state) - tmp << " in Solve Forward complete, convergence not achieved in block " << block_num+1 << ", after " << iter << " iterations\n"; + throw FatalExceptionHandling(" in Solve Forward complete, convergence not achieved in block " + + to_string(block_num+1) + ", after " + to_string(iter) + + " iterations\n"); else - tmp << " in Solve Forward complete, convergence not achieved in block " << block_num+1 << ", at time " << it_ << ", after " << iter << " iterations\n"; - throw FatalExceptionHandling(tmp.str()); + throw FatalExceptionHandling(" in Solve Forward complete, convergence not achieved in block " + + to_string(block_num+1) + ", at time " + to_string(it_) + + ", after " + to_string(iter) + " iterations\n"); } } void -dynSparseMatrix::Simulate_Newton_One_Boundary(const bool forward) +dynSparseMatrix::Simulate_Newton_One_Boundary(bool forward) { g1 = static_cast(mxMalloc(size*size*sizeof(double))); test_mxMalloc(g1, __LINE__, __FILE__, __func__, size*size*sizeof(double)); r = static_cast(mxMalloc(size*sizeof(double))); test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double)); iter = 0; - if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state)) + if ((solve_algo == 6 && steady_state) + || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state)) { Ap_save = static_cast(mxMalloc((size + 1) * sizeof(SuiteSparse_long))); test_mxMalloc(Ap_save, __LINE__, __FILE__, __func__, (size + 1) * sizeof(SuiteSparse_long)); @@ -4264,30 +4036,23 @@ dynSparseMatrix::Simulate_Newton_One_Boundary(const bool forward) else if (forward) { if (!is_linear) - { - for (it_ = y_kmin; it_ < periods+y_kmin; it_++) - solve_non_linear(block_num, y_size, y_kmin, y_kmax, size); - } + for (it_ = y_kmin; it_ < periods+y_kmin; it_++) + solve_non_linear(block_num, y_size, y_kmin, y_kmax, size); else - { - for (int it_ = y_kmin; it_ < periods+y_kmin; it_++) - solve_linear(block_num, y_size, y_kmin, y_kmax, size, 0); - } + for (int it_ = y_kmin; it_ < periods+y_kmin; it_++) + solve_linear(block_num, y_size, y_kmin, y_kmax, size, 0); } else { if (!is_linear) - { - for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--) - solve_non_linear(block_num, y_size, y_kmin, y_kmax, size); - } + for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--) + solve_non_linear(block_num, y_size, y_kmin, y_kmax, size); else - { - for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--) - solve_linear(block_num, y_size, y_kmin, y_kmax, size, 0); - } + for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--) + solve_linear(block_num, y_size, y_kmin, y_kmax, size, 0); } - if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state)) + if ((solve_algo == 6 && steady_state) + || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state)) { mxFree(Ap_save); mxFree(Ai_save); @@ -4329,7 +4094,7 @@ dynSparseMatrix::preconditioner_print_out(string s, int preconditioner, bool ss) } void -dynSparseMatrix::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) +dynSparseMatrix::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, const char *P_endo_names, const vector_table_conditional_local_type &vector_table_conditional_local) { double top = 0.5; double bottom = 0.1; @@ -4375,14 +4140,16 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin else mexPrintf(" variable %s (%d) at time %d = %f direction = %f\n", res.str().c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); } - ostringstream Error; if (iter == 0) - Error << " in Simulate_Newton_Two_Boundaries, the initial values of endogenous variables are too far from the solution.\nChange them!\n"; + throw FatalExceptionHandling(" in Simulate_Newton_Two_Boundaries, the initial values of endogenous variables are too far from the solution.\nChange them!\n"); else - Error << " in Simulate_Newton_Two_Boundaries, dynare cannot improve the simulation in block " << blck+1 << " at time " << it_+1 << " (variable " << index_vara[max_res_idx]+1 << " = " << max_res << ")\n"; - throw FatalExceptionHandling(Error.str()); + throw FatalExceptionHandling(" in Simulate_Newton_Two_Boundaries, dynare cannot improve the simulation in block " + + to_string(blck+1) + " at time " + to_string(it_+1) + + " (variable " + to_string(index_vara[max_res_idx]+1) + + " = " + to_string(max_res) + ")\n"); } - if (!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)) && (stack_solve_algo == 4 || stack_solve_algo == 5)) + if (!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)) + && (stack_solve_algo == 4 || stack_solve_algo == 5)) { if (try_at_iteration == 0) { @@ -4393,10 +4160,15 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin { double t1 = res2 - gp0 * slowc_save - g0; double t2 = glambda2 - gp0 * prev_slowc_save - g0; - double a = (1/(slowc_save * slowc_save) * t1 - 1/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save); - double b = (-prev_slowc_save/(slowc_save * slowc_save) * t1 + slowc_save/(prev_slowc_save * prev_slowc_save) * t2) / (slowc_save - prev_slowc_save); + double a = (1/(slowc_save * slowc_save) * t1 + - 1/(prev_slowc_save * prev_slowc_save) * t2) + / (slowc_save - prev_slowc_save); + double b = (-prev_slowc_save/(slowc_save * slowc_save) * t1 + + slowc_save/(prev_slowc_save * prev_slowc_save) * t2) + / (slowc_save - prev_slowc_save); prev_slowc_save = slowc_save; - slowc_save = max(min(-b + sqrt(b*b - 3 * a * gp0) / (3 * a), top * slowc_save), bottom * slowc_save); + slowc_save = max(min(-b + sqrt(b*b - 3 * a * gp0) / (3 * a), + top * slowc_save), bottom * slowc_save); } glambda2 = res2; try_at_iteration++; @@ -4439,7 +4211,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin markowitz_c = markowitz_c_s; alt_symbolic_count++; } - if (((res1/res1a-1) > -0.3) && symbolic && iter > 0) + if (res1/res1a-1 > -0.3 && symbolic && iter > 0) { if (restart > 2) { @@ -4503,9 +4275,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin mexEvalString("drawnow;"); } if (cvg) - { - return; - } + return; else { if (stack_solve_algo == 5) @@ -4514,27 +4284,15 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin { b_m = mxCreateDoubleMatrix(periods*Size, 1, mxREAL); if (!b_m) - { - ostringstream tmp; - tmp << " in Simulate_Newton_Two_Boundaries, can't allocate b_m vector\n"; - throw FatalExceptionHandling(tmp.str()); - } + throw FatalExceptionHandling(" in Simulate_Newton_Two_Boundaries, can't allocate b_m vector\n"); x0_m = mxCreateDoubleMatrix(periods*Size, 1, mxREAL); if (!x0_m) - { - ostringstream tmp; - tmp << " in Simulate_Newton_Two_Boundaries, can't allocate x0_m vector\n"; - throw FatalExceptionHandling(tmp.str()); - } + throw FatalExceptionHandling(" in Simulate_Newton_Two_Boundaries, can't allocate x0_m vector\n"); if (stack_solve_algo != 0 && stack_solve_algo != 4 && stack_solve_algo != 7) { A_m = mxCreateSparse(periods*Size, periods*Size, IM_i.size()* periods*2, mxREAL); if (!A_m) - { - ostringstream tmp; - tmp << " in Simulate_Newton_Two_Boundaries, can't allocate A_m matrix\n"; - throw FatalExceptionHandling(tmp.str()); - } + throw FatalExceptionHandling(" in Simulate_Newton_Two_Boundaries, can't allocate A_m matrix\n"); } if (stack_solve_algo == 0 || stack_solve_algo == 4) Init_UMFPACK_Sparse(periods, y_kmin, y_kmax, Size, IM_i, &Ap, &Ai, &Ax, &b, x0_m, vector_table_conditional_local, blck); @@ -4559,14 +4317,13 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin mexPrintf("(** %f milliseconds **)\n", 1000.0*(static_cast(t2) - static_cast(t1))/static_cast(CLOCKS_PER_SEC)); mexEvalString("drawnow;"); } - if ((!steady_state && (stack_solve_algo == 4 /*|| stack_solve_algo == 0*/)) /* || steady_state*/) + if (!steady_state && stack_solve_algo == 4) { clock_t t2 = clock(); double ax = -0.1, bx = 1.1, cx = 0.5, fa, fb, fc, xmin; if (!mnbrak(&ax, &bx, &cx, &fa, &fb, &fc)) return; - //mexPrintf("ax= %f, bx=%f, cx=%f, fa=%f, fb=%f, fc=%d\n", ax, bx, cx, fa, fb, fc); if (!golden(ax, bx, cx, 1e-1, solve_tolf, &xmin)) return; slowc = xmin; @@ -4577,7 +4334,6 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin time00 = clock(); if (tbreak_g == 0) tbreak_g = periods; - return; } void @@ -4588,11 +4344,11 @@ dynSparseMatrix::fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_p #ifdef DEBUG mexPrintf("fixe_u : alloc(%d double)\n", u_count_alloc); #endif - (*u) = static_cast(mxMalloc(u_count_alloc*sizeof(double))); + *u = static_cast(mxMalloc(u_count_alloc*sizeof(double))); test_mxMalloc(*u, __LINE__, __FILE__, __func__, u_count_alloc*sizeof(double)); #ifdef DEBUG mexPrintf("*u=%d\n", *u); #endif - memset((*u), 0, u_count_alloc*sizeof(double)); + fill_n(*u, u_count_alloc, 0); u_count_init = max_lag_plus_max_lead_plus_1; } diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index 2721652b6..a25955883 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/mex/sources/bytecode/SparseMatrix.hh @@ -20,16 +20,19 @@ #ifndef SPARSEMATRIX_HH_INCLUDED #define SPARSEMATRIX_HH_INCLUDED -#include -#include +#include +#include #include -#include #include -#include "dynblas.h" +#include +#include +#include +#include + #include "dynumfpack.h" +#include "dynmex.h" #include "Mem_Mngr.hh" -#include "ErrorHandling.hh" #include "Evaluate.hh" using namespace std; @@ -50,7 +53,7 @@ class dynSparseMatrix : public Evaluate public: 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); - 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_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, const char *P_endo_names, const 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(const string &file_name, int Size, int periods, int y_kmin, int y_kmax, bool two_boundaries, int stack_solve_algo, int solve_algo); @@ -93,7 +96,6 @@ private: 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 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(int r, int c, int u_index, int lag_index); void Delete(int r, int c); int At_Row(int r, NonZeroElem **first) const; diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 53764e59c..e9296004f 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -16,13 +16,13 @@ * You should have received a copy of the GNU General Public License * along with Dynare. If not, see . */ -#include -#include "Interpreter.hh" -#include "ErrorHandling.hh" + #include #include +#include -void (*prev_fn)(int); +#include "Interpreter.hh" +#include "ErrorHandling.hh" string Get_Argument(const mxArray *prhs) @@ -39,9 +39,6 @@ Get_Argument(const mxArray *prhs) return f; } -//#include -#include - string deblank(string x) { @@ -254,82 +251,42 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (extended_path) { - if (extended_path_struct == nullptr) - { - string tmp = "The 'extended_path' option must be followed by the extended_path descriptor"; - mexErrMsgTxt(tmp.c_str()); - } + if (!extended_path_struct) + mexErrMsgTxt("The 'extended_path' option must be followed by the extended_path descriptor"); mxArray *date_str = mxGetField(extended_path_struct, 0, "date_str"); - if (date_str == nullptr) - { - string tmp = "date_str"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!date_str) + mexErrMsgTxt("The extended_path description structure does not contain the member: date_str"); int nb_periods = mxGetM(date_str) * mxGetN(date_str); mxArray *constrained_vars_ = mxGetField(extended_path_struct, 0, "constrained_vars_"); - if (constrained_vars_ == nullptr) - { - string tmp = "constrained_vars_"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!constrained_vars_) + mexErrMsgTxt("The extended_path description structure does not contain the member: constrained_vars_"); mxArray *constrained_paths_ = mxGetField(extended_path_struct, 0, "constrained_paths_"); - if (constrained_paths_ == nullptr) - { - string tmp = "constrained_paths_"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!constrained_paths_) + mexErrMsgTxt("The extended_path description structure does not contain the member: constrained_paths_"); mxArray *constrained_int_date_ = mxGetField(extended_path_struct, 0, "constrained_int_date_"); - if (constrained_int_date_ == nullptr) - { - string tmp = "constrained_int_date_"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!constrained_int_date_) + mexErrMsgTxt("The extended_path description structure does not contain the member: constrained_int_date_"); mxArray *constrained_perfect_foresight_ = mxGetField(extended_path_struct, 0, "constrained_perfect_foresight_"); - if (constrained_perfect_foresight_ == nullptr) - { - string tmp = "constrained_perfect_foresight_"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } - + if (!constrained_perfect_foresight_) + mexErrMsgTxt("The extended_path description structure does not contain the member: constrained_perfect_foresight_"); mxArray *shock_var_ = mxGetField(extended_path_struct, 0, "shock_vars_"); - if (shock_var_ == nullptr) - { - string tmp = "shock_vars_"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!shock_var_) + mexErrMsgTxt("The extended_path description structure does not contain the member: shock_vars_"); mxArray *shock_paths_ = mxGetField(extended_path_struct, 0, "shock_paths_"); - if (shock_paths_ == nullptr) - { - string tmp = "shock_paths_"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!shock_paths_) + mexErrMsgTxt("The extended_path description structure does not contain the member: shock_paths_"); mxArray *shock_int_date_ = mxGetField(extended_path_struct, 0, "shock_int_date_"); - if (shock_int_date_ == nullptr) - { - string tmp = "shock_int_date_"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!shock_int_date_) + mexErrMsgTxt("The extended_path description structure does not contain the member: shock_int_date_"); mxArray *shock_str_date_ = mxGetField(extended_path_struct, 0, "shock_str_date_"); - if (shock_str_date_ == nullptr) - { - string tmp = "shock_str_date_"; - tmp.insert(0, "The extended_path description structure does not contain the member: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!shock_str_date_) + mexErrMsgTxt("The extended_path description structure does not contain the member: shock_str_date_"); int nb_constrained = mxGetM(constrained_vars_) * mxGetN(constrained_vars_); int nb_controlled = 0; mxArray *options_cond_fcst_ = mxGetField(extended_path_struct, 0, "options_cond_fcst_"); mxArray *controlled_varexo = nullptr; - if (options_cond_fcst_ != nullptr) + if (options_cond_fcst_) { controlled_varexo = mxGetField(options_cond_fcst_, 0, "controlled_varexo"); nb_controlled = mxGetM(controlled_varexo) * mxGetN(controlled_varexo); @@ -337,7 +294,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) 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) + if (controlled_varexo) controlled_varexo_value = mxGetPr(controlled_varexo); double *constrained_var_value = mxGetPr(constrained_vars_); sconditional_extended_path.resize(nb_constrained); @@ -372,7 +329,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) int *constrained_int_date = static_cast(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) - mexErrMsgTxt((string{"The total number of simulation periods ("} + to_string(nb_periods) + mexErrMsgTxt(("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()); @@ -426,10 +383,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) char *buf = static_cast(mxCalloc(buflen, sizeof(char))); int info = mxGetString(mxGetCell(date_str, i), buf, buflen); if (info) - { - string tmp = "Can not allocated memory to store the date_str in the extended path descriptor"; - mexErrMsgTxt(tmp.c_str()); - } + mexErrMsgTxt("Can not allocated memory to store the date_str in the extended path descriptor"); dates.emplace_back(buf); //string(Dates[i]); mxFree(buf); } @@ -437,12 +391,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (plan.length() > 0) { mxArray *plan_struct = mexGetVariable("base", plan.c_str()); - if (plan_struct == nullptr) - { - string tmp = plan; - tmp.insert(0, "Can't find the plan: "); - mexErrMsgTxt(tmp.c_str()); - } + if (!plan_struct) + mexErrMsgTxt(("Can't find the plan: " + plan).c_str()); size_t n_plan = mxGetN(plan_struct); splan.resize(n_plan); for (int i = 0; i < static_cast(n_plan); i++) @@ -460,12 +410,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (variable_type == SymbolType::exogenous || variable_type == SymbolType::exogenousDet) splan[i].var_num = exo_num; else - { - string tmp = name; - tmp.insert(0, "the variable '"); - tmp.append("' defined as var in plan is not an exogenous or a deterministic exogenous\n"); - mexErrMsgTxt(tmp.c_str()); - } + mexErrMsgTxt(("The variable '" + string{name} + "' defined as var in plan is not an exogenous or a deterministic exogenous\n").c_str()); } tmp = mxGetField(plan_struct, i, "var"); if (tmp) @@ -478,12 +423,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (variable_type == SymbolType::endogenous) splan[i].exo_num = exo_num; else - { - string tmp = name; - tmp.insert(0, "the variable '"); - tmp.append("' defined as exo in plan is not an endogenous variable\n"); - mexErrMsgTxt(tmp.c_str()); - } + mexErrMsgTxt(("The variable '" + string{name} + "' defined as exo in plan is not an endogenous variable\n").c_str()); } tmp = mxGetField(plan_struct, i, "per_value"); if (tmp) @@ -532,12 +472,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (variable_type == SymbolType::exogenous || variable_type == SymbolType::exogenousDet) splan[i].var_num = exo_num; else - { - string tmp = name; - tmp.insert(0, "the variable '"); - tmp.append("' defined as var in pfplan is not an exogenous or a deterministic exogenous\n"); - mexErrMsgTxt(tmp.c_str()); - } + mexErrMsgTxt(("The variable '" + string{name} + "' defined as var in pfplan is not an exogenous or a deterministic exogenous\n").c_str()); } tmp = mxGetField(pfplan_struct, i, "exo"); if (tmp) @@ -550,12 +485,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (variable_type == SymbolType::endogenous) spfplan[i].exo_num = exo_num; else - { - string tmp = name; - tmp.insert(0, "the variable '"); - tmp.append("' defined as exo in pfplan is not an endogenous variable\n"); - mexErrMsgTxt(tmp.c_str()); - } + mexErrMsgTxt(("The variable '" + string{name} + "' defined as exo in pfplan is not an endogenous variable\n").c_str()); } tmp = mxGetField(pfplan_struct, i, "per_value"); if (tmp) @@ -828,11 +758,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) } else { - int out_periods; - if (extended_path) - out_periods = max_periods + y_kmin; - else - out_periods = row_y; + int out_periods = extended_path ? max_periods + y_kmin : row_y; plhs[0] = mxCreateDoubleMatrix(out_periods, static_cast(col_y), mxREAL); pind = mxGetPr(plhs[0]); for (i = 0; i < out_periods*col_y; i++) @@ -841,11 +767,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) } else { - int out_periods; - if (extended_path) - out_periods = max_periods + y_kmin; - else - out_periods = col_y; + int out_periods = extended_path ? max_periods + y_kmin : col_y; plhs[0] = mxCreateDoubleMatrix(static_cast(row_y), out_periods, mxREAL); pind = mxGetPr(plhs[0]); if (evaluate) @@ -862,7 +784,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (evaluate) { - int jacob_field_number = 0, jacob_exo_field_number = 0, jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0; + int jacob_field_number = 0, jacob_exo_field_number = 0, + jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0; if (!block_structur) { const char *field_names[] = {"g1", "g1_x", "g1_xd", "g1_o"}; @@ -895,18 +818,16 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mexErrMsgTxt("Fatal error in bytecode: in main, cannot add extra field jacob_other_endo to the structArray\n"); } if (!dont_store_a_structure) - { - for (int i = 0; i < nb_blocks; i++) - { - mxSetFieldByNumber(plhs[1], i, jacob_field_number, interprete.get_jacob(i)); - if (!steady_state) - { - mxSetFieldByNumber(plhs[1], i, jacob_exo_field_number, interprete.get_jacob_exo(i)); - mxSetFieldByNumber(plhs[1], i, jacob_exo_det_field_number, interprete.get_jacob_exo_det(i)); - mxSetFieldByNumber(plhs[1], i, jacob_other_endo_field_number, interprete.get_jacob_other_endo(i)); - } - } - } + for (int i = 0; i < nb_blocks; i++) + { + mxSetFieldByNumber(plhs[1], i, jacob_field_number, interprete.get_jacob(i)); + if (!steady_state) + { + mxSetFieldByNumber(plhs[1], i, jacob_exo_field_number, interprete.get_jacob_exo(i)); + mxSetFieldByNumber(plhs[1], i, jacob_exo_det_field_number, interprete.get_jacob_exo_det(i)); + mxSetFieldByNumber(plhs[1], i, jacob_other_endo_field_number, interprete.get_jacob_other_endo(i)); + } + } } else {