diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh index 760d0ec7c..39fc01ec3 100644 --- a/mex/sources/bytecode/ErrorHandling.hh +++ b/mex/sources/bytecode/ErrorHandling.hh @@ -24,8 +24,20 @@ #include #include #include +#include #define BYTE_CODE #include "CodeInterpreter.hh" + +#define _USE_MATH_DEFINES +#include +#ifndef M_PI +#define M_PI (3.14159265358979323846) +#endif + +#ifndef M_SQRT2 +#define M_SQRT2 1.41421356237309504880 +#endif + #ifdef DEBUG_EX # include # include "mex_interface.hh" @@ -37,15 +49,6 @@ # define CHAR_LENGTH 2 #endif -//Work around for: https://sourceware.org/bugzilla/show_bug.cgi?id=19439 -#ifndef __builtin_isnan -# define isnan(x) std::isnan(x) -#endif - -#ifndef __builtin_isinf -# define isinf(x) std::isinf(x) -#endif - #ifdef _MSC_VER #include #define M_E 2.71828182845904523536 @@ -250,7 +253,7 @@ public: value2(value2_arg) { ostringstream tmp; - if (abs(value1) > 1e-10 ) + if (fabs(value1) > 1e-10 ) tmp << " with X=" << value1 << "\n"; else tmp << " with X=" << value1 << " and a=" << value2 << "\n"; @@ -2231,6 +2234,17 @@ public: it_code_ret = it_code; return (tmp_out.str()); } + void + +inline test_mxMalloc(void* z, int line, string file, string func, int amount) +{ + if (!z && (amount > 0)) + { + ostringstream tmp; + tmp << " mxMalloc: out of memory " << amount << " bytes required at line " << line << " in function " << func << " (file " << file; + throw FatalExceptionHandling(tmp.str()); + } +} }; diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc index b0208bc0e..0c7e242a4 100644 --- a/mex/sources/bytecode/Evaluate.cc +++ b/mex/sources/bytecode/Evaluate.cc @@ -140,7 +140,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int if ( utIsInterruptPending() ) throw UserExceptionHandling(); #endif - + while (go_on) { #ifdef DEBUG @@ -1172,6 +1172,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int case ExternalFunctionWithFirstandSecondDerivative: { input_arguments = (mxArray **) mxMalloc(nb_input_arguments * sizeof(mxArray *)); + test_mxMalloc(input_arguments, __LINE__, __FILE__, __func__, nb_input_arguments * sizeof(mxArray *)); #ifdef DEBUG mexPrintf("Stack.size()=%d\n", Stack.size()); mexEvalString("drawnow;"); @@ -1188,7 +1189,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int tmp << " external function: " << function_name << " not found"; throw FatalExceptionHandling(tmp.str()); } - + double *rr = mxGetPr(output_arguments[0]); Stack.push(*rr); if (function_type == ExternalFunctionWithFirstDerivative || function_type == ExternalFunctionWithFirstandSecondDerivative) @@ -1215,6 +1216,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int case ExternalFunctionNumericalFirstDerivative: { input_arguments = (mxArray **) mxMalloc((nb_input_arguments+1+nb_add_input_arguments) * sizeof(mxArray *)); + test_mxMalloc(input_arguments, __LINE__, __FILE__, __func__, (nb_input_arguments+1+nb_add_input_arguments) * sizeof(mxArray *)); mxArray *vv = mxCreateString(arg_func_name.c_str()); input_arguments[0] = vv; vv = mxCreateDoubleScalar(fc->get_row()); @@ -1254,6 +1256,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int case ExternalFunctionFirstDerivative: { input_arguments = (mxArray **) mxMalloc(nb_input_arguments * sizeof(mxArray *)); + test_mxMalloc(input_arguments, __LINE__, __FILE__, __func__, nb_input_arguments * sizeof(mxArray *)); for (unsigned int i = 0; i < nb_input_arguments; i++) { mxArray *vv = mxCreateDoubleScalar(Stack.top()); @@ -1277,6 +1280,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int case ExternalFunctionNumericalSecondDerivative: { input_arguments = (mxArray **) mxMalloc((nb_input_arguments+1+nb_add_input_arguments) * sizeof(mxArray *)); + test_mxMalloc(input_arguments, __LINE__, __FILE__, __func__, (nb_input_arguments+1+nb_add_input_arguments) * sizeof(mxArray *)); mxArray *vv = mxCreateString(arg_func_name.c_str()); input_arguments[0] = vv; vv = mxCreateDoubleScalar(fc->get_row()); @@ -1315,6 +1319,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int case ExternalFunctionSecondDerivative: { input_arguments = (mxArray **) mxMalloc(nb_input_arguments * sizeof(mxArray *)); + test_mxMalloc(input_arguments, __LINE__, __FILE__, __func__, nb_input_arguments * sizeof(mxArray *)); for (unsigned int i = 0; i < nb_input_arguments; i++) { mxArray *vv = mxCreateDoubleScalar(Stack.top()); @@ -1583,7 +1588,9 @@ void Evaluate::solve_simple_over_periods(const bool forward) { g1 = (double *) mxMalloc(sizeof(double)); + test_mxMalloc(g1, __LINE__, __FILE__, __func__, sizeof(double)); r = (double *) mxMalloc(sizeof(double)); + test_mxMalloc(r, __LINE__, __FILE__, __func__, sizeof(double)); start_code = it_code; if (steady_state) { diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh index 33621b7a4..e5fa778e3 100644 --- a/mex/sources/bytecode/Evaluate.hh +++ b/mex/sources/bytecode/Evaluate.hh @@ -63,13 +63,13 @@ protected: double solve_tolf; bool GaussSeidel; map, int>, int> IM_i; - int equation, derivative_equation, derivative_variable; + int equation, derivative_equation, derivative_variable; string filename; int stack_solve_algo, solve_algo; bool global_temporary_terms; bool print, print_error; double res1, res2, max_res; - int max_res_idx; + int max_res_idx; vector Block_Contain; int size; diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index f85c9e709..48c600190 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -77,7 +77,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub GlobalTemporaryTerms = GlobalTemporaryTerms_arg; print_error = print_error_arg; //steady_state = steady_state_arg; - print_it = print_it_arg; + print_it = print_it_arg; } @@ -117,8 +117,10 @@ Interpreter::evaluate_a_block(bool initialization) } break; case SOLVE_FORWARD_SIMPLE: - g1 = (double *) mxMalloc(size*size*sizeof(double)); - r = (double *) mxMalloc(size*sizeof(double)); + g1 = (double *) mxMalloc(size*size*sizeof(double)); + test_mxMalloc(g1, __LINE__, __FILE__, __func__, size*size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double)); if (steady_state) { compute_block_time(0, true, /*block_num, size, steady_state,*/ false); @@ -144,8 +146,8 @@ Interpreter::evaluate_a_block(bool initialization) for (int j = 0; j < size; j++) residual[it_*size+j] = r[j]; } - } - mxFree(g1); + } + mxFree(g1); mxFree(r); break; case SOLVE_FORWARD_COMPLETE: @@ -157,7 +159,8 @@ Interpreter::evaluate_a_block(bool initialization) #ifdef DEBUG mexPrintf("in SOLVE_FORWARD_COMPLETE r = mxMalloc(%d*sizeof(double))\n", size); #endif - r = (double *) mxMalloc(size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double)); if (steady_state) { compute_block_time(0, true, /*block_num, size, steady_state,*/ false); @@ -215,8 +218,10 @@ Interpreter::evaluate_a_block(bool initialization) } break; case SOLVE_BACKWARD_SIMPLE: - g1 = (double *) mxMalloc(size*size*sizeof(double)); - r = (double *) mxMalloc(size*sizeof(double)); + g1 = (double *) mxMalloc(size*size*sizeof(double)); + test_mxMalloc(g1, __LINE__, __FILE__, __func__, size*size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double)); if (steady_state) { compute_block_time(0, true, /*block_num, size, steady_state,*/ false); @@ -252,7 +257,8 @@ Interpreter::evaluate_a_block(bool initialization) fixe_u(&u, u_count_int, u_count_int); Read_SparseMatrix(bin_base_name, size, 1, 0, 0, false, stack_solve_algo, solve_algo); } - r = (double *) mxMalloc(size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double)); if (steady_state) { compute_block_time(0, true, /*block_num, size, steady_state,*/ false); @@ -289,7 +295,8 @@ Interpreter::evaluate_a_block(bool initialization) Read_SparseMatrix(bin_base_name, size, periods, y_kmin, y_kmax, true, stack_solve_algo, solve_algo); } u_count = u_count_int*(periods+y_kmax+y_kmin); - r = (double *) mxMalloc(size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double)); begining = it_code; for (it_ = y_kmin; it_ < periods+y_kmin; it_++) { @@ -430,9 +437,12 @@ Interpreter::simulate_a_block(vector_table_conditional_local_type vector_table_c Read_SparseMatrix(bin_base_name, size, periods, y_kmin, y_kmax, true, stack_solve_algo, solve_algo); } u_count = u_count_int*(periods+y_kmax+y_kmin); - r = (double *) mxMalloc(size*sizeof(double)); - res = (double *) mxMalloc(size*periods*sizeof(double)); - y_save = (double *) mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); + r = (double *) mxMalloc(size*sizeof(double)); + test_mxMalloc(r, __LINE__, __FILE__, __func__, size*sizeof(double)); + res = (double *) mxMalloc(size*periods*sizeof(double)); + test_mxMalloc(res, __LINE__, __FILE__, __func__, size*periods*sizeof(double)); + y_save = (double *) mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); + test_mxMalloc(y_save, __LINE__, __FILE__, __func__, y_size*sizeof(double)*(periods+y_kmax+y_kmin)); start_code = it_code; iter = 0; if (!is_linear) @@ -457,7 +467,7 @@ Interpreter::simulate_a_block(vector_table_conditional_local_type vector_table_c //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; } - + } } compute_complete_2b(false, &res1, &res2, &max_res, &max_res_idx); @@ -498,13 +508,19 @@ Interpreter::simulate_a_block(vector_table_conditional_local_type vector_table_c 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); max_res = 0; max_res_idx = 0; } - it_code = end_code; - mxFree(r); - mxFree(y_save); - mxFree(u); - mxFree(index_vara); - mxFree(index_equa); - mxFree(res); + it_code = end_code; + if (r) + mxFree(r); + if (y_save) + mxFree(y_save); + if (u) + mxFree(u); + if (index_vara) + mxFree(index_vara); + if (index_equa) + mxFree(index_equa); + if (res) + mxFree(res); memset(direction, 0, size_of_direction); End_Solver(); break; @@ -558,14 +574,14 @@ Interpreter::print_a_block() } } -void +void Interpreter::ReadCodeFile(string file_name, CodeLoad &code) { if (steady_state) file_name += "_static"; else file_name += "_dynamic"; - + //First read and store in memory the code code_liste = code.get_op_code(file_name); EQN_block_number = code.get_block_number(); @@ -653,9 +669,7 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo Block_Contain = fb->get_Block_Contain(); it_code++; if (constrained) - { - check_for_controlled_exo_validity(fb,sconstrained_extended_path); - } + check_for_controlled_exo_validity(fb,sconstrained_extended_path); set_block(fb->get_size(), fb->get_type(), file_name, bin_basename, Block_Count, fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int(), block); if (print) print_a_block(); @@ -690,7 +704,7 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo #ifdef DEBUG mexPrintf("endo in Block_Count=%d, block=%d, type=%d, steady_state=%d, print_it=%d, Block_Count=%d, fb->get_is_linear()=%d, fb->get_endo_nbr()=%d, fb->get_Max_Lag()=%d, fb->get_Max_Lead()=%d, fb->get_u_count_int()=%d\n", Block_Count, fb->get_size(), fb->get_type(), steady_state, print_it, Block_Count, fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); -#endif +#endif bool result; if (sconstrained_extended_path.size()) { @@ -702,7 +716,7 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo jacobian_other_endo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_other_endo_jacob(), mxREAL)); residual = vector(fb->get_size()*(periods+y_kmin)); result = simulate_a_block(vector_table_conditional_local); - + mxDestroyArray(jacobian_block.back()); jacobian_block.pop_back(); mxDestroyArray(jacobian_exo_block.back()); @@ -713,10 +727,7 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo jacobian_other_endo_block.pop_back(); } else - { - result = simulate_a_block(vector_table_conditional_local); - } - + result = simulate_a_block(vector_table_conditional_local); //mexPrintf("OKe\n"); if (max_res > max_res_local) { @@ -749,8 +760,9 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo #endif var = ((FDIMT_ *) it_code->second)->get_size(); if (T) - mxFree(T); - T = (double *) mxMalloc(var*(periods+y_kmin+y_kmax)*sizeof(double)); + mxFree(T); + T = (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); @@ -775,8 +787,12 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo GlobalTemporaryTerms = mxCreateDoubleMatrix(var, 1, mxREAL); T = mxGetPr(GlobalTemporaryTerms); } - else - T = (double *) mxMalloc(var*sizeof(double)); + else + { + T = (double *) mxMalloc(var*sizeof(double)); + test_mxMalloc(T, __LINE__, __FILE__, __func__, var*sizeof(double)); + } + if (block >= 0) it_code = code_liste.begin() + code.get_begin_block(block); @@ -842,9 +858,18 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate, 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); + /*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)); + double *x_save = (double *) mxMalloc((periods + y_kmax + y_kmin) * col_x *sizeof(double));*/ + int max_periods = max(periods, nb_periods); + size_t size_of_direction = y_size*(max_periods + y_kmax + y_kmin)*sizeof(double); + double *y_save = (double *) mxMalloc(size_of_direction); + test_mxMalloc(y_save, __LINE__, __FILE__, __func__, size_of_direction); + + + double *x_save = (double *) mxMalloc((max_periods + y_kmax + y_kmin) * col_x *sizeof(double)); + test_mxMalloc(x_save, __LINE__, __FILE__, __func__, (max_periods + y_kmax + y_kmin) * col_x *sizeof(double)); + vector_table_conditional_local_type vector_table_conditional_local; vector_table_conditional_local.clear(); @@ -896,7 +921,7 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate, vector_table_conditional_local = table_conditional_global[t]; if (t < nb_periods) MainLoop(bin_basename, code, evaluate, block, false, true, sconstrained_extended_path, vector_table_conditional_local); - else + else MainLoop(bin_basename, code, evaluate, block, true, true, sconstrained_extended_path, vector_table_conditional_local); for (int j = 0; j < y_size; j++) { @@ -909,7 +934,7 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate, x_save[t + y_kmin + j * nb_row_x] = x[y_kmin + j * nb_row_x]; x[y_kmin + j * nb_row_x] = x_save[t + 1 + y_kmin + j * nb_row_x]; } - + if (old_print_it) { ostringstream res, res1; @@ -921,21 +946,23 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate, mexPrintf(line.c_str()); mexEvalString("drawnow;"); } - } + } 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*(periods + y_kmax + y_kmin)); i++) - y[i] = y_save[i]; + y[i] = y_save[i]; for (int j = 0; j < col_x* nb_row_x; j++) x[j] = x_save[j]; - - mxFree(Init_Code->second); - mxFree(y_save); - mxFree(x_save); + if (Init_Code->second) + mxFree(Init_Code->second); + if (y_save) + mxFree(y_save); + if (x_save) + mxFree(x_save); nb_blocks = Block_Count+1; if (T && !global_temporary_terms) mxFree(T); @@ -947,7 +974,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool evaluate { CodeLoad code; ReadCodeFile(file_name, code); - + //The big loop on intructions it_code = code_liste.begin(); it_code_type Init_Code = it_code; @@ -955,8 +982,8 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool evaluate vector_table_conditional_local_type vector_table_conditional_local_junk; MainLoop(bin_basename, code, evaluate, block, true, false, s_plan_junk, vector_table_conditional_local_junk); - - + + mxFree(Init_Code->second); nb_blocks = Block_Count+1; if (T && !global_temporary_terms) diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh index a3d9973ef..8e8a00e0a 100644 --- a/mex/sources/bytecode/Interpreter.hh +++ b/mex/sources/bytecode/Interpreter.hh @@ -27,7 +27,7 @@ #define BYTE_CODE #include "CodeInterpreter.hh" #include "SparseMatrix.hh" -#include "Evaluate.hh" +#include "Evaluate.hh" #ifdef LINBCG # include "linbcg.hh" #endif diff --git a/mex/sources/bytecode/Mem_Mngr.cc b/mex/sources/bytecode/Mem_Mngr.cc index b239b60df..d9d20f7a8 100644 --- a/mex/sources/bytecode/Mem_Mngr.cc +++ b/mex/sources/bytecode/Mem_Mngr.cc @@ -24,7 +24,7 @@ Mem_Mngr::Mem_Mngr() swp_f = false; swp_f_b = 0; } -void +/*void Mem_Mngr::Print_heap() { unsigned int i; @@ -32,7 +32,8 @@ Mem_Mngr::Print_heap() for (i = 0; i < CHUNK_SIZE; i++) mexPrintf("%3d ", i); mexPrintf("\n"); -} +} +*/ void Mem_Mngr::init_Mem() @@ -49,7 +50,7 @@ Mem_Mngr::init_Mem() void Mem_Mngr::fixe_file_name(string filename_arg) { - filename = filename_arg; + filename_mem = filename_arg; } void @@ -77,26 +78,26 @@ Mem_Mngr::mxMalloc_NZE() { CHUNK_SIZE += CHUNK_BLCK_SIZE; Nb_CHUNK++; - NZE_Mem = (NonZeroElem *) mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); /*The block of memory allocated*/ + NZE_Mem = (NonZeroElem *) mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); /*The block of memory allocated*/ + error_msg.test_mxMalloc(NZE_Mem, __LINE__, __FILE__, __func__, CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); NZE_Mem_Allocated.push_back(NZE_Mem); if (!NZE_Mem) - { - mexPrintf("Not enough memory available\n"); - mexEvalString("drawnow;"); + mexPrintf("Not enough memory available\n"); + if (NZE_Mem_add) + { + NZE_Mem_add = (NonZeroElem **) mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to redefine the size of pointer on the memory*/ + error_msg.test_mxMalloc(NZE_Mem_add , __LINE__, __FILE__, __func__, CHUNK_SIZE*sizeof(NonZeroElem *)); } - if (NZE_Mem_add) - NZE_Mem_add = (NonZeroElem **) mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to redefine the size of pointer on the memory*/ - else - NZE_Mem_add = (NonZeroElem **) mxMalloc(CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to define the size of pointer on the memory*/ + else + { + NZE_Mem_add = (NonZeroElem **) mxMalloc(CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to define the size of pointer on the memory*/ + error_msg.test_mxMalloc(NZE_Mem_add , __LINE__, __FILE__, __func__, CHUNK_SIZE*sizeof(NonZeroElem *)); + } + if (!NZE_Mem_add) - { - mexPrintf("Not enough memory available\n"); - mexEvalString("drawnow;"); - } + mexPrintf("Not enough memory available\n"); for (i = CHUNK_heap_pos; i < CHUNK_SIZE; i++) - { - NZE_Mem_add[i] = (NonZeroElem *) (NZE_Mem+(i-CHUNK_heap_pos)); - } + NZE_Mem_add[i] = (NonZeroElem *) (NZE_Mem+(i-CHUNK_heap_pos)); i = CHUNK_heap_pos++; return (NZE_Mem_add[i]); } diff --git a/mex/sources/bytecode/Mem_Mngr.hh b/mex/sources/bytecode/Mem_Mngr.hh index be3c3c404..f57a4d312 100644 --- a/mex/sources/bytecode/Mem_Mngr.hh +++ b/mex/sources/bytecode/Mem_Mngr.hh @@ -19,15 +19,16 @@ #ifndef MEM_MNGR_HH_INCLUDED #define MEM_MNGR_HH_INCLUDED - + +#include "ErrorHandling.hh" #include #include #ifndef DEBUG_EX # include #else # include "mex_interface.hh" -#endif -using namespace std; +#endif +//using namespace std; struct NonZeroElem { @@ -41,7 +42,7 @@ typedef vector v_NonZeroElem; class Mem_Mngr { public: - void Print_heap(); + //void Print_heap(); void init_Mem(); void mxFree_NZE(void *pos); NonZeroElem *mxMalloc_NZE(); @@ -49,7 +50,8 @@ public: void Free_All(); Mem_Mngr(); void fixe_file_name(string filename_arg); - bool swp_f; + bool swp_f; + ErrorMsg error_msg; private: v_NonZeroElem Chunk_Stack; unsigned int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK; @@ -59,7 +61,7 @@ private: vector NZE_Mem_Allocated; int swp_f_b; fstream SaveCode_swp; - string filename; + string filename_mem; }; #endif diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 840f302a6..8cd42ca09 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -183,7 +183,7 @@ dynSparseMatrix::dynSparseMatrix() } 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 + const int minimal_solving_periods_arg, const double slowc_arg #ifdef CUDA , const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg #endif @@ -206,7 +206,7 @@ dynSparseMatrix::dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, con IM_i.clear(); lu_inc_tol = 1e-10; Symbolic = NULL; - Numeric = NULL; + Numeric = NULL; #ifdef CUDA CUDA_device = CUDA_device_arg; cublas_handle = cublas_handle_arg; @@ -621,7 +621,8 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods } } } - index_vara = (int *) mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int)); + index_vara = (int *) mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int)); + test_mxMalloc(index_vara, __LINE__, __FILE__, __func__, Size*(periods+y_kmin+y_kmax)*sizeof(int)); for (int j = 0; j < Size; j++) SaveCode.read(reinterpret_cast(&index_vara[j]), sizeof(*index_vara)); if (periods+y_kmin+y_kmax > 1) @@ -630,7 +631,8 @@ dynSparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods for (int j = 0; j < Size; j++) index_vara[j+Size*i] = index_vara[j+Size*(i-1)] + y_size; } - index_equa = (int *) mxMalloc(Size*sizeof(int)); + index_equa = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(index_equa, __LINE__, __FILE__, __func__, Size*sizeof(int)); for (int j = 0; j < Size; j++) SaveCode.read(reinterpret_cast(&index_equa[j]), sizeof(*index_equa)); } @@ -641,25 +643,38 @@ dynSparseMatrix::Simple_Init(int Size, map, int>, int> &IM, int i, eq, var, lag; map, int>, int>::iterator it4; NonZeroElem *first; - pivot = (int *) mxMalloc(Size*sizeof(int)); - pivot_save = (int *) mxMalloc(Size*sizeof(int)); - pivotk = (int *) mxMalloc(Size*sizeof(int)); - pivotv = (double *) mxMalloc(Size*sizeof(double)); - pivotva = (double *) mxMalloc(Size*sizeof(double)); - b = (int *) mxMalloc(Size*sizeof(int)); + pivot = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(pivot, __LINE__, __FILE__, __func__, Size*sizeof(int)); + pivot_save = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(pivot_save, __LINE__, __FILE__, __func__, Size*sizeof(int)); + pivotk = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(pivotk, __LINE__, __FILE__, __func__, Size*sizeof(int)); + pivotv = (double *) mxMalloc(Size*sizeof(double)); + test_mxMalloc(pivotv, __LINE__, __FILE__, __func__, Size*sizeof(double)); + pivotva = (double *) mxMalloc(Size*sizeof(double)); + test_mxMalloc(pivotva, __LINE__, __FILE__, __func__, Size*sizeof(double)); + b = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(b, __LINE__, __FILE__, __func__, Size*sizeof(int)); line_done = (bool *) mxMalloc(Size*sizeof(bool)); + test_mxMalloc(line_done, __LINE__, __FILE__, __func__, Size*sizeof(bool)); mem_mngr.init_CHUNK_BLCK_SIZE(u_count); g_save_op = NULL; g_nop_all = 0; i = Size*sizeof(NonZeroElem *); - FNZE_R = (NonZeroElem **) mxMalloc(i); - FNZE_C = (NonZeroElem **) mxMalloc(i); - NonZeroElem **temp_NZE_R = (NonZeroElem **) mxMalloc(i); - NonZeroElem **temp_NZE_C = (NonZeroElem **) mxMalloc(i); + FNZE_R = (NonZeroElem **) mxMalloc(i); + test_mxMalloc(FNZE_R, __LINE__, __FILE__, __func__, i); + FNZE_C = (NonZeroElem **) mxMalloc(i); + test_mxMalloc(FNZE_C, __LINE__, __FILE__, __func__, i); + NonZeroElem **temp_NZE_R = (NonZeroElem **) mxMalloc(i); + test_mxMalloc(*temp_NZE_R, __LINE__, __FILE__, __func__, i); + NonZeroElem **temp_NZE_C = (NonZeroElem **) mxMalloc(i); + test_mxMalloc(*temp_NZE_C, __LINE__, __FILE__, __func__, i); i = Size*sizeof(int); - NbNZRow = (int *) mxMalloc(i); - NbNZCol = (int *) mxMalloc(i); + NbNZRow = (int *) mxMalloc(i); + test_mxMalloc(NbNZRow, __LINE__, __FILE__, __func__, i); + NbNZCol = (int *) mxMalloc(i); + test_mxMalloc(NbNZCol, __LINE__, __FILE__, __func__, i); it4 = IM.begin(); eq = -1; for (i = 0; i < Size; i++) @@ -838,7 +853,8 @@ void dynSparseMatrix::Init_UMFPACK_Sparse_Simple(int Size, map, int>, int> &IM, SuiteSparse_long **Ap, SuiteSparse_long **Ai, double **Ax, double **b, bool &zero_solution, mxArray *x0_m) { int eq, var; - *b = (double*)mxMalloc(Size * sizeof(double)); + *b = (double*)mxMalloc(Size * sizeof(double)); + test_mxMalloc(*b, __LINE__, __FILE__, __func__, Size * sizeof(double)); if (!(*b)) { ostringstream tmp; @@ -852,7 +868,8 @@ dynSparseMatrix::Init_UMFPACK_Sparse_Simple(int Size, map, i tmp << " in Init_UMFPACK_Sparse_Simple, can't retrieve x0 vector\n"; throw FatalExceptionHandling(tmp.str()); } - *Ap = (SuiteSparse_long*)mxMalloc((Size+1) * sizeof(SuiteSparse_long)); + *Ap = (SuiteSparse_long*)mxMalloc((Size+1) * sizeof(SuiteSparse_long)); + test_mxMalloc(*Ap, __LINE__, __FILE__, __func__, (Size+1) * sizeof(SuiteSparse_long)); if (!(*Ap)) { ostringstream tmp; @@ -860,14 +877,16 @@ dynSparseMatrix::Init_UMFPACK_Sparse_Simple(int Size, map, i throw FatalExceptionHandling(tmp.str()); } size_t prior_nz = IM.size(); - *Ai = (SuiteSparse_long*)mxMalloc(prior_nz * sizeof(SuiteSparse_long)); + *Ai = (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()); } - *Ax = (double*)mxMalloc(prior_nz * sizeof(double)); + *Ax = (double*)mxMalloc(prior_nz * sizeof(double)); + test_mxMalloc(*Ax, __LINE__, __FILE__, __func__, prior_nz * sizeof(double)); if (!(*Ax)) { ostringstream tmp; @@ -985,12 +1004,12 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si { int t, eq, var, lag, ti_y_kmin, ti_y_kmax; double* jacob_exo ; - int row_x; + int row_x = 0; #ifdef DEBUG int col_x; -#endif +#endif int n = periods * Size; - *b = (double*)mxMalloc(n * sizeof(double)); + *b = (double*)mxMalloc(n * sizeof(double)); if (!(*b)) { ostringstream tmp; @@ -1004,22 +1023,25 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si tmp << " in Init_UMFPACK_Sparse_Simple, can't retrieve x0 vector\n"; throw FatalExceptionHandling(tmp.str()); } - *Ap = (SuiteSparse_long*)mxMalloc((n+1) * sizeof(SuiteSparse_long)); + *Ap = (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()); - } + } size_t prior_nz = IM.size() * periods; - *Ai = (SuiteSparse_long*)mxMalloc(prior_nz * sizeof(SuiteSparse_long)); + *Ai = (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()); } - *Ax = (double*)mxMalloc(prior_nz * sizeof(double)); + *Ax = (double*)mxMalloc(prior_nz * sizeof(double)); + test_mxMalloc(*Ax, __LINE__, __FILE__, __func__, prior_nz * sizeof(double)); if (!(*Ax)) { ostringstream tmp; @@ -1046,11 +1068,15 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si #ifdef DEBUG col_x = mxGetN(jacobian_exo_block[block_num]); #endif + } + else + { + jacob_exo = NULL; } #ifdef DEBUG int local_index; #endif - + bool fliped = false; bool fliped_exogenous_derivatives_updated = false; int flip_exo; @@ -1060,8 +1086,8 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si last_var = -1; it4 = IM.begin(); var = 0; - while (it4 != IM.end()) - { + while (it4 != IM.end()) + { var = it4->first.first.first; #ifdef DEBUG if (var < 0 || var >= Size) @@ -1070,7 +1096,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si tmp << " in Init_UMFPACK_Sparse, var (" << var << ") out of range\n"; throw FatalExceptionHandling(tmp.str()); } -#endif +#endif eq = it4->first.second+Size*t; #ifdef DEBUG if (eq < 0 || eq >= Size) @@ -1117,7 +1143,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si (*Ax)[NZE] = jacob_exo[k + row_x*flip_exo]; (*Ai)[NZE] = k; NZE++; - + #ifdef DEBUG if (local_index < 0 || local_index >= Size * periods) { @@ -1137,7 +1163,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si 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()); } -#endif +#endif u[k] -= jacob_exo[k + row_x*flip_exo] * x[t+y_kmin+flip_exo*nb_row_x]; } } @@ -1151,7 +1177,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si if (max_lag < lag) max_lag = lag; }*/ - + if (var < (periods+y_kmax)*Size) { ti_y_kmin = -min(t, y_kmin); @@ -1204,7 +1230,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si #endif (*b)[eq - lag * Size] += u[index] * y[index_vara[var+Size*(y_kmin+t/*+lag*/)]]; } - + } if (lag > ti_y_kmax || lag < ti_y_kmin) { @@ -1276,7 +1302,8 @@ dynSparseMatrix::Init_CUDA_Sparse_Simple(int Size, map, int> { int eq, var; - *b = (double*)mxMalloc(Size * sizeof(double)); + *b = (double*)mxMalloc(Size * sizeof(double)); + test_mxMalloc(*b, __LINE__, __FILE__, __func__, Size * sizeof(double)); if (!(*b)) { ostringstream tmp; @@ -1290,7 +1317,8 @@ dynSparseMatrix::Init_CUDA_Sparse_Simple(int Size, map, int> tmp << " in Init_CUDA_Sparse_Simple, can't retrieve x0 vector\n"; throw FatalExceptionHandling(tmp.str()); } - *Ap = (SuiteSparse_long*)mxMalloc((Size+1) * sizeof(SuiteSparse_long)); + *Ap = (SuiteSparse_long*)mxMalloc((Size+1) * sizeof(SuiteSparse_long)); + test_mxMalloc(*Ap, __LINE__, __FILE__, __func__, (Size+1) * sizeof(SuiteSparse_long)); if (!(*Ap)) { ostringstream tmp; @@ -1298,14 +1326,16 @@ dynSparseMatrix::Init_CUDA_Sparse_Simple(int Size, map, int> throw FatalExceptionHandling(tmp.str()); } size_t prior_nz = IM.size(); - *Ai = (SuiteSparse_long*)mxMalloc(prior_nz * sizeof(SuiteSparse_long)); + *Ai = (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_CUDA_Sparse, can't allocate Ai index vector\n"; throw FatalExceptionHandling(tmp.str()); } - *Ax = (double*)mxMalloc(prior_nz * sizeof(double)); + *Ax = (double*)mxMalloc(prior_nz * sizeof(double)); + test_mxMalloc(*Ax, __LINE__, __FILE__, __func__, prior_nz * sizeof(double)); if (!(*Ax)) { ostringstream tmp; @@ -1411,6 +1441,7 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, double *Host_b = (double*)mxMalloc(n * sizeof(double)); + test_mxMalloc(Host_b, __LINE__, __FILE__, __func__, n * sizeof(double)); cudaChk(cudaMalloc((void**)b, n * sizeof(double)), " in Init_Cuda_Sparse, not enought memory to allocate b vector on the graphic card\n"); double *Host_x0 = mxGetPr(x0_m); @@ -1423,18 +1454,23 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, cudaChk(cudaMalloc((void**)x0, n * sizeof(double)), " in Init_Cuda_Sparse, not enought memory to allocate x0 vector on the graphic card\n"); int* Host_Ap = (int*)mxMalloc((n+1) * sizeof(int)); + test_mxMalloc(Host_Ap, __LINE__, __FILE__, __func__, (n+1) * sizeof(int)); int* Host_Ai = (int*)mxMalloc(prior_nz * sizeof(int)); + test_mxMalloc(Host_Ai, __LINE__, __FILE__, __func__, prior_nz * sizeof(int)); double* Host_Ax = (double*)mxMalloc(prior_nz * sizeof(double)); + test_mxMalloc(Host_Ax, __LINE__, __FILE__, __func__, prior_nz * sizeof(double)); int* Host_Ai_tild, * Host_Ap_tild; if (preconditioner == 3) { Host_Ap_tild = (int*) mxMalloc((n+1)*sizeof(int)); + test_mxMalloc(Host_Ap_tild, __LINE__, __FILE__, __func__, (n+1)*sizeof(int)); Host_Ai_tild = (int*) mxMalloc(prior_nz*sizeof(int)); + test_mxMalloc(Host_Ai_tild, __LINE__, __FILE__, __func__, prior_nz*sizeof(int)); Host_Ap_tild[0] = 0; } @@ -1445,6 +1481,7 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, preconditioner_size = prior_nz; double *Host_A_tild = (double*)mxMalloc(preconditioner_size * sizeof(double)); + test_mxMalloc(Host_A_tild, __LINE__, __FILE__, __func__, preconditioner_size * sizeof(double)); map, int>, int>::iterator it4; @@ -1595,8 +1632,11 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, if (preconditioner == 3) { int* tmp_Ap_tild = (int*) mxMalloc((Size + 1) * sizeof(int) ); + test_mxMalloc(tmp_Ap_tild, __LINE__, __FILE__, __func__, (Size + 1) * sizeof(int)) ; int* tmp_Ai_tild = (int*) mxMalloc(NZE_tild * sizeof(int) ); - double* tmp_A_tild = (double*) mxMalloc(NZE_tild * sizeof(double) ); + test_mxMalloc(tmp_Ai_tild, __LINE__, __FILE__, __func__, NZE_tild * sizeof(int)); + double* tmp_A_tild = (double*) mxMalloc(NZE_tild * sizeof(double)); + test_mxMalloc(tmp_A_tild, __LINE__, __FILE__, __func__, NZE_tild * sizeof(double)); memcpy(tmp_Ap_tild, Host_Ap_tild, (Size + 1) * sizeof(int)); memcpy(tmp_Ai_tild, Host_Ai_tild, NZE_tild * sizeof(int)); memcpy(tmp_A_tild, Host_A_tild, NZE_tild * sizeof(double)); @@ -1668,10 +1708,11 @@ dynSparseMatrix::Init_CUDA_Sparse(int periods, int y_kmin, int y_kmax, int Size, void -PrintM(int n, double* Ax, mwIndex *Ap, mwIndex *Ai) +dynSparseMatrix::PrintM(int n, double* Ax, mwIndex *Ap, mwIndex *Ai) { int nnz = Ap[n]; double *A = (double*)mxMalloc(n * n * sizeof(double)); + test_mxMalloc(A, __LINE__, __FILE__, __func__, n * n * sizeof(double)); memset(A,0,n * n * sizeof(double)); int k = 0; for (int i = 0; i< n; i++) @@ -1849,23 +1890,36 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map, int>, int>::iterator it4; NonZeroElem *first; pivot = (int *) mxMalloc(Size*periods*sizeof(int)); - pivot_save = (int *) mxMalloc(Size*periods*sizeof(int)); - pivotk = (int *) mxMalloc(Size*periods*sizeof(int)); - pivotv = (double *) mxMalloc(Size*periods*sizeof(double)); - pivotva = (double *) mxMalloc(Size*periods*sizeof(double)); - b = (int *) mxMalloc(Size*periods*sizeof(int)); - line_done = (bool *) mxMalloc(Size*periods*sizeof(bool)); + test_mxMalloc(pivot, __LINE__, __FILE__, __func__, Size*periods*sizeof(int)); + pivot_save = (int *) mxMalloc(Size*periods*sizeof(int)); + test_mxMalloc(pivot_save, __LINE__, __FILE__, __func__, Size*periods*sizeof(int)); + pivotk = (int *) mxMalloc(Size*periods*sizeof(int)); + test_mxMalloc(pivotk, __LINE__, __FILE__, __func__, Size*periods*sizeof(int)); + pivotv = (double *) mxMalloc(Size*periods*sizeof(double)); + test_mxMalloc(pivotv, __LINE__, __FILE__, __func__, Size*periods*sizeof(double)); + pivotva = (double *) mxMalloc(Size*periods*sizeof(double)); + test_mxMalloc(pivotva, __LINE__, __FILE__, __func__, Size*periods*sizeof(double)); + b = (int *) mxMalloc(Size*periods*sizeof(int)); + test_mxMalloc(b, __LINE__, __FILE__, __func__, Size*periods*sizeof(int)); + line_done = (bool *) mxMalloc(Size*periods*sizeof(bool)); + test_mxMalloc(line_done, __LINE__, __FILE__, __func__, Size*periods*sizeof(bool)); mem_mngr.init_CHUNK_BLCK_SIZE(u_count); g_save_op = NULL; g_nop_all = 0; i = (periods+y_kmax+1)*Size*sizeof(NonZeroElem *); - FNZE_R = (NonZeroElem **) mxMalloc(i); - FNZE_C = (NonZeroElem **) mxMalloc(i); - NonZeroElem **temp_NZE_R = (NonZeroElem **) mxMalloc(i); - NonZeroElem **temp_NZE_C = (NonZeroElem **) mxMalloc(i); + FNZE_R = (NonZeroElem **) mxMalloc(i); + test_mxMalloc(FNZE_R, __LINE__, __FILE__, __func__, i); + FNZE_C = (NonZeroElem **) mxMalloc(i); + test_mxMalloc(FNZE_C, __LINE__, __FILE__, __func__, i); + NonZeroElem **temp_NZE_R = (NonZeroElem **) mxMalloc(i); + test_mxMalloc(*temp_NZE_R, __LINE__, __FILE__, __func__, i); + NonZeroElem **temp_NZE_C = (NonZeroElem **) mxMalloc(i); + test_mxMalloc(*temp_NZE_C, __LINE__, __FILE__, __func__, i); i = (periods+y_kmax+1)*Size*sizeof(int); - NbNZRow = (int *) mxMalloc(i); - NbNZCol = (int *) mxMalloc(i); + NbNZRow = (int *) mxMalloc(i); + test_mxMalloc(NbNZRow, __LINE__, __FILE__, __func__, i); + NbNZCol = (int *) mxMalloc(i); + test_mxMalloc(NbNZCol, __LINE__, __FILE__, __func__, i); for (int i = 0; i < periods*Size; i++) { @@ -2029,7 +2083,9 @@ dynSparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, t_save_op_s *save_op_s, *save_opa_s, *save_opaa_s; int *diff1, *diff2; diff1 = (int *) mxMalloc(nop*sizeof(int)); + test_mxMalloc(diff1, __LINE__, __FILE__, __func__, nop*sizeof(int)); diff2 = (int *) mxMalloc(nop*sizeof(int)); + test_mxMalloc(diff2, __LINE__, __FILE__, __func__, nop*sizeof(int)); int max_save_ops_first = -1; j = i = 0; while (i < nop4 && OK) @@ -2180,8 +2236,10 @@ dynSparseMatrix::complete(int beg_t, int Size, int periods, int *b) int size_of_save_code = (1+y_kmax)*Size*(Size+1+4)/2*4; save_code = (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 = (int *) mxMalloc(size_of_diff*sizeof(int)); + test_mxMalloc(diff, __LINE__, __FILE__, __func__, size_of_diff*sizeof(int)); cal_y = y_size*y_kmin; i = (beg_t+1)*Size-1; @@ -2439,6 +2497,7 @@ dynSparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int perio mexPrintf("row(2)=%d\n", row); double *B; B = (double *) mxMalloc(row*sizeof(double)); + test_mxMalloc(B, __LINE__, __FILE__, __func__, row*sizeof(double)); for (int i = 0; i < row; i++) SaveResult >> B[i]; SaveResult.close(); @@ -3293,8 +3352,11 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do #else double *Control, *Info, *res; Control = (double*)mxMalloc(UMFPACK_CONTROL * sizeof(double)); + test_mxMalloc(Control, __LINE__, __FILE__, __func__, UMFPACK_CONTROL * sizeof(double)); Info = (double*)mxMalloc(UMFPACK_INFO * sizeof(double)); + test_mxMalloc(Info, __LINE__, __FILE__, __func__, UMFPACK_INFO * sizeof(double)); res = (double*)mxMalloc(n * sizeof(double)); + test_mxMalloc(res, __LINE__, __FILE__, __func__, n * sizeof(double)); #endif umfpack_dl_defaults(Control); @@ -3397,7 +3459,7 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do y[eq+it_*y_size] += slowc_l * yy; } } - + mxFree(Ap); mxFree(Ai); mxFree(Ax); @@ -3419,8 +3481,11 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do #else double *Control, *Info, *res; Control = (double*)mxMalloc(UMFPACK_CONTROL * sizeof(double)); + test_mxMalloc(Control, __LINE__, __FILE__, __func__, UMFPACK_CONTROL * sizeof(double)); Info = (double*)mxMalloc(UMFPACK_INFO * sizeof(double)); + test_mxMalloc(Info, __LINE__, __FILE__, __func__, UMFPACK_INFO * sizeof(double)); res = (double*)mxMalloc(n * sizeof(double)); + test_mxMalloc(res, __LINE__, __FILE__, __func__, n * sizeof(double)); #endif umfpack_dl_defaults(Control); @@ -3503,8 +3568,11 @@ dynSparseMatrix::Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double s #else double *Control, *Info, *res; Control = (double*)mxMalloc(UMFPACK_CONTROL * sizeof(double)); + test_mxMalloc(Control, __LINE__, __FILE__, __func__, UMFPACK_CONTROL * sizeof(double)); Info = (double*)mxMalloc(UMFPACK_INFO * sizeof(double)); + test_mxMalloc(Info, __LINE__, __FILE__, __func__, UMFPACK_INFO * sizeof(double)); res = (double*)mxMalloc(n * sizeof(double)); + test_mxMalloc(res, __LINE__, __FILE__, __func__, n * sizeof(double)); #endif void *Symbolic, *Numeric ; umfpack_dl_defaults (Control) ; @@ -3561,6 +3629,7 @@ printM(int n,double *Ax, int* Ap, int* Ai, cusparseMatDescr_t descrA, cusparseH cusparseChk(cusparseDcsr2dense(cusparse_handle, n, n, descrA, Ax, Ap,Ai, A_dense, n), "cusparseDcsr2dense has failed\n"); double *A_dense_hoste = (double*)mxMalloc(n * n * sizeof(double)); + test_mxMalloc(A_dense_hoste, __LINE__, __FILE__, __func__, n * n * sizeof(double)); cudaChk(cudaMemcpy(A_dense_hoste, A_dense, n * n * sizeof(double),cudaMemcpyDeviceToHost), " cudaMemcpy(A_dense_hoste, A_dense) has failed\n"); mexPrintf("----------------------\n"); mexPrintf("FillMode=%d, IndexBase=%d, MatType=%d, DiagType=%d\n",cusparseGetMatFillMode(descrA), cusparseGetMatIndexBase(descrA), cusparseGetMatType(descrA), cusparseGetMatDiagType(descrA)); @@ -3731,6 +3800,7 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, int periods = n / Size; double * tmp_vect_host = (double*)mxMalloc(n * sizeof(double)); + test_mxMalloc(tmp_vect_host, __LINE__, __FILE__, __func__, n * sizeof(double)); cublasChk(cublasDnrm2(cublas_handle, n,b, 1, &bnorm), " in Solve_Cuda_BiCGStab, cublasDnrm2(b) has failed\n"); @@ -3856,10 +3926,15 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, // we have to transpose it to get a CSR format used by CUDA mwIndex* Awi, *Awp; double* A_tild_host = (double*)mxMalloc(nnz*sizeof(double)); + test_mxMalloc(A_tild_host, __LINE__, __FILE__, __func__, nnz*sizeof(double)); Awi = (mwIndex*)mxMalloc(nnz * sizeof(mwIndex)); + test_mxMalloc(Awi, __LINE__, __FILE__, __func__, nnz * sizeof(mwIndex)); Awp = (mwIndex*)mxMalloc((n + 1) * sizeof(mwIndex)); + test_mxMalloc(Awp, __LINE__, __FILE__, __func__, (n + 1) * sizeof(mwIndex)); int* Aii = (int*)mxMalloc(nnz * sizeof(int)); + test_mxMalloc(Aii, __LINE__, __FILE__, __func__, nnz * sizeof(int)); int* Aip = (int*)mxMalloc((n + 1) * sizeof(int)); + test_mxMalloc(Aip, __LINE__, __FILE__, __func__, (n + 1) * sizeof(int)); cudaChk(cudaMemcpy(A_tild_host, A_tild, nnz*sizeof(double), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy A_tild_host = A_tild has failed\n"); cudaChk(cudaMemcpy(Aii, Ai, nnz*sizeof(int), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy Aii = Ai has failed\n"); cudaChk(cudaMemcpy(Aip, Ap, (n+1)*sizeof(int), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy Aip = Ai has failed\n"); @@ -3941,7 +4016,9 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, mwIndex* Wi = mxGetIr(W); mwIndex* Wp = mxGetJc(W); int *Wii = (int*)mxMalloc(nnz * sizeof(int)); + test_mxMalloc(Wii, __LINE__, __FILE__, __func__, nnz * sizeof(int)); int *Wip = (int*)mxMalloc((n + 1) * sizeof(int)); + test_mxMalloc(Wip, __LINE__, __FILE__, __func__, (n + 1) * sizeof(int)); for (int i = 0; i < nnz; i++) Wii[i] = Wi[i]; for (int i = 0; i < n + 1; i++) @@ -3968,10 +4045,15 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, { mwIndex* Aowi, *Aowp; double* A_host = (double*)mxMalloc(nnz*sizeof(double)); + test_mxMalloc(A_host, __LINE__, __FILE__, __func__, nnz*sizeof(double)); Aowi = (mwIndex*)mxMalloc(nnz * sizeof(mwIndex)); + test_mxMalloc(Aowi, __LINE__, __FILE__, __func__, nnz * sizeof(mwIndex)); Aowp = (mwIndex*)mxMalloc((n + 1) * sizeof(mwIndex)); + test_mxMalloc(Aowp, __LINE__, __FILE__, __func__, (n + 1) * sizeof(mwIndex)); int* Aoii = (int*)mxMalloc(nnz * sizeof(int)); + test_mxMalloc(Aoii, __LINE__, __FILE__, __func__, nnz * sizeof(int)); int* Aoip = (int*)mxMalloc((n + 1) * sizeof(int)); + test_mxMalloc(Aoip, __LINE__, __FILE__, __func__, (n + 1) * sizeof(int)); cudaChk(cudaMemcpy(A_host, Ax, nnz*sizeof(double), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy A_tild_host = A_tild has failed\n"); cudaChk(cudaMemcpy(Aoii, Ai, nnz*sizeof(int), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy Aii = Ai_tild has failed\n"); cudaChk(cudaMemcpy(Aoip, Ap, (n+1)*sizeof(int), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy Aip = Ap_tild has failed\n"); @@ -3996,10 +4078,15 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, // we have to transpose it to get a CSR format used by CUDA mwIndex* Awi, *Awp; double* A_tild_host = (double*)mxMalloc(nnz_tild*sizeof(double)); + test_mxMalloc(A_tild_host, __LINE__, __FILE__, __func__, nnz_tild*sizeof(double)); Awi = (mwIndex*)mxMalloc(nnz_tild * sizeof(mwIndex)); + test_mxMalloc(Awi, __LINE__, __FILE__, __func__, nnz_tild * sizeof(mwIndex)); Awp = (mwIndex*)mxMalloc((Size + 1) * sizeof(mwIndex)); + test_mxMalloc(Awp, __LINE__, __FILE__, __func__, (Size + 1) * sizeof(mwIndex)); int* Aii = (int*)mxMalloc(nnz_tild * sizeof(int)); + test_mxMalloc(Aii, __LINE__, __FILE__, __func__, nnz_tild * sizeof(int)); int* Aip = (int*)mxMalloc((Size + 1) * sizeof(int)); + test_mxMalloc(Aip, __LINE__, __FILE__, __func__, (Size + 1) * sizeof(int)); cudaChk(cudaMemcpy(A_tild_host, A_tild, nnz_tild*sizeof(double), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy A_tild_host = A_tild has failed\n"); cudaChk(cudaMemcpy(Aii, Ai_tild, nnz_tild*sizeof(int), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy Aii = Ai_tild has failed\n"); cudaChk(cudaMemcpy(Aip, Ap_tild, (Size+1)*sizeof(int), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy Aip = Ap_tild has failed\n"); @@ -4060,8 +4147,11 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, Q_nnz = Qjw_host[Size]; mexPrintf("Q_nnz=%d\n",Q_nnz); int *Qi_host = (int*)mxMalloc(Q_nnz * periods * sizeof(int)); + test_mxMalloc(Qi_host, __LINE__, __FILE__, __func__, Q_nnz * periods * sizeof(int)); double *Q_x_host = (double*)mxMalloc(Q_nnz * periods * sizeof(double)); + test_mxMalloc(Q_x_host, __LINE__, __FILE__, __func__, Q_nnz * periods * sizeof(double)); int *Qj_host = (int*)mxMalloc((n + 1) * sizeof(int)); + test_mxMalloc(Qj_host, __LINE__, __FILE__, __func__, (n + 1) * sizeof(int)); for (int t = 0; t < periods; t++) { for (int i = 0; i < Q_nnz; i++) @@ -4120,8 +4210,11 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, double* Px_host = mxGetPr(P); P_nnz = Pjw_host[Size]; int *Pi_host = (int*)mxMalloc(P_nnz * periods * sizeof(int)); + test_mxMalloc(Pi_host, __LINE__, __FILE__, __func__, P_nnz * periods * sizeof(int)); double *P_x_host = (double*)mxMalloc(P_nnz * periods * sizeof(double)); + test_mxMalloc(P_x_host, __LINE__, __FILE__, __func__, P_nnz * periods * sizeof(double)); int *Pj_host = (int*)mxMalloc((n + 1) * sizeof(int)); + test_mxMalloc(Pj_host, __LINE__, __FILE__, __func__, (n + 1) * sizeof(int)); for (int t = 0; t < periods; t++) { for (int i = 0; i < P_nnz; i++) @@ -4202,8 +4295,11 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild, int U_nnz = Ujw_host[Size]; double *pW = (double*)mxMalloc((L_nnz + U_nnz - Size) * periods * sizeof(double)); + test_mxMalloc(pW, __LINE__, __FILE__, __func__, (L_nnz + U_nnz - Size) * periods * sizeof(double)); int *Wi = (int*)mxMalloc((L_nnz + U_nnz - Size) * periods * sizeof(int)); + test_mxMalloc(Wi, __LINE__, __FILE__, __func__, (L_nnz + U_nnz - Size) * periods * sizeof(int)); int *Wj = (int*)mxMalloc((n + 1) * sizeof(int)); + test_mxMalloc(Wj, __LINE__, __FILE__, __func__, (n + 1) * sizeof(int)); Wj[0] = 0; W_nnz = 0; for (int t = 0; t < periods; t++) @@ -4865,7 +4961,10 @@ 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; + mxArray *L1, *U1, *Diag; + L1 = NULL; + U1 = NULL; + Diag = NULL; mxArray *rhs0[4]; if (preconditioner == 0) @@ -4919,7 +5018,8 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou mxDestroyArray(Setup); } double flags = 2; - mxArray *z; + mxArray *z; + z = NULL; 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); @@ -5151,10 +5251,15 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i double piv_abs; NonZeroElem **bc; bc = (NonZeroElem **) mxMalloc(Size*sizeof(*bc)); + test_mxMalloc(bc, __LINE__, __FILE__, __func__, Size*sizeof(*bc)); piv_v = (double *) mxMalloc(Size*sizeof(double)); + test_mxMalloc(piv_v, __LINE__, __FILE__, __func__, Size*sizeof(double)); pivj_v = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(pivj_v, __LINE__, __FILE__, __func__, Size*sizeof(int)); pivk_v = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(pivk_v, __LINE__, __FILE__, __func__, Size*sizeof(int)); NR = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(NR, __LINE__, __FILE__, __func__, Size*sizeof(int)); for (int i = 0; i < Size; i++) { @@ -5419,12 +5524,17 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo int tbreak = 0, last_period = periods; piv_v = (double *) mxMalloc(Size*sizeof(double)); + test_mxMalloc(piv_v, __LINE__, __FILE__, __func__, Size*sizeof(double)); pivj_v = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(pivj_v, __LINE__, __FILE__, __func__, Size*sizeof(int)); pivk_v = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(pivk_v, __LINE__, __FILE__, __func__, Size*sizeof(int)); NR = (int *) mxMalloc(Size*sizeof(int)); + test_mxMalloc(NR, __LINE__, __FILE__, __func__, Size*sizeof(int)); //clock_t time00 = clock(); NonZeroElem **bc; bc = (NonZeroElem **) mxMalloc(Size*sizeof(first)); + test_mxMalloc(bc, __LINE__, __FILE__, __func__, Size*sizeof(first)); for (int t = 0; t < periods; t++) { @@ -5443,6 +5553,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo save_op = NULL; }*/ save_op = (int *) mxMalloc(nop*sizeof(int)); + test_mxMalloc(save_op, __LINE__, __FILE__, __func__, nop*sizeof(int)); nopa = nop; } nop = 0; @@ -6111,6 +6222,7 @@ dynSparseMatrix::Check_and_Correct_Previous_Iteration(int block_num, int y_size, { double *p = (double*)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++) @@ -6337,7 +6449,9 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in mxFree(Ai_save); mxFree(Ax_save); Ai_save = (SuiteSparse_long*)mxMalloc(Ap[size] * sizeof(SuiteSparse_long)); + test_mxMalloc(Ai_save, __LINE__, __FILE__, __func__, Ap[size] * sizeof(SuiteSparse_long)); Ax_save = (double*)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)); @@ -6420,15 +6534,21 @@ void dynSparseMatrix::Simulate_Newton_One_Boundary(const bool forward) { g1 = (double *) mxMalloc(size*size*sizeof(double)); + test_mxMalloc(g1, __LINE__, __FILE__, __func__, size*size*sizeof(double)); r = (double *) 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)) { Ap_save = (SuiteSparse_long*)mxMalloc((size + 1) * sizeof(SuiteSparse_long)); + test_mxMalloc(Ap_save, __LINE__, __FILE__, __func__, (size + 1) * sizeof(SuiteSparse_long)); Ap_save[size] = 0; Ai_save = (SuiteSparse_long*)mxMalloc(1 * sizeof(SuiteSparse_long)); + test_mxMalloc(Ai_save, __LINE__, __FILE__, __func__, 1 * sizeof(SuiteSparse_long)); Ax_save = (double*)mxMalloc(1 * sizeof(double)); + test_mxMalloc(Ax_save, __LINE__, __FILE__, __func__, 1 * sizeof(double)); b_save = (double*)mxMalloc((size) * sizeof(SuiteSparse_long)); + test_mxMalloc(b_save, __LINE__, __FILE__, __func__, (size) * sizeof(SuiteSparse_long)); } if (steady_state) { @@ -6491,7 +6611,7 @@ dynSparseMatrix::preconditioner_print_out(string s, int preconditioner, bool ss) case 1: if (ss) tmp.append("incomplet lutp on static jacobian"); - else + else tmp.append("incomplet lu0 on dynamic jacobian"); break; case 2: @@ -6540,7 +6660,8 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter > 0)) { if (iter == 0 || fabs(slowc_save) < 1e-8) - { + { + mexPrintf("res1 = %f, res2 = %f g0 = %f iter = %d\n", res1, res2, g0, iter); for (int j = 0; j < y_size; j++) { ostringstream res; @@ -6781,6 +6902,7 @@ dynSparseMatrix::fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_p mexPrintf("fixe_u : alloc(%d double)\n", u_count_alloc); #endif (*u) = (double *) 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 diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index 82fd5a3ab..23090797a 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/mex/sources/bytecode/SparseMatrix.hh @@ -146,7 +146,7 @@ private: 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); 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_); @@ -245,7 +245,7 @@ protected: int u_count_alloc, u_count_alloc_save; vector jac; double *jcb; - double slowc_save, prev_slowc_save, markowitz_c; + double slowc_save, prev_slowc_save, markowitz_c; int y_decal; int *index_equa; int u_count, tbreak_g; diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 25260c82b..eee79bf45 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -91,9 +91,9 @@ Get_Argument(const mxArray *prhs) #endif -//#include -#include - +//#include +#include + #ifdef CUDA int @@ -449,7 +449,8 @@ main(int nrhs, const char *prhs[]) char *plhs[1]; load_global((char *) prhs[1]); #endif - mxArray *pfplan_struct = NULL; + mxArray *pfplan_struct = NULL; + ErrorMsg error_msg; size_t i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0; size_t steady_row_y, steady_col_y; int y_kmin = 0, y_kmax = 0, y_decal = 0; @@ -527,7 +528,7 @@ main(int nrhs, const char *prhs[]) DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str()); } int nb_periods = mxGetM(date_str) * mxGetN(date_str); - + mxArray* constrained_vars_ = mxGetField(extended_path_struct, 0, "constrained_vars_"); if (constrained_vars_ == NULL) { @@ -556,7 +557,7 @@ main(int nrhs, const char *prhs[]) tmp.insert(0,"The extended_path description structure does not contain the member: "); DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str()); } - + mxArray* shock_var_ = mxGetField(extended_path_struct, 0, "shock_vars_"); if (shock_var_ == NULL) { @@ -621,7 +622,7 @@ main(int nrhs, const char *prhs[]) table_conditional_global[i] = vector_conditional_local; } } - + vector_table_conditional_local_type vv3 = table_conditional_global[0]; for (int i = 0; i < nb_constrained; i++) { @@ -632,6 +633,7 @@ main(int nrhs, const char *prhs[]) double *specific_constrained_int_date_ = mxGetPr(mxGetCell(constrained_int_date_, i)); int nb_local_periods = mxGetM(Array_constrained_paths_) * mxGetN(Array_constrained_paths_); int* constrained_int_date = (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; @@ -642,7 +644,7 @@ main(int nrhs, const char *prhs[]) oss << nb_local_periods; string tmp1 = oss.str(); tmp.append(tmp1); - tmp.append(")"); + tmp.append(")"); DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str()); } (sconditional_extended_path[i]).per_value.resize(nb_local_periods); @@ -685,7 +687,7 @@ main(int nrhs, const char *prhs[]) oss << nb_local_periods; string tmp1 = oss.str(); tmp.append(tmp1); - tmp.append(")"); + tmp.append(")"); DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str()); } (sextended_path[i]).per_value.resize(nb_local_periods); @@ -713,7 +715,7 @@ main(int nrhs, const char *prhs[]) dates.push_back(string(buf));//string(Dates[i]); mxFree(buf); } - } + } if (plan.length()>0) { mxArray* plan_struct = mexGetVariable("base", plan.c_str()); @@ -1057,14 +1059,19 @@ main(int nrhs, const char *prhs[]) #else if (stack_solve_algo == 7 && !steady_state) DYN_MEX_FUNC_ERR_MSG_TXT("bytecode has not been compiled with CUDA option. Bytecode Can't use options_.stack_solve_algo=7\n"); -#endif +#endif size_t size_of_direction = col_y*row_y*sizeof(double); double *y = (double *) mxMalloc(size_of_direction); + error_msg.test_mxMalloc(y, __LINE__, __FILE__, __func__, size_of_direction); double *ya = (double *) mxMalloc(size_of_direction); + error_msg.test_mxMalloc(ya, __LINE__, __FILE__, __func__, size_of_direction); direction = (double *) mxMalloc(size_of_direction); - memset(direction, 0, 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);*/ double *x = (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]); @@ -1075,7 +1082,7 @@ main(int nrhs, const char *prhs[]) ya[i] = double (yd[i]); } size_t y_size = row_y; - size_t nb_row_x = row_x; + size_t nb_row_x = row_x; clock_t t0 = clock(); Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms, steady_state,