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

- Correction of several bugs 

- Reduction of the SparseMatrix.cc code

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2835 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2009-07-10 15:10:11 +00:00
parent d19ccced9a
commit d2dfda62d6
4 changed files with 357 additions and 1223 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

295
mex/sources/simulate/Interpreter.cc
View File

@ -40,9 +40,11 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
solve_tolf=solve_tolf_arg;
size_of_direction=size_of_direction_arg;
slowc=slowc_arg;
slowc_save = slowc;
y_decal=y_decal_arg;
markowitz_c=markowitz_c_arg;
filename=filename_arg;
T=NULL;
//GaussSeidel=true;
}
@ -52,6 +54,7 @@ Interpreter::pow1(double a, double b)
double r=pow_(a,b);
if (isnan(r) || isinf(r))
{
//mexPrintf("pow(%f, %f)=%f\n",a, b, r);
max_res=res1=res2=r;
return(r);
}
@ -65,7 +68,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
{
int var, lag, op;
double v1, v2;
char cc;
char/*uint8_t*/ cc;
bool go_on=true;
double *ll;
while (go_on)
@ -85,7 +88,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
switch (get_code_char)
{
case eParameter :
var=get_code_int
var=get_code_int;
#ifdef DEBUGC
if(Block_Count==2)
{
@ -96,8 +99,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
Stack.push(params[var]);
break;
case eEndogenous :
var=get_code_int
lag=get_code_int
var=get_code_int;
lag=get_code_int;
#ifdef DEBUGC
if(Block_Count==2)
{
@ -110,8 +113,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
#ifdef DEBUGC
mexPrintf("Exogenous\n");
#endif
var=get_code_int
lag=get_code_int
var=get_code_int;
lag=get_code_int;
#ifdef DEBUGC
if(Block_Count==2)
{
@ -125,8 +128,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
mexPrintf("ExogenousDet\n");
#endif
var=get_code_int
lag=get_code_int
var=get_code_int;
lag=get_code_int;
#ifdef DEBUGC
mexPrintf(" x(det)[%d]=%f\n",it_+lag+var*nb_row_xd,x[it_+lag+var*nb_row_xd]);
mexEvalString("drawnow;");
@ -134,12 +137,12 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
Stack.push(x[it_+lag+var*nb_row_xd]);
break;
default:
mexPrintf("Unknown vraibale type\n");
mexPrintf("Unknown variable type\n");
}
break;
case FLDT :
//load a temporary variable in the processor
var=get_code_int
var=get_code_int;
#ifdef DEBUGC
mexPrintf("FLDT %d\n",var);
mexEvalString("drawnow;");
@ -152,7 +155,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
mexPrintf("FLDU\n");
mexEvalString("drawnow;");
#endif
var=get_code_int
var=get_code_int;
var+=Per_u_;
Stack.push(u[var]);
break;
@ -162,7 +165,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
mexPrintf("FLDR\n");
mexEvalString("drawnow;");
#endif
var=get_code_int
var=get_code_int;
Stack.push(r[var]);
break;
case FLDZ :
@ -177,7 +180,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
//load a numerical constant in the processor
/*asm("fld\n\t"
"fstp %%st" : "=t" (ll) : "0" ((double)(*Code)));*/
ll=get_code_pdouble
ll=get_code_pdouble;
#ifdef DEBUGC
mexPrintf("FLDC %f\n",*ll);
mexEvalString("drawnow;");
@ -193,7 +196,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
switch (get_code_char)
{
case eParameter :
var=get_code_int
var=get_code_int;
params[var] = Stack.top();
#ifdef DEBUGC
mexPrintf("FSTP params[%d]=%f\n", var, params[var]);
@ -205,8 +208,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
#ifdef DEBUGC
mexPrintf("FSTP Endogenous\n");
#endif
var=get_code_int
lag=get_code_int
var=get_code_int;
lag=get_code_int;
#ifdef DEBUGC
//mexPrintf("y[%d(it_=%d, lag=%d, y_size=%d, var=%d)](%d)=",(it_+lag)*y_size+var,it_, lag, y_size, var, Stack.size());
mexPrintf("y[it_=%d, lag=%d, y_size=%d, var=%d, block=%d)]=",it_, lag, y_size, var+1, Block_Count+1);
@ -223,9 +226,9 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
#ifdef DEBUGC
mexPrintf("Exogenous\n");
#endif
var=get_code_int
var=get_code_int
lag=get_code_int
var=get_code_int;
var=get_code_int;
lag=get_code_int;
x[it_+lag+var*nb_row_x] = Stack.top();
Stack.pop();
break;
@ -233,10 +236,9 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
#ifdef DEBUGC
mexPrintf("ExogenousDet\n");
#endif
var=get_code_int
var=get_code_int
lag=get_code_int
var=get_code_int;
var=get_code_int;
lag=get_code_int;
x[it_+lag+var*nb_row_xd] = Stack.top();
Stack.pop();
break;
@ -246,7 +248,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
break;
case FSTPT :
//load a temporary variable in the processor
var=get_code_int
var=get_code_int;
#ifdef DEBUGC
mexPrintf("FSTPT T[(var=%d, it_=%d, periods=%d, y_kmin=%d, y_kmax=%d)%d]=", var, it_, periods, y_kmin, y_kmax, var*(periods+y_kmin+y_kmax)+it_);
mexEvalString("drawnow;");
@ -260,7 +262,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
break;
case FSTPU :
//load u variable in the processor
var=get_code_int
var=get_code_int;
var+=Per_u_;
#ifdef DEBUGC
mexPrintf("FSTPU u[%d]",var);
@ -275,7 +277,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
break;
case FSTPR :
//load u variable in the processor
var=get_code_int
var=get_code_int;
r[var] = Stack.top();
#ifdef DEBUGC
mexPrintf("FSTPR residual[%d]=%f\n",var,r[var]);
@ -285,7 +287,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
break;
case FSTPG :
//load u variable in the processor
var=get_code_int
var=get_code_int;
g1[var] = Stack.top();
#ifdef DEBUGC
mexPrintf("FSTPG g1[%d)=%f\n",var,g1[var]);
@ -298,8 +300,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
mexPrintf("FBINARY\n");
mexEvalString("drawnow;");
#endif
op=get_code_int
v2=Stack.top();
op=get_code_int;
v2=Stack.top();
Stack.pop();
v1=Stack.top();
Stack.pop();
@ -376,11 +378,11 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
#endif
break;
case oPower:
Stack.push(pow1(v1, v2));
#ifdef DEBUGC
mexPrintf("pow(%f, %f)\n",v1,v2);
mexEvalString("drawnow;");
#endif
Stack.push(pow1(v1, v2));
break;
case oMax:
Stack.push(max(v1, v2));
@ -407,7 +409,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
mexPrintf("FUNARY\n");
mexEvalString("drawnow;");
#endif
op=get_code_int
op=get_code_int;
v1=Stack.top();
Stack.pop();
switch (op)
@ -503,7 +505,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
mexEvalString("drawnow;");
#endif
break;
case oAcosh:
/* case oAcosh:
Stack.push(acosh(v1));
#ifdef DEBUGC
mexPrintf("acosh\n");
@ -523,7 +525,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
mexPrintf("atanh\n");
mexEvalString("drawnow;");
#endif
break;
break;*/
case oSqrt:
Stack.push(sqrt(v1));
#ifdef DEBUGC
@ -550,8 +552,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
case FENDBLOCK :
//it's the block end
#ifdef DEBUGC
mexPrintf("FENDBLOCK\n");
mexEvalString("drawnow;");
//mexPrintf("FENDBLOCK\n");
//mexEvalString("drawnow;");
#endif
//Block[Block_Count].end=get_code_pos;
go_on=false;
@ -569,7 +571,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
mexPrintf("FOK\n");
mexEvalString("drawnow;");
#endif
op=get_code_int
op=get_code_int;
#ifdef DEBUGC
mexPrintf("var=%d\n",op);
#endif
@ -626,10 +628,12 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
{
case EVALUATE_FORWARD :
//case EVALUATE_FORWARD_R :
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("EVALUATE_FORWARD\n");
#endif
mexEvalString("drawnow;");
//#endif
begining=get_code_pointer;
//mexPrintf("y_kmin=%d periods=%d",y_kmin, periods);
for (it_=y_kmin;it_<periods+y_kmin;it_++)
{
//mexPrintf("begining=%x\n",begining);
@ -637,26 +641,14 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Per_y_=it_*y_size;
//mexPrintf("bef compute_block_time()\n");
compute_block_time(0);
#ifdef PRINT_OUT
for (int j = 34; j</*size*/37; j++)
mexPrintf("y[%d, %d] = %f\n",j+1, it_, y[j+it_*y_size]);
//mexPrintf("y[%d, %d] = %f\n", Block_Contain[j].Variable, it_, y[Per_y_ + Block_Contain[j].Variable]);
#endif
}
/*mexPrintf("Evaluate Forward\n");
for (it_=y_kmin;it_<periods+y_kmin;it_++)
{
mexPrintf("it_=%d ",it_);
for(i=0; i<y_size; i++)
mexPrintf(" %f",y[i+it_*y_size]);
mexPrintf("\n");
}*/
break;
case EVALUATE_BACKWARD :
//case EVALUATE_BACKWARD_R :
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("EVALUATE_BACKWARD\n");
#endif
mexEvalString("drawnow;");
//#endif
begining=get_code_pointer;
for (it_=periods+y_kmin-1;it_>=y_kmin;it_--)
{
@ -670,9 +662,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
}
break;
case SOLVE_FORWARD_SIMPLE :
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("SOLVE_FORWARD_SIMPLE\n");
#endif
mexEvalString("drawnow;");
//#endif
g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer;
@ -707,9 +700,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mxFree(r);
break;
case SOLVE_BACKWARD_SIMPLE :
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("SOLVE_BACKWARD_SIMPLE\n");
#endif
mexEvalString("drawnow;");
//#endif
g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer;
@ -744,85 +738,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mxFree(g1);
mxFree(r);
break;
/*case SOLVE_TWO_BOUNDARIES_SIMPLE :
#ifdef DEBUGC
mexPrintf("SOLVE_TWO_BOUNDARIES_SIMPLE\n");
#endif
is_linear=get_code_bool;
max_lag_plus_max_lead_plus_1=get_code_int;
symbol_table_endo_nbr=get_code_int;
Block_List_Max_Lag=get_code_int;
Block_List_Max_Lead=get_code_int;
Read_file(file_name, periods, max_lag_plus_max_lead_plus_1, symbol_table_endo_nbr, Block_List_Max_Lag, Block_List_Max_Lead, nb_endo, u_count, u_count_init, u);
//sparse_matrix.initialize(periods, nb_endo, y_kmin, y_kmax, y_size, u_count, u_count_init, u, y, ya, slowc, y_decal, markowitz_c, res1, res2, max_res);
g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer;
if (!is_linear)
{
cvg=false;
iter=0;
while (!(cvg||(iter>maxit_)))
{
res2=0;
res1=0;
max_res=0;
for (it_=y_kmin;it_<periods+y_kmin;it_++)
{
Per_u_=(it_-y_kmin)*max_lag_plus_max_lead_plus_1;
set_code_pointer(begining);
Per_y_=it_*y_size;
compute_block_time();
for (i=0; i<size; i++)
{
double rr;
rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]);
if (max_res<fabs(rr))
max_res=fabs(rr);
res2+=rr*rr;
res1+=fabs(rr);
}
}
iter++;
cvg=(max_res<solve_tolf);
Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax,size, periods, true, iter);
}
if (!cvg)
{
mexPrintf("Convergence not achieved in block %d, after %d iterations\n",Block_Count,iter);
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("End of simulate");
}
}
else
{
for (it_=y_kmin;it_<periods+y_kmin;it_++)
{
Per_u_=(it_-y_kmin)*max_lag_plus_max_lead_plus_1;
set_code_pointer(begining);
Per_y_=it_*y_size;
compute_block_time();
for (i=0; i<size; i++)
{
double rr;
rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]);
if (max_res<fabs(rr))
max_res=fabs(rr);
res2+=rr*rr;
res1+=fabs(rr);
}
}
iter++;
cvg=(max_res<solve_tolf);
Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax,size, periods, true, iter);
}
break;*/
case SOLVE_FORWARD_COMPLETE :
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("SOLVE FORWARD_COMPLETE\n");
mexPrintf("confirmation!\n");
mexEvalString("drawnow;");
#endif
//#endif
is_linear=get_code_bool;
//mexPrintf("is_linear=%d\n",is_linear);
//mexEvalString("drawnow;");
@ -838,7 +758,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Block_List_Max_Lead=get_code_int;
//mexPrintf("Block_List_Max_Lead=%d\n",Block_List_Max_Lead);
//mexEvalString("drawnow;");
u_count_int=get_code_int
u_count_int=get_code_int;
//mexPrintf("u_count_int=%d\n",u_count_int);
//mexEvalString("drawnow;");
fixe_u(&u, u_count_int, u_count_int);
@ -863,7 +783,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
{
set_code_pointer(begining);
compute_block_time(0);
/*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK\n");
/*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK it_%d\n",it_);
mexEvalString("drawnow;");*/
//Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter);
simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter);
@ -911,7 +831,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
//Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter);
//simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, /*true*/false, cvg, iter);
cvg=false;
/*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK it_%d\n",it_);
mexEvalString("drawnow;");*/
simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter);
/*mexPrintf("End of simulate_NG\n");
mexEvalString("drawnow;");*/
//simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter);
}
/*mexPrintf("solve forward complete simulation\n");
@ -923,15 +847,19 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mexPrintf("\n");
}*/
}
/*mexPrintf("end of forward solve\n");
mexEvalString("drawnow;");*/
mxFree(index_equa);
memset(direction,0,size_of_direction);
mxFree(g1);
mxFree(r);
mxFree(u);
break;
case SOLVE_BACKWARD_COMPLETE :
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("SOLVE_BACKWARD_COMPLETE\n");
#endif
mexEvalString("drawnow;");
//#endif
is_linear=get_code_bool;
max_lag_plus_max_lead_plus_1=get_code_int;
symbol_table_endo_nbr=get_code_int;
@ -939,7 +867,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Block_List_Max_Lead=get_code_int;
//Read_file(file_name, periods, 0, symbol_table_endo_nbr, Block_List_Max_Lag, Block_List_Max_Lead, nb_endo, u_count, u_count_init, u);
//sparse_matrix.initialize(periods, nb_endo, y_kmin, y_kmax, y_size, u_count, u_count_init, u, y, ya, slowc, y_decal, markowitz_c, res1, res2, max_res);
u_count_int=get_code_int
u_count_int=get_code_int;
//mexPrintf("u_count_int=%d\n",u_count_int);
//mexEvalString("drawnow;");
fixe_u(&u, u_count_int, u_count_int);
@ -1011,36 +939,36 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
#if GNUVER >= 432
//mexPrintf("omp_get_max_threads=%d\n",omp_get_max_threads());
#endif
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("SOLVE_TWO_BOUNDARIES_COMPLETE\n");
mexEvalString("drawnow;");
#endif
//#endif
is_linear=get_code_bool;
#ifdef DEBUGC
mexPrintf("is_linear=%d\n",is_linear);
mexEvalString("drawnow;");
#endif
max_lag_plus_max_lead_plus_1=get_code_int
max_lag_plus_max_lead_plus_1=get_code_int;
#ifdef DEBUGC
mexPrintf("max_lag_plus_max_lead_plus_1=%d\n",max_lag_plus_max_lead_plus_1);
mexEvalString("drawnow;");
#endif
symbol_table_endo_nbr=get_code_int
symbol_table_endo_nbr=get_code_int;
#ifdef DEBUGC
mexPrintf("symbol_table_endo_nbr=%d\n",symbol_table_endo_nbr);
mexEvalString("drawnow;");
#endif
Block_List_Max_Lag=get_code_int
Block_List_Max_Lag=get_code_int;
#ifdef DEBUGC
mexPrintf("Block_List_Max_Lag=%d\n",Block_List_Max_Lag);
mexEvalString("drawnow;");
#endif
Block_List_Max_Lead=get_code_int
Block_List_Max_Lead=get_code_int;
#ifdef DEBUGC
mexPrintf("Block_List_Max_Lead=%d\n",Block_List_Max_Lead);
mexEvalString("drawnow;");
#endif
u_count_int=get_code_int
u_count_int=get_code_int;
#ifdef DEBUGC
mexPrintf("u_count_int=%d\n",u_count_int);
mexPrintf("periods=%d\n",periods);
@ -1113,7 +1041,6 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
if (isnan(res1)||isinf(res1))
{
memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
//GaussSeidel=false;
break;
}
for (i=0; i< size; i++)
@ -1153,6 +1080,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
/*mexPrintf("u_count=%d &u_count=%x\n",u_count,&u_count);
mexEvalString("drawnow;");*/
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter);
/*mexPrintf("after simulate_NG1\n");
mexEvalString("drawnow;");*/
/*mexPrintf("aft simulate_NG1\n");
mexEvalString("drawnow;");*/
}
@ -1195,7 +1124,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
res1=res2=max_res=0;max_res_idx=0;
cvg = false;
if(Gaussian_Elimination)
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter);
{
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter);
/*mexPrintf("after simulate_NG1\n");
mexEvalString("drawnow;");*/
}
else
{
#ifdef LINBCG
@ -1217,11 +1150,16 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
#endif
//mxFree(g1);
/*mexPrintf("end of simulate_a_block\n");
mexEvalString("drawnow;");*/
mxFree(r);
mxFree(y_save);
mxFree(u);
mxFree(index_vara);
mxFree(index_equa);
memset(direction,0,size_of_direction);
/*mexPrintf("after free\n");
mexEvalString("drawnow;");*/
//GaussSeidel=false;
break;
default:
@ -1230,6 +1168,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mexEvalString("drawnow;");
mexErrMsgTxt("End of simulate");
}
/*mexPrintf("finish simulate_a_block\n");
mexEvalString("drawnow;");*/
}
void
@ -1237,6 +1177,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename)
{
ifstream CompiledCode;
int Code_Size, var;
//printf("open(%s)\n",(file_name + ".cod").c_str());
//First read and store inn memory the code
CompiledCode.open((file_name + ".cod").c_str(),std::ios::in | std::ios::binary| std::ios::ate);
if (!CompiledCode.is_open())
@ -1279,77 +1220,89 @@ Interpreter::compute_blocks(string file_name, string bin_basename)
{
case FBEGINBLOCK :
//it's a new block
Block_Count++;
Block_type lBlock;
Block.clear();
Block_Contain.clear();
Block_contain_type lBlock_Contain;
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("FBEGINBLOCK\n");
mexEvalString("drawnow;");
#endif
lBlock.begin=get_code_pos-(long int)(Init_Code);
#ifdef DEBUGC
//#endif
/*uint64_t *Init_Code_64;
memcpy(Init_Code_64, Init_Code, sizeof(*Init_Code_64));*/
lBlock.begin=get_code_pos-(uint64_t*)Init_Code;
//#ifdef DEBUGC
mexPrintf("Block[%d].begin=%d\n",Block_Count, lBlock.begin);
mexEvalString("drawnow;");
#endif
lBlock.size=get_code_int
#ifdef DEBUGC
//#endif
lBlock.size=get_code_int;
//#ifdef DEBUGC
mexPrintf("Block[Block_Count].size=%d\n",lBlock.size);
mexEvalString("drawnow;");
#endif
lBlock.type=get_code_int
#ifdef DEBUGC
//#endif
lBlock.type=get_code_int;
//#ifdef DEBUGC
mexPrintf("Block[Block_Count].type=%d\n",lBlock.type);
mexEvalString("drawnow;");
#endif
//#endif
Block.push_back(lBlock);
for (int i=0;i</*Block[Block_Count].size*/lBlock.size;i++)
{
lBlock_Contain.Variable=get_code_int
#ifdef DEBUGC
lBlock_Contain.Variable=get_code_int;
//#ifdef DEBUGC
mexPrintf("Block_Contain[%d].Variable=%d\n",i,lBlock_Contain.Variable);
mexEvalString("drawnow;");
#endif
lBlock_Contain.Equation=get_code_int
#ifdef DEBUGC
//#endif
lBlock_Contain.Equation=get_code_int;
//#ifdef DEBUGC
mexPrintf("Block_Contain[%d].Equation=%d\n",i,lBlock_Contain.Equation);
mexEvalString("drawnow;");
#endif
lBlock_Contain.Own_Derivative=get_code_int
//#endif
lBlock_Contain.Own_Derivative=get_code_int;
//mexPrintf("Block_Contain[%d].Own_Derivative=%d\n",i,lBlock_Contain.Own_Derivative);
Block_Contain.push_back(lBlock_Contain);
}
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("Block Completed\n");
mexEvalString("drawnow;");
#endif
//#endif
simulate_a_block(lBlock.size,lBlock.type, file_name, bin_basename,true);
/*mexPrintf("after simulate_a_block\n");
mexEvalString("drawnow;");*/
//simulate_a_block(lBlock.size,lBlock.type, file_name, bin_basename,false);
break;
case FEND :
#ifdef DEBUGC
//#ifdef DEBUGC
mexPrintf("FEND\n");
mexEvalString("drawnow;");
#endif
//#endif
go_on=false;
break;
case FDIMT :
var=get_code_int
var=get_code_int;
#ifdef DEBUGC
mexPrintf("FDIMT var=%d mxMalloc(%d)\n",var,var*(periods+y_kmin+y_kmax)*sizeof(double));
mexEvalString("drawnow;");
#endif
if(T)
mxFree(T);
T=(double*)mxMalloc(var*(periods+y_kmin+y_kmax)*sizeof(double));
break;
default :
mexPrintf("Unknow command : %d at pos %d !!\n",(long int)(code),(long int)(get_code_pos)-(long int)(Init_Code));
mexPrintf("Unknow command : %d at pos %d !!\n",(long int)(code),(uint64_t*)(get_code_pos)-(uint64_t*)(Init_Code));
mexEvalString("st=fclose('all');clear all;");
mexEvalString("drawnow;");
mexErrMsgTxt("End of simulate");
break;
}
}
mxFree(Init_Code);
if(T)
mxFree(T);
/*mexPrintf("compute_blocks\n");
mexEvalString("drawnow;");*/
}

