Skip to content
GitLab
Commit 77b60969 authored by Ferhat Mihoubi's avatar Ferhat Mihoubi Committed by Sébastien Villemot
Browse files

- extends the mex_interface: bytecode's debuging purpose 

- correction of memory leaks in bytecode
parent 60f76786
Expand all Hide whitespace changes
Inline Side-by-side
mex/sources/bytecode/Interpreter.cc
View file @ 77b60969
This diff is collapsed.
mex/sources/bytecode/Interpreter.hh
View file @ 77b60969
......@@ -47,22 +47,23 @@ typedef code_liste_type::const_iterator it_code_type;
class Interpreter : SparseMatrix
{
private:
#ifndef DEBUG_EX
vector<mxArray*> jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block;
#endif
ExpressionType EQN_type;
char *P_endo_names, *P_exo_names, *P_param_names;
unsigned int nb_endo, nb_exo, nb_param;
unsigned int endo_name_length, exo_name_length, param_name_length;
unsigned int EQN_equation, EQN_block, EQN_block_number;
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
int EQN_lag1, EQN_lag2, EQN_lag3;
it_code_type it_code_expr;
protected:
double pow1(double a, double b, bool evaluate);
double log1(double a, bool evaluate);
double log10_1(double a, bool evaluate);
string remove_white(string str);
string add_underscore_to_fpe(string str);
string get_variable(SymbolType variable_type, unsigned int variable_num);
string error_location(bool evaluate);
double pow1(double a, double b, bool evaluate, bool steady_state);
double log1(double a, bool evaluate, bool steady_state);
double log10_1(double a, bool evaluate, bool steady_state);
/*string remove_white(string str);*/
string add_underscore_to_fpe(const string &str);
string get_variable(const SymbolType variable_type, const unsigned int variable_num);
string error_location(bool evaluate, bool steady_state);
void compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state);
string print_expression(it_code_type it_code, bool evaluate);
void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
......@@ -73,7 +74,6 @@ private:
vector<Block_contain_type> Block_Contain;
code_liste_type code_liste;
it_code_type it_code;
stack<double> Stack;
int Block_Count, Per_u_, Per_y_;
int it_, nb_row_x, nb_row_xd, maxit_, size_of_direction;
double *g1, *r;
......@@ -93,12 +93,10 @@ public:
int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg);
bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks);
#ifndef DEBUG_EX
inline mxArray* get_jacob(int block_num) {return jacobian_block[block_num];};
inline mxArray* get_jacob_exo(int block_num) {return jacobian_exo_block[block_num];};
inline mxArray* get_jacob_exo_det(int block_num) {return jacobian_det_exo_block[block_num];};
inline mxArray* get_jacob_other_endo(int block_num) {return jacobian_other_endo_block[block_num];};
#endif
};
#endif
mex/sources/bytecode/Mem_Mngr.cc
View file @ 77b60969
......@@ -84,7 +84,10 @@ Mem_Mngr::mxMalloc_NZE()
mexPrintf("Not enough memory available\n");
mexEvalString("drawnow;");
}
NZE_Mem_add = (NonZeroElem **) mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to redefine the size of pointer on the memory*/
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*/
if (!NZE_Mem_add)
{
mexPrintf("Not enough memory available\n");
......@@ -121,6 +124,7 @@ Mem_Mngr::Free_All()
mxFree(NZE_Mem_Allocated.back());
NZE_Mem_Allocated.pop_back();
}
mxFree(NZE_Mem_add);
if (NZE_Mem_add)
mxFree(NZE_Mem_add);
init_Mem();
}
mex/sources/bytecode/SparseMatrix.cc
View file @ 77b60969
......@@ -478,9 +478,9 @@ SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>,
map<pair<pair<int, int>, int>, int>::iterator it4;
for (i = 0; i < y_size*(periods+y_kmin); i++)
ya[i] = y[i];
#if DEBUG
//#if DEBUG
unsigned int max_nze = mxGetNzmax(A_m);
#endif
//#endif
unsigned int NZE = 0;
int last_var = 0;
for (i = 0; i < Size; i++)
......@@ -499,9 +499,9 @@ SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>,
eq = it4->first.second;
int index = it4->second;
#if DEBUG
if (index<0 || index >= u_count_alloc)
if (index<0 || index >= u_count_alloc || index > Size + Size*Size)
{
mexPrintf("index (%d) out of range for u vector (0)\n",index);
mexPrintf("index (%d) out of range for u vector (0) max %d allocated %d\n",index, Size+Size*Size, u_count_alloc);
mexErrMsgTxt("end of bytecode\n");
}
if (NZE >= max_nze)
......@@ -514,24 +514,19 @@ SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>,
Ai[NZE] = eq;
NZE++;
#if DEBUG
if ((index+lag*u_count_init) < 0 || (index+lag*u_count_init) >= u_count_alloc)
{
mexPrintf("index (%d) out of range for u vector (1)\n",index+lag*u_count_init);
mexErrMsgTxt("end of bytecode\n");
}
if (eq < 0 || eq >= (Size*periods))
if (eq < 0 || eq >= Size)
{
mexPrintf("index (%d) out of range for b vector (0)\n",eq);
mexErrMsgTxt("end of bytecode\n");
}
if (var+Size*(y_kmin+t+lag) < 0 || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
if (var < 0 || var >= Size)
{
mexPrintf("index (%d) out of range for index_vara vector (0)\n",var+Size*(y_kmin+t+lag));
mexPrintf("index (%d) out of range for index_vara vector (0)\n",var);
mexErrMsgTxt("end of bytecode\n");
}
if (index_vara[var+Size*(y_kmin+t+lag)] < 0 || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
if (index_vara[var] < 0 || index_vara[var] >= y_size)
{
mexPrintf("index (%d) out of range for y vector max=%d (0)\n",index_vara[var+Size*(y_kmin+t+lag)], y_size*(periods+y_kmin+y_kmax));
mexPrintf("index (%d) out of range for y vector max=%d (0)\n",index_vara[var], y_size);
mexErrMsgTxt("end of bytecode\n");
}
#endif
......@@ -1358,7 +1353,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
mexPrintf("Can't allocate b_m matrix in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
A_m = mxCreateSparse(Size, Size, IM_i.size()*2, mxREAL);
A_m = mxCreateSparse(Size, Size, min(int(IM_i.size()*2), Size*Size), mxREAL);
if (!A_m)
{
mexPrintf("Can't allocate A_m matrix in LU solver\n");
......@@ -1612,11 +1607,11 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
mxFree(bc);
}
else if (solve_algo == 1 || solve_algo == 4)
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc);
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
else if (solve_algo == 2)
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck);
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_);
else if (solve_algo == 3)
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck);
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_);
return true;
}
......@@ -1732,29 +1727,38 @@ SparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size,
}
void
SparseMatrix::Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l)
SparseMatrix::Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l, bool is_two_boundaries, int it_)
{
int n = mxGetM(A_m);
mxArray *z;
mxArray *rhs[2];
rhs[0] = A_m;
rhs[1] = b_m;
mexCallMATLAB(1,&z,2, rhs, "mldivide");
mexCallMATLAB(1, &z, 2, rhs, "mldivide");
double *res = mxGetPr(z);
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc_l * yy;
}
if (is_two_boundaries)
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc_l * yy;
}
else
for (int i = 0; i < n; i++)
{
int eq = index_vara[i/*+Size*it_*/];
double yy = - (res[i] + y[eq+it_*y_size]);
direction[eq] = yy;
y[eq+it_*y_size] += slowc_l * yy;
}
mxDestroyArray(A_m);
mxDestroyArray(b_m);
mxDestroyArray(z);
}
void
SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block)
SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_)
{
int n = mxGetM(A_m);
/*[L1, U1]=luinc(g1a,luinc_tol);*/
......@@ -1809,13 +1813,22 @@ SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double sl
else
{
double *res = mxGetPr(z);
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc * yy;
}
if (is_two_boundaries)
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc * yy;
}
else
for (int i = 0; i < n; i++)
{
int eq = index_vara[i];
double yy = - (res[i] + y[eq+it_*y_size]);
direction[eq] = yy;
y[eq] += slowc * yy;
}
}
mxDestroyArray(A_m);
mxDestroyArray(b_m);
......@@ -1825,7 +1838,7 @@ SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double sl
void
SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block)
SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_)
{
int n = mxGetM(A_m);
/*[L1, U1]=luinc(g1a,luinc_tol);*/
......@@ -1878,13 +1891,22 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double
else
{
double *res = mxGetPr(z);
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc * yy;
}
if (is_two_boundaries)
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc * yy;
}
else
for (int i = 0; i < n; i++)
{
int eq = index_vara[i];
double yy = - (res[i] + y[eq+it_*y_size]);
direction[eq] = yy;
y[eq] += slowc * yy;
}
}
mxDestroyArray(A_m);
mxDestroyArray(b_m);
......@@ -1894,7 +1916,7 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double
int
SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number, int stack_solve_algo)
SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number, int stack_solve_algo, int endo_name_length, char *P_endo_names)
{
/*Triangularisation at each period of a block using a simple gaussian Elimination*/
t_save_op_s *save_op_s;
......@@ -1925,8 +1947,12 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
{
if (iter == 0)
{
for (j = 0; j < y_size; j++)
for (j = 0; j < y_size; j++)
{
ostringstream res;
for (unsigned int i = 0; i < endo_name_length; i++)
if (P_endo_names[2*(j+i*nb_endo)] != ' ')
res << P_endo_names[2*(j+i*nb_endo)];
bool select = false;
for (int i = 0; i < Size; i++)
if (j == index_vara[i])
......@@ -1935,9 +1961,9 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
break;
}
if (select)
mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
mexPrintf("-> variable %s (%d) at time %d = %f direction = %f\n", res.str().c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
else
mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
mexPrintf(" variable %s (%d) at time %d = %f direction = %f\n", res.str().c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
}
mexPrintf("res1=%5.10f\n", res1);
mexPrintf("The initial values of endogenous variables are too far from the solution.\n");
......@@ -1950,6 +1976,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
{
for (j = 0; j < y_size; j++)
{
ostringstream res;
for (unsigned int i = 0; i < endo_name_length; i++)
if (P_endo_names[2*(j+i*nb_endo)] != ' ')
res << P_endo_names[2*(j+i*nb_endo)];
bool select = false;
for (int i = 0; i < Size; i++)
if (j == index_vara[i])
......@@ -1958,16 +1988,16 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
break;
}
if (select)
mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
mexPrintf("-> variable %s (%d) at time %d = %f direction = %f\n", res.str().c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
else
mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
mexPrintf(" variable %s (%d) at time %d = %f direction = %f\n", res.str().c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
}
mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx);
mexEvalString("drawnow;");
filename += " stopped";
mexErrMsgTxt(filename.c_str());
}
if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)) && (stack_solve_algo == 1 || stack_solve_algo == 5))
if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)) && (stack_solve_algo == 4 || stack_solve_algo == 5))
{
if (try_at_iteration == 0)
{
......@@ -2002,7 +2032,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
slowc_save /= 1.1;
}
if (print_it)
mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save);
{
if (isnan(res1) || isinf(res1))
mexPrintf("Error: The model cannot be evaluated, trying to correct it using slowc=%f\n", slowc_save);
else
mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save);
}
for (i = 0; i < y_size*(periods+y_kmin); i++)
y[i] = ya[i]+slowc_save*direction[i];
iter--;
......@@ -2609,11 +2645,11 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
End_GE(Size);
}
else if (stack_solve_algo == 1 || stack_solve_algo == 4)
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc);
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, true, 0);
else if (stack_solve_algo == 2)
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck);
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, true, 0);
else if (stack_solve_algo == 3)
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck);
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, true, 0);
}
if (print_it)
{
......
mex/sources/bytecode/SparseMatrix.hh
View file @ 77b60969
......@@ -103,7 +103,7 @@ class SparseMatrix
{
public:
SparseMatrix();
int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number, int stack_solve_algo);
int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number, int stack_solve_algo, int endo_name_length, char *P_endo_names);
bool simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int Block_number, int solve_algo);
void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter);
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
......@@ -117,9 +117,9 @@ private:
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m);
void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM);
void End_GE(int Size);
void Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l);
void Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block);
void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block);
void Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_);
void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_);
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
#ifdef PROFILER
, long int *ndiv, long int *nsub
......
mex/sources/bytecode/bytecode.cc
View file @ 77b60969
......@@ -16,24 +16,8 @@
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
////////////////////////////////////////////////////////////////////////
// simulate.cc //
// simulate file designed for GNU GCC C++ compiler //
////////////////////////////////////////////////////////////////////////
//#define _GLIBCXX_USE_C99_FENV_TR1 1
//#include <cfenv>
#include <cstring>
#include "Interpreter.hh"
#ifndef DEBUG_EX
# include "mex.h"
#else
# include "mex_interface.hh"
#endif
#include "Mem_Mngr.hh"
#ifdef DEBUG_EX
......@@ -49,168 +33,7 @@ Get_Argument(const char *argv)
}
int
main(int argc, const char *argv[])
{
FILE *fid;
bool steady_state = false;
bool evaluate = false;
bool block = false;
if (argc < 2)
{
mexPrintf("model filename expected\n");
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
float f_tmp;
ostringstream tmp_out("");
tmp_out << argv[1] << "_options.txt";
int nb_params;
int i, row_y, col_y, row_x, col_x;
double *yd, *xd;
double *direction;
int minimal_solving_periods;
string file_name(argv[1]);
int count_array_argument = 0;
for (i = 2; i < argc; i++)
{
if (Get_Argument(argv[i]) == "static")
steady_state = true;
else if (Get_Argument(argv[i]) == "dynamic")
steady_state = false;
else if (Get_Argument(argv[i]) == "evaluate")
evaluate = true;
else if (Get_Argument(argv[i]) == "block")
block = true;
else
{
mexPrintf("Unknown argument : ");
mexEvalString("st=fclose('all');clear all;");
string f;
f = Get_Argument(argv[i]);
f.append("\n");
mexErrMsgTxt(f.c_str());
}
}
fid = fopen(tmp_out.str().c_str(), "r");
int periods = 1;
if (!steady_state)
fscanf(fid, "%d", &periods);
int maxit_;
fscanf(fid, "%d", &maxit_);
fscanf(fid, "%f", &f_tmp);
double slowc = f_tmp;
fscanf(fid, "%f", &f_tmp);
double markowitz_c = f_tmp;
fscanf(fid, "%f", &f_tmp);
double solve_tolf = f_tmp;
fscanf(fid, "%d", &minimal_solving_periods);
fclose(fid);
tmp_out.str("");
tmp_out << argv[1] << "_M.txt";
fid = fopen(tmp_out.str().c_str(), "r");
int y_kmin;
fscanf(fid, "%d", &y_kmin);
int y_kmax;
fscanf(fid, "%d", &y_kmax);
int y_decal;
fscanf(fid, "%d", &y_decal);
fscanf(fid, "%d", &nb_params);
fscanf(fid, "%d", &row_x);
fscanf(fid, "%d", &col_x);
fscanf(fid, "%d", &row_y);
fscanf(fid, "%d", &col_y);
int steady_row_y, steady_col_y;
int steady_row_x, steady_col_x;
fscanf(fid, "%d", &steady_row_y);
fscanf(fid, "%d", &steady_col_y);
fscanf(fid, "%d", &steady_row_x);
fscanf(fid, "%d", &steady_col_x);
int nb_row_xd;
fscanf(fid, "%d", &nb_row_xd);
double *params = (double *) malloc(nb_params*sizeof(params[0]));
for (i = 0; i < nb_params; i++)
{
fscanf(fid, "%f", &f_tmp);
params[i] = f_tmp;
}
fclose(fid);
yd = (double *) malloc(row_y*col_y*sizeof(yd[0]));
xd = (double *) malloc(row_x*col_x*sizeof(xd[0]));
tmp_out.str("");
tmp_out << argv[1] << "_oo.txt";
fid = fopen(tmp_out.str().c_str(), "r");
for (i = 0; i < col_y*row_y; i++)
{
fscanf(fid, "%f", &f_tmp);
yd[i] = f_tmp;
}
for (i = 0; i < col_x*row_x; i++)
{
fscanf(fid, "%f", &f_tmp);
xd[i] = f_tmp;
}
double *steady_yd, *steady_xd;
steady_yd = (double *) malloc(steady_row_y*steady_col_y*sizeof(steady_yd[0]));
steady_xd = (double *) malloc(steady_row_x*steady_col_x*sizeof(steady_xd[0]));
for (i = 0; i < steady_row_y*steady_col_y; i++)
{
fscanf(fid, "%f", &f_tmp);
steady_yd[i] = f_tmp;
}
for (i = 0; i < steady_row_x*steady_col_x; i++)
{
fscanf(fid, "%f", &f_tmp);
steady_xd[i] = f_tmp;
}
fclose(fid);
int size_of_direction = col_y*row_y*sizeof(double);
double *y = (double *) mxMalloc(size_of_direction);
double *ya = (double *) mxMalloc(size_of_direction);
direction = (double *) mxMalloc(size_of_direction);
memset(direction, 0, size_of_direction);
double *x = (double *) mxMalloc(col_x*row_x*sizeof(double));
for (i = 0; i < row_x*col_x; i++)
x[i] = double (xd[i]);
for (i = 0; i < row_y*col_y; i++)
y[i] = double (yd[i]);
free(yd);
free(xd);
int y_size = row_y;
int 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);
string f(file_name);
int nb_blocks = 0;
interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks);
clock_t t1 = clock();
if (!evaluate)
mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
if (x)
mxFree(x);
if (y)
mxFree(y);
if (ya)
mxFree(ya);
if (direction)
mxFree(direction);
if (steady_yd)
mxFree(steady_yd);