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

Corrects several bugs related to bytecode: 

- Memory allocation is checked
- The amount of memory allocated for conditional forecast is corrected
time-shift
Ferhat Mihoubi 2017-03-24 23:15:25 +01:00 committed by Stéphane Adjemian (Charybdis)
parent 8a18e87d98
commit 857fc3c4f4
10 changed files with 342 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -24,8 +24,20 @@
#include <iostream>
#include <sstream>
#include <map>
#include <stack>
#define BYTE_CODE
#include "CodeInterpreter.hh"
#define _USE_MATH_DEFINES
#include <math.h>
#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif
#ifndef M_SQRT2
#define M_SQRT2 1.41421356237309504880
#endif
#ifdef DEBUG_EX
# include <math.h>
# 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 <limits>
#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());
}
}
};

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

@ -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)
{

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

@ -63,13 +63,13 @@ protected:
double solve_tolf;
bool GaussSeidel;
map<pair<pair<int, int>, 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_type> Block_Contain;
int size;

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

@ -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<double>(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)

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

@ -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

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

@ -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]);
}

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

@ -19,15 +19,16 @@
#ifndef MEM_MNGR_HH_INCLUDED
#define MEM_MNGR_HH_INCLUDED
#include "ErrorHandling.hh"
#include <vector>
#include <fstream>
#ifndef DEBUG_EX
# include <dynmex.h>
#else
# include "mex_interface.hh"
#endif
using namespace std;
#endif
//using namespace std;
struct NonZeroElem
{
@ -41,7 +42,7 @@ typedef vector<NonZeroElem *> 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<NonZeroElem *> NZE_Mem_Allocated;
int swp_f_b;
fstream SaveCode_swp;
string filename;
string filename_mem;
};
#endif

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

@ -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<char *>(&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<char *>(&index_equa[j]), sizeof(*index_equa));
}
@ -641,25 +643,38 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
int i, eq, var, lag;
map<pair<pair<int, int>, 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<pair<pair<int, int>, 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<pair<pair<int, int>, 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<pair<pair<int, int>, 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<pair<pair<int, int>, 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<pair<pair<int, int>, 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<pair<pair<int, int>, 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<pair<pair<int, int>, 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<pair
map<pair<pair<int, int>, 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

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

@ -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<double *> 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;

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

@ -91,9 +91,9 @@ Get_Argument(const mxArray *prhs)
#endif
//#include <windows.h>
#include <stdio.h>
//#include <windows.h>
#include <stdio.h>
#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,