9
mex/sources/simulate/Mem_Mngr.cc
View File

@ -16,7 +16,7 @@
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdint.h>
#include "Mem_Mngr.hh"
Mem_Mngr::Mem_Mngr()
@ -82,6 +82,7 @@ Mem_Mngr::mxMalloc_NZE()
mexPrintf("CHUNK_BLCK_SIZE=%d\n",CHUNK_BLCK_SIZE);
#endif
NZE_Mem=(NonZeroElem*)mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem));
//mexPrintf("in mxMalloc NZE_Mem=%x CHUNK_heap_pos=%d CHUNK_BLCK_SIZE=%d Nb_CHUNK=%d\n",NZE_Mem, CHUNK_heap_pos, CHUNK_BLCK_SIZE, Nb_CHUNK);
if(!NZE_Mem)
{
mexPrintf("Not enough memory available\n");
@ -123,7 +124,7 @@ Mem_Mngr::mxFree_NZE(void* pos)
mexPrintf("NZE_Mem_add[i*CHUNK_BLCK_SIZE]=%d\n",NZE_Mem_add[i*CHUNK_BLCK_SIZE]);
mexEvalString("drawnow;");
}*/
gap=((long int)(pos)-(long int)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem);
gap=((uint64_t)(pos)-(uint64_t)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem);
if ((gap<CHUNK_BLCK_SIZE) && (gap>=0))
break;
}
@ -277,8 +278,12 @@ void
Mem_Mngr::Free_All()
{
int i;
/*mexPrintf("Nb_CHUNK=%d\n",Nb_CHUNK);
mexEvalString("drawnow;");*/
for (i=0;i<Nb_CHUNK;i++)
{
/*mexPrintf("NZE_Mem_add[%d]=%x\n",i*CHUNK_BLCK_SIZE,NZE_Mem_add[i*CHUNK_BLCK_SIZE]);
mexEvalString("drawnow;");*/
mxFree(NZE_Mem_add[i*CHUNK_BLCK_SIZE]);
}
mxFree(NZE_Mem_add);

1069
mex/sources/simulate/SparseMatrix.cc
View File

File diff suppressed because it is too large Load Diff

207
mex/sources/simulate/simulate.cc
View File

@ -42,28 +42,157 @@ max(int a, int b)
#ifdef DEBUG_EX
/*The Matlab c++ interface*/
using namespace std;
#include <sstream>
#include "mex_interface.hh"
/*class EvalException
int
main( int argc, const char* argv[] )
{
};*/
FILE *fid;
printf("argc=%d\n",argc);
float f_tmp;
ostringstream tmp_out("");
tmp_out << argv[1] << "_options.txt";
cout << tmp_out.str().c_str() << "\n";
int nb_params;
int i, row_y, col_y, row_x, col_x;
double *yd, *xd;
double *direction;
string file_name(argv[1]);
//mexPrintf("file_name=%s\n",file_name.c_str());
fid = fopen(tmp_out.str().c_str(),"r");
fscanf(fid,"%d",&periods);
fscanf(fid,"%d",&maxit_);
fscanf(fid,"%f",&f_tmp);
slowc = f_tmp;
//mexPrintf("slowc_save=%f\n",slowc_save);
fscanf(fid,"%f",&f_tmp);
markowitz_c = f_tmp;
fscanf(fid,"%f",&f_tmp);
solve_tolf = f_tmp;
fclose(fid);
tmp_out.str("");
tmp_out << argv[1] << "_M.txt";
//printf("%s\n",tmp_out.str().c_str());
fid = fopen(tmp_out.str().c_str(),"r");
fscanf(fid,"%d",&y_kmin);
//printf("y_kmin=%d\n",y_kmin);
fscanf(fid,"%d",&y_kmax);
//printf("y_kmax=%d\n",y_kmax);
fscanf(fid,"%d",&y_decal);
//printf("y_decal=%d\n",y_decal);
fscanf(fid,"%d",&nb_params);
//printf("nb_params=%d\n",nb_params);
fscanf(fid,"%d",&row_x);
//printf("row_x=%d\n",row_x);
fscanf(fid,"%d",&col_x);
//printf("col_x=%d\n",col_x);
fscanf(fid,"%d",&row_y);
//printf("row_y=%d\n",row_y);
fscanf(fid,"%d",&col_y);
//printf("col_y=%d\n",col_y);
fscanf(fid,"%d",&nb_row_xd);
//printf("nb_row_xd=%d\n",nb_row_xd);
params = (double*)malloc(nb_params*sizeof(params[0]));
//printf("OK1\n");
for(i=0; i < nb_params; i++)
{
fscanf(fid,"%f",&f_tmp);
params[i] = f_tmp;
//printf("param[%d]=%f\n",i,params[i]);
}
//printf("OK2\n");
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;
}
fclose(fid);
size_of_direction=col_y*row_y*sizeof(double);
y=(double*)mxMalloc(size_of_direction);
ya=(double*)mxMalloc(size_of_direction);
direction=(double*)mxMalloc(size_of_direction);
memset(direction,0,size_of_direction);
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);
y_size=row_y;
x_size=col_x/*row_x*/;
nb_row_x=row_x;
/*for(int i=0; i<y_size; i++)
{
for(int it_=0; it_<8;it_++)
mexPrintf("y[t=%d, var=%d]=%f ",it_+1, i+1, y[(it_)*y_size+i]);
mexPrintf("\n");
}
for(int i=0; i<col_x; i++)
{
for(int it_=0; it_<8;it_++)
mexPrintf("x[t=%d, var=%d]=%f ",it_, i+1, x[it_+i*nb_row_x]);
mexPrintf("\n");
}*/
t0= clock();
Interpreter interprete(params, y, ya, x, 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);
string f(file_name);
interprete.compute_blocks(f+"_dynamic", f);
t1= clock();
mexPrintf("Simulation Time=%f milliseconds\n",1000.0*(double(t1)-double(t0))/double(CLOCKS_PER_SEC));
/*if (nlhs>0)
{
plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
pind = mxGetPr(plhs[0]);
for (i=0;i<row_y*col_y;i++)
pind[i]=y[i];
}*/
if(x)
mxFree(x);
if(y)
mxFree(y);
if(ya)
mxFree(ya);
if(direction)
mxFree(direction);
//#define DEBUG
/* The gateway routine */
free(params);
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
}
#else
/* The gateway routine */
void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mxArray *M_, *oo_, *options_;
int i, row_y, col_y, row_x, col_x;
double * pind ;
double *direction;
//mexPrintf("mexFunction\n");
//mexEvalString("drawnow;");
/* Gets model parameters from global workspace of Matlab */
//mexPrintf("starting simulation\n");
M_ = mexGetVariable("global","M_");
if (M_ == NULL )
{
@ -100,7 +229,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
periods=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"periods"))))));
maxit_=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"maxit_"))))));
slowc=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"slowc")))));
slowc_save=slowc;
//slowc_save=slowc;
markowitz_c=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"markowitz")))));
nb_row_xd=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"exo_det_nbr"))))));
mxArray *mxa=mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"fname"));
@ -111,15 +240,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int status = mxGetString(mxa, fname, buflen);
if (status != 0)
mexWarnMsgTxt("Not enough space. Filename is truncated.");
//mexPrintf("fname=%s\n",fname);
col_y=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"endo_simul")));;
/*if (col_y<row_x)
{
row_y=row_y/row_x;
col_y=row_x;
}*/
solve_tolf=*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"dynatol"))));
//mexPrintf("col_y=%d row_y=%d\n",col_y, row_y);
size_of_direction=col_y*row_y*sizeof(double);
y=(double*)mxMalloc(size_of_direction);
ya=(double*)mxMalloc(size_of_direction);
@ -128,46 +250,32 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
x=(double*)mxMalloc(col_x*row_x*sizeof(double));
for (i=0;i<row_x*col_x;i++)
x[i]=double(xd[i]);
/*mexPrintf("1 ");
for (i=0;i<row_x*col_x;i++)
{
x[i]=double(xd[i]);
if(!(i % row_x) && i>0)
mexPrintf("\n%d %f ",i/row_x+1,x[i]);
else
mexPrintf("%f ",x[i]);
}
for(i=0;i<row_x;i++)
{
for(j=0;j<col_x;j++)
{
x[j*row_x+i]=xd[i*col_x+j];
mexPrintf("%f ",xd[i*col_x+j]);
}
mexPrintf("\n");
}
*/
for (i=0;i<row_y*col_y;i++)
y[i]=double(yd[i]);
y_size=row_y;
x_size=row_x;
x_size=col_x/*row_x*/;
nb_row_x=row_x;
/* Call the C subroutines. */
//mexPrintf("Call subroutines\n");
//mexEvalString("drawnow;");
//t0= pctimer();
/*for(int i=0; i<y_size; i++)
{
for(int it_=0; it_<8;it_++)
mexPrintf("y[t=%d, var=%d]=%f ",it_+1, i+1, y[(it_)*y_size+i]);
mexPrintf("\n");
}
for(int i=0; i<col_x; i++)
{
for(int it_=0; it_<8;it_++)
mexPrintf("x[t=%d, var=%d]=%f ",it_, i+1, x[it_+i*nb_row_x]);
mexPrintf("\n");
}*/
t0= clock();
Interpreter interprete(params, y, ya, x, 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);
string f(fname);
interprete.compute_blocks(f+"_dynamic", f);
//t1= pctimer();
t1= clock();
mexPrintf("Simulation Time=%f milliseconds\n",1000.0*(double(t1)-double(t0))/double(CLOCKS_PER_SEC));
//mexPrintf("SaveCode.is_open()=%d nlhs=%d \n",SaveCode.is_open(),nlhs);
//interprete.sparse_matrix.close_SaveCode();
//mexPrintf("End all-1\n");
if (nlhs>0)
{
plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
@ -175,7 +283,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
for (i=0;i<row_y*col_y;i++)
pind[i]=y[i];
}
//mexPrintf("End all0\n");
if(x)
mxFree(x);
if(y)
@ -184,5 +291,5 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mxFree(ya);
if(direction)
mxFree(direction);
//mexPrintf("End all\n");
}
#endif