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

- SymbolGaussElim.cc is discarded from the preprocessor (all simulations are now implemented in simulate) 

- Correction in model_info.m

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2358 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2009-01-20 22:04:37 +00:00
parent b1e3185592
commit 4fe7e3e155
17 changed files with 707 additions and 2023 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
matlab/model_info.m
View File

@ -72,17 +72,20 @@ function model_info;
function ret=Sym_type(type);
EVALUATE_FORWARD=0;
EVALUATE_BACKWARD=1;
SOLVE_FORWARD_SIMPLE=2;
SOLVE_BACKWARD_SIMPLE=3;
SOLVE_TWO_BOUNDARIES_SIMPLE=4;
SOLVE_FORWARD_COMPLETE=5;
SOLVE_BACKWARD_COMPLETE=6;
SOLVE_TWO_BOUNDARIES_COMPLETE=7;
EVALUATE_FORWARD_R=8;
EVALUATE_BACKWARD_R=9;
UNKNOWN=0;
EVALUATE_FORWARD=1;
EVALUATE_BACKWARD=2;
SOLVE_FORWARD_SIMPLE=3;
SOLVE_BACKWARD_SIMPLE=4;
SOLVE_TWO_BOUNDARIES_SIMPLE=5;
SOLVE_FORWARD_COMPLETE=6;
SOLVE_BACKWARD_COMPLETE=7;
SOLVE_TWO_BOUNDARIES_COMPLETE=8;
EVALUATE_FORWARD_R=9;
EVALUATE_BACKWARD_R=10;
switch (type)
case (UNKNOWN),
ret='UNKNOWN ';
case {EVALUATE_FORWARD,EVALUATE_FORWARD_R},
ret='EVALUATE FORWARD ';
case {EVALUATE_BACKWARD,EVALUATE_BACKWARD_R},

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

@ -62,7 +62,7 @@ Interpreter::pow1(double a, double b)
void
Interpreter::compute_block_time() /*throw(EvalException)*/
Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
{
int var, lag, op;
double v1, v2;
@ -586,7 +586,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
#endif
//SparseMatrix sparse_matrix;
int nb_endo, u_count_init;
//int nb_endo, u_count_init;
//mexPrintf("simulate_a_block\n");
@ -608,12 +608,20 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
{
set_code_pointer(begining);
Per_y_=it_*y_size;
compute_block_time();
compute_block_time(0);
#ifdef PRINT_OUT
for (j = 0; j<size; j++)
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 :
@ -625,7 +633,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
{
set_code_pointer(begining);
Per_y_=it_*y_size;
compute_block_time();
compute_block_time(0);
#ifdef PRINT_OUT
for (j = 0; j<size; j++)
mexPrintf("y[%d, %d] = %f\n", Block_Contain[j].Variable, it_, y[Per_y_ + Block_Contain[j].Variable]);
@ -648,7 +656,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
{
set_code_pointer(begining);
Per_y_=it_*y_size;
compute_block_time();
compute_block_time(0);
y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0];
//mexPrintf("y[%d] += -r[0] (%f) / g1[0] (%f) = %f\n",Per_y_+Block_Contain[0].Variable,r[0],g1[0],y[Per_y_+Block_Contain[0].Variable]);
double rr;
@ -685,9 +693,9 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
{
set_code_pointer(begining);
Per_y_=it_*y_size;
compute_block_time();
mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_SIMPLE : OK\n");
mexEvalString("drawnow;");
compute_block_time(0);
/*mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_SIMPLE : OK\n");
mexEvalString("drawnow;");*/
y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0];
double rr;
rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]);
@ -781,19 +789,39 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
}
break;*/
case SOLVE_FORWARD_COMPLETE :
//#ifdef DEBUGC
mexPrintf("SOLVE_FORWARD_COMPLETE\n");
//#endif
#ifdef DEBUGC
mexPrintf("SOLVE FORWARD_COMPLETE\n");
mexPrintf("confirmation!\n");
mexEvalString("drawnow;");
#endif
is_linear=get_code_bool;
//mexPrintf("is_linear=%d\n",is_linear);
//mexEvalString("drawnow;");
max_lag_plus_max_lead_plus_1=get_code_int;
//mexPrintf("max_lag_plus_max_lead_plus_1=%d\n",max_lag_plus_max_lead_plus_1);
//mexEvalString("drawnow;");
symbol_table_endo_nbr=get_code_int;
//mexPrintf("symbol_table_endo_nbr=%d\n",symbol_table_endo_nbr);
//mexEvalString("drawnow;");
Block_List_Max_Lag=get_code_int;
//mexPrintf("Block_List_Max_Lag=%d\n",Block_List_Max_Lag);
//mexEvalString("drawnow;");
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
//mexPrintf("u_count_int=%d\n",u_count_int);
//mexEvalString("drawnow;");
fixe_u(&u, u_count_int, u_count_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);
Read_SparseMatrix(bin_basename, size, 1, 0, 0);
//mexPrintf("size=%d\n",size);
g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer;
Per_u_ = 0;
if (!is_linear)
{
for (it_=y_kmin;it_<periods+y_kmin;it_++)
@ -804,10 +832,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
while (!(cvg||(iter>maxit_)))
{
set_code_pointer(begining);
compute_block_time();
mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK\n");
mexEvalString("drawnow;");
Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter);
compute_block_time(0);
/*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK\n");
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);
res2=0;
res1=0;
max_res=0;
@ -841,15 +870,26 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
res1=res2=max_res=0;
/*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE before compute_block_time OK\n");
mexEvalString("drawnow;");*/
compute_block_time();
compute_block_time(0);
//mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK\n");
/*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : %d OK\n",it_);
mexEvalString("drawnow;");*/
//Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter);
//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);
//simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, /*true*/false, cvg, iter);
simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter);
//simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter);
}
/*mexPrintf("solve forward complete simulation\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");
}*/
}
memset(direction,0,size_of_direction);
mxFree(g1);
mxFree(r);
mxFree(u);
@ -863,8 +903,13 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
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, 0, symbol_table_endo_nbr, Block_List_Max_Lag, Block_List_Max_Lead, nb_endo, u_count, u_count_init, u);
//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
//mexPrintf("u_count_int=%d\n",u_count_int);
//mexEvalString("drawnow;");
fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, 1, 0, 0);
g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer;
@ -878,10 +923,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
while (!(cvg||(iter>maxit_)))
{
set_code_pointer(begining);
compute_block_time();
mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_COMPLETE : OK\n");
mexEvalString("drawnow;");
Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter);
compute_block_time(0);
/*mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_COMPLETE : OK\n");
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);
res2=0;
res1=0;
max_res=0;
@ -911,16 +957,19 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
{
set_code_pointer(begining);
Per_y_=it_*y_size;
compute_block_time();
Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 1, false, iter);
compute_block_time(0);
//Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 1, false, iter);
simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter);
}
}
memset(direction,0,size_of_direction);
mxFree(g1);
mxFree(r);
mxFree(u);
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE :
case SOLVE_TWO_BOUNDARIES_COMPLETE:
mexPrintf("omp_get_max_threads=%d\n",omp_get_max_threads());
#ifdef DEBUGC
mexPrintf("SOLVE_TWO_BOUNDARIES_COMPLETE\n");
mexEvalString("drawnow;");
@ -982,7 +1031,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Per_u_=0;
Per_y_=it_*y_size;
set_code_pointer(begining);
compute_block_time();
compute_block_time(0);
linbcg.Initialize(filename, res1, res2, max_res, slowc, ya, direction, iter);
linbcg.Preconditioner(periods, y_kmin, y_kmax, size, IM_i, index_vara, index_equa, y_size, y, true, 0, a, indx);
#endif
@ -1012,7 +1061,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
//mexPrintf("ok\n");
//mexPrintf("compute_block_time\n");
set_code_pointer(begining);
compute_block_time();
compute_block_time(Per_u_);
//mexPrintf("end of compute_block_time\n");
/*if(Gaussian_Elimination)
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);*/
@ -1084,7 +1133,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Per_u_=(it_-y_kmin)*max_lag_plus_max_lead_plus_1;
Per_y_=it_*y_size;
set_code_pointer(begining);
compute_block_time();
compute_block_time(Per_u_);
#ifdef PRINT_OUT
for (j=0; j<max_lag_plus_max_lead_plus_1; j++)
{
@ -1092,6 +1141,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
}
mexPrintf("\n");
#endif
/*mexPrintf("it_=%d ",it_);
for(i=0; i<y_size; i++)
mexPrintf(" %f",y[i]);
mexPrintf("\n");*/
}
res1=res2=max_res=0;
if(Gaussian_Elimination)
@ -1103,6 +1156,14 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
linbcg.SolveLinear(periods, y_kmin, y_kmax, size, IM_i, index_vara, index_equa, y_size, y, true, cvg, a, indx);
#endif
}
/*mexPrintf("Two boundaries simulation\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");
}*/
}
#ifdef DEBUGC
//mexErrMsgTxt("End of simulate");

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

@ -63,7 +63,7 @@ class Interpreter : SparseMatrix
{
protected :
double pow1(double a, double b);
void compute_block_time();
void compute_block_time(int Per_u_);
void simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination);
double *T;
vector<Block_contain_type> Block_Contain;

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

@ -277,7 +277,7 @@ void
Mem_Mngr::Free_All()
{
int i;
for (int i=0;i<Nb_CHUNK;i++)
for (i=0;i<Nb_CHUNK;i++)
{
mxFree(NZE_Mem_add[i*CHUNK_BLCK_SIZE]);
}

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

@ -410,8 +410,8 @@ void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, in
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j));
/*mexPrintf("0i=%d eq=%d var=%d lag=%d j=%d\n",i,eq,var,lag,j );
mexEvalString("drawnow;");*/
//mexPrintf("i=%d eq=%d var=%d lag=%d j=%d\n",i,eq,var,lag,j );
//mexEvalString("drawnow;");
IM_i[std::make_pair(std::make_pair(eq, var), lag)] = j;
/*mexPrintf("1i=%d\n",i);
mexEvalString("drawnow;");*/
@ -419,8 +419,6 @@ void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, in
#ifdef MEM_ALLOC_CHK
mexPrintf("index_vara=(int*)mxMalloc(%d*sizeof(int))\n",Size*(periods+y_kmin+y_kmax));
#endif
/*mexPrintf("Size=%d periods=%d y_kmin=%d y_kmax=%d\n",Size,periods, y_kmin, y_kmax);
mexEvalString("drawnow;");*/
index_vara=(int*)mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int));
#ifdef MEM_ALLOC_CHK
mexPrintf("ok\n");
@ -430,18 +428,21 @@ void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, in
SaveCode.read(reinterpret_cast<char *>(&index_vara[j]), sizeof(*index_vara));
//mexPrintf("index_vara[%d]=%d\n",j,index_vara[j]);
}
for (i=1;i<periods+y_kmin+y_kmax;i++)
if(periods+y_kmin+y_kmax>1)
{
for (j=0;j<Size;j++)
for (i=1;i<periods+y_kmin+y_kmax;i++)
{
#ifdef PRINT_OUT
mexPrintf("index_vara[%d]=index_vara[%d]+y_size=",j+Size*i,j+Size*(i-1));
#endif
index_vara[j+Size*i]=index_vara[j+Size*(i-1)]+y_size;
#ifdef PRINT_OUT
mexPrintf("%d\n",index_vara[j+Size*i]);
#endif
}
for (j=0;j<Size;j++)
{
#ifdef PRINT_OUT
mexPrintf("index_vara[%d]=index_vara[%d]+y_size=",j+Size*i,j+Size*(i-1));
#endif
index_vara[j+Size*i]=index_vara[j+Size*(i-1)]+y_size;
#ifdef PRINT_OUT
mexPrintf("%d\n",index_vara[j+Size*i]);
#endif
}
}
}
index_equa=(int*)mxMalloc(Size*sizeof(int));
for(j=0;j<Size;j++)
@ -452,13 +453,109 @@ void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, in
}
void SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int> ,int>, int> IM)
{
int i, eq, var, lag;
//double tmp_b=0.0;
std::map<std::pair<std::pair<int, int> ,int>, int>::iterator it4;
NonZeroElem* first;
//mexPrintf("periods=%d, y_kmin=%d, y_kmax=%d, SizeInit=%d, IM.size()=%d\n",periods, y_kmin, y_kmax, Size, IM.size());
pivot=(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));
line_done=(bool*)mxMalloc(Size*sizeof(bool));
memset(line_done, 0, Size*sizeof(*line_done));
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);
memset(FNZE_R, 0, i);
memset(FNZE_C, 0, i);
NonZeroElem** temp_NZE_R=(NonZeroElem**)mxMalloc(i);
NonZeroElem** temp_NZE_C=(NonZeroElem**)mxMalloc(i);
memset(temp_NZE_R, 0, i);
memset(temp_NZE_C, 0, i);
i=Size*sizeof(int);
NbNZRow=(int*)mxMalloc(i);
NbNZCol=(int*)mxMalloc(i);
memset(NbNZRow, 0, i);
memset(NbNZCol, 0, i);
i=Size*sizeof(*b);
memset(b,0,i);
it4=IM.begin();
eq=-1;
double tmp_b[Size];
for(i=0; i< Size;i++);
tmp_b[i]=0;//u[i];
int u_count1=Size;
while (it4!=IM.end())
{
var=it4->first.first.second;
/*if (eq!=it4->first.first.first)
tmp_b=0;*/
eq=it4->first.first.first;
lag=it4->first.second;
if (lag==0) /*Build the index for sparse matrix containing the jacobian : u*/
{
//mexPrintf("Add eq=%d, var=%d, lag=%d at it_=%d u=%f\n",eq,var,lag, it_, u[u_count1]);
//mexPrintf(" u_index=%d\n",/*it4->second+u_count_init*it_*/u_count1);
NbNZRow[eq]++;
NbNZCol[var]++;
#ifdef NEW_ALLOC
first=mem_mngr.mxMalloc_NZE();
#else
first=(NonZeroElem*)mxMalloc(sizeof(NonZeroElem));
#endif
first->NZE_C_N=NULL;
first->NZE_R_N=NULL;
first->u_index=u_count1/*it4->second+u_count_init*it_*/;
first->r_index=eq;
first->c_index=var;
first->lag_index=lag;
//mexPrintf(" u[%d](%f)*y[%d](%f)=%f\n",u_count1, u[u_count1], index_vara[var]+it_*y_size, y[index_vara[var]+it_*y_size], u[u_count1]*y[index_vara[var]+it_*y_size]);
tmp_b[eq] += u[u_count1]*y[index_vara[var]+it_*y_size];
if (FNZE_R[eq]==NULL)
{
FNZE_R[eq]=first;
}
if (FNZE_C[var]==NULL)
FNZE_C[var]=first;
if (temp_NZE_R[eq]!=NULL)
temp_NZE_R[eq]->NZE_R_N=first;
if (temp_NZE_C[var]!=NULL)
temp_NZE_C[var]->NZE_C_N=first;
temp_NZE_R[eq]=first;
temp_NZE_C[var]=first;
u_count1++;
}
it4++;
}
for(i=0;i<Size;i++)
{
b[i]=u_count1+i;
u[b[i]]=-tmp_b[i];
//mexPrintf("b[%d]=%f\n",i,u[b[i]]);
}
//mexEvalString("Init");
mxFree(temp_NZE_R);
mxFree(temp_NZE_C);
/*mexPrintf("end of Simple_Init\n");
mexEvalString("drawnow;");*/
}
void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int> ,int>, int> IM)
{
int t,i, eq, var, lag;
double tmp_b=0.0;
std::map<std::pair<std::pair<int, int> ,int>, int>::iterator it4;
NonZeroElem* first;
//mexPrintf("periods=%d, y_kmin=%d, y_kmax=%d, Size=%d, IM.size()=%d\n",periods, y_kmin, y_kmax, Size, IM.size());
//mexPrintf("periods=%d, y_kmin=%d, y_kmax=%d, SizeInit=%d, IM.size()=%d\n",periods, y_kmin, y_kmax, Size, IM.size());
#ifdef MEM_ALLOC_CHK
mexPrintf("pivot=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods);
#endif
@ -1319,73 +1416,6 @@ SparseMatrix::complete(int beg_t, int Size, int periods, int *b)
return(beg_t);
}
double
SparseMatrix::bksub( int tbreak, int last_period, int Size, double slowc_l
#ifdef PROFILER
, /*NonZeroElem *first,*/ long int *nmul
#endif
)
{
NonZeroElem *first;
int i, j, k;
double yy;
res1 = res2 = max_res = 0;
for (i=0;i<y_size*(periods+y_kmin);i++)
y[i]=ya[i];
if (symbolic && tbreak)
last_period=complete(tbreak, Size, periods, b);
else
last_period=periods;
for (int t=last_period+y_kmin-1;t>=y_kmin;t--)
{
int ti=(t-y_kmin)*Size;
#ifdef PRINT_OUT
mexPrintf("t=%d ti=%d\n",t,ti);
#endif
int cal=y_kmin*Size;
int cal_y=y_size*y_kmin;
for (i=ti-1;i>=ti-Size;i--)
{
j=i+cal;
#ifdef PRINT_OUT_y
mexPrintf("t=%d, ti=%d j=%d i+Size=%d\n",t,ti,j,i+Size);
#endif
int pos=pivot[/*j*/i+Size];
#ifdef PRINT_OUT_y
mexPrintf("i-ti+Size=%d pos=%d j=%d\n",i-ti+Size,pos,j);
#endif
int nb_var=At_Row(pos,&first);
first=first->NZE_R_N;
nb_var--;
int eq=index_vara[j]+y_size;
#ifdef PRINT_OUT_y1
mexPrintf("y[index_vara[%d]=%d]=",j,index_vara[j]+y_size);
#endif
yy=0;
for (k=0;k<nb_var;k++)
{
yy+=y[index_vara[first->c_index]+cal_y]*u[first->u_index];
#ifdef PROFILER
(*nmul)++;
#endif
first=first->NZE_R_N;
}
#ifdef PRINT_OUT_y1
mexPrintf("|u[%d](%f)|",b[pos],double(u[b[pos]]));
#endif
yy=-(yy+y[eq]+u[b[pos]]);
direction[eq]=yy;
//mexPrintf("direction[%d] = %f\n",eq,yy);
y[eq] += slowc_l*yy;
#ifdef PRINT_OUT_y1
mexPrintf("=%f (%f)\n",double(yy),double(y[eq]));
#endif
}
}
return res1;
}
void
SparseMatrix::close_swp_file()
{
@ -1955,6 +1985,384 @@ SparseMatrix::Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_k
double
SparseMatrix::bksub( int tbreak, int last_period, int Size, double slowc_l
#ifdef PROFILER
, /*NonZeroElem *first,*/ long int *nmul
#endif
)
{
NonZeroElem *first;
int i, j, k;
double yy;
res1 = res2 = max_res = 0;
for (i=0;i<y_size*(periods+y_kmin);i++)
y[i]=ya[i];
if (symbolic && tbreak)
last_period=complete(tbreak, Size, periods, b);
else
last_period=periods;
for (int t=last_period+y_kmin-1;t>=y_kmin;t--)
{
int ti=(t-y_kmin)*Size;
#ifdef PRINT_OUT
mexPrintf("t=%d ti=%d\n",t,ti);
#endif
int cal=y_kmin*Size;
int cal_y=y_size*y_kmin;
for (i=ti-1;i>=ti-Size;i--)
{
j=i+cal;
#ifdef PRINT_OUT_y
mexPrintf("t=%d, ti=%d j=%d i+Size=%d\n",t,ti,j,i+Size);
#endif
int pos=pivot[/*j*/i+Size];
#ifdef PRINT_OUT_y
mexPrintf("i-ti+Size=%d pos=%d j=%d\n",i-ti+Size,pos,j);
#endif
int nb_var=At_Row(pos,&first);
first=first->NZE_R_N;
nb_var--;
int eq=index_vara[j]+y_size;
#ifdef PRINT_OUT_y1
mexPrintf("y[index_vara[%d]=%d]=",j,index_vara[j]+y_size);
#endif
yy=0;
for (k=0;k<nb_var;k++)
{
yy+=y[index_vara[first->c_index]+cal_y]*u[first->u_index];
#ifdef PROFILER
(*nmul)++;
#endif
first=first->NZE_R_N;
}
#ifdef PRINT_OUT_y1
mexPrintf("|u[%d](%f)|",b[pos],double(u[b[pos]]));
#endif
yy=-(yy+y[eq]+u[b[pos]]);
direction[eq]=yy;
//mexPrintf("direction[%d] = %f\n",eq,yy);
y[eq] += slowc_l*yy;
#ifdef PRINT_OUT_y1
mexPrintf("=%f (%f)\n",double(yy),double(y[eq]));
#endif
}
}
return res1;
}
double
SparseMatrix::simple_bksub(int it_, int Size, double slowc_l)
{
int i,k;
double yy;
NonZeroElem *first;
res1 = res2 = max_res = 0;
//mexPrintf("simple_bksub\n");
for (i=0;i<y_size;i++)
y[i+it_*y_size]=ya[i+it_*y_size];
//mexPrintf("setp 1\n");
//mexPrintf("Size=%d\n",Size);
for(i=Size-1;i>=0;i--)
{
//mexPrintf("i=%d\n",i);
int pos=pivot[i];
//mexPrintf("pos=%d\n",pos);
int nb_var=At_Row(pos,&first);
//mexPrintf("nb_var=%d\n",nb_var);
first=first->NZE_R_N;
nb_var--;
//mexPrintf("i=%d\n",i);
int eq=index_vara[i];
//mexPrintf("eq=%d\n",eq);
yy = 0;
for (k=0;k<nb_var;k++)
{
//mexPrintf("y[index_vara[%d]=%d]=%f\n",first->c_index, index_vara[first->c_index], y[index_vara[first->c_index]+it_*y_size]);
//mexPrintf("u[first->u_index=%d]=%f\n",first->u_index, u[first->u_index]);
yy+=y[index_vara[first->c_index]+it_*y_size]*u[first->u_index];
first=first->NZE_R_N;
}
yy=-(yy+y[eq+it_*y_size]+u[b[pos]]);
//mexPrintf("yy=%f\n",yy);
direction[eq+it_*y_size]=yy;
y[eq+it_*y_size] += slowc_l*yy;
}
return res1;
}
int
SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter)
{
int i, j, k;
int pivj=0, pivk=0;
double piv_abs, first_elem;
NonZeroElem *first, *firsta, *first_sub, *first_piv, *first_suba;
#ifdef MARKOVITZ
double piv_v[Size];
int pivj_v[Size], pivk_v[Size], NR[Size], l, N_max;
bool one;
#endif
/*mexPrintf("In simulate_NG it_=%d\n",it_);
mexEvalString("drawnow;");*/
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
if (isnan(res1) || isinf(res1))
{
if (slowc_save<1e-8)
{
mexPrintf("Dynare cannot improve the simulation\n");
mexEvalString("drawnow;");
mexEvalString("st=fclose('all');clear all;");
filename+=" stopped";
mexErrMsgTxt(filename.c_str());
}
slowc_save/=2;
mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n",slowc_save);
for (i=0;i<y_size;i++)
y[i+it_*y_size]=ya[i+it_*y_size]+slowc_save*direction[i+it_*y_size];
iter--;
return(0);
}
//mexPrintf("before the main loop y_size=%d\n",y_size);
for (i=0;i<Size;i++)
{
/*finding the max-pivot*/
double piv=piv_abs=0;
//mexPrintf("i=%d\n",i);
int nb_eq=At_Col(i, 0, &first);
//mexPrintf("nb_eq=%d\n",nb_eq);
#ifdef MARKOVITZ
l=0; N_max=0;
one=false;
piv_abs=0;
#endif
//mexPrintf("Size=%d\n",Size);
//#pragma omp parallel for
for (j=0;j<Size;j++)
{
//mexPrintf("j=%d ",j);
//mexPrintf("first->r_index=%d\n", first->r_index);
if (!line_done[first->r_index])
{
k=first->u_index;
//mexPrintf("k=%d\n",k);
int jj=first->r_index;
//mexPrintf("jj=%d\n",jj);
int NRow_jj=NRow(jj);
//mexPrintf("NRow_jj=%d\n",NRow_jj);
#ifdef PROFILER
nbpivot++;
nbpivot_it++;
#endif
#ifdef MARKOVITZ
//mexPrintf("l=%d\n",l);
piv_v[l]=u[k];
double piv_fabs=fabs(u[k]);
pivj_v[l]=jj;
pivk_v[l]=k;
NR[l]=NRow_jj;
//mexPrintf("one=%d\n",one);
if (NRow_jj==1 && !one)
{
one=true;
piv_abs=piv_fabs;
N_max=NRow_jj;
}
if (!one)
{
if (piv_fabs>piv_abs)
piv_abs=piv_fabs;
if (NRow_jj>N_max)
N_max=NRow_jj;
}
else
{
if (NRow_jj==1)
{
if (piv_fabs>piv_abs)
piv_abs=piv_fabs;
if (NRow_jj>N_max)
N_max=NRow_jj;
}
}
l++;
#else
if (piv_abs<fabs(u[k])||NRow_jj==1)
{
piv=u[k];
piv_abs=fabs(piv);
pivj=jj; //Line number
pivk=k; //position in u
if (NRow_jj==1)
break;
}
#endif
}
//mexPrintf("j=%d first->NZE_C_N=%x\n",j,first->NZE_C_N);
first=first->NZE_C_N;
}
//mexPrintf("Before Markovitz\n");
#ifdef MARKOVITZ
double markovitz=0, markovitz_max=-9e70;
//mexPrintf("l=%d one=%d\n",l, one);
if (!one)
{
for (j=0;j<l;j++)
{
/*mexPrintf("j=%d\n",j);
mexPrintf("piv_abs=%f\n",piv_abs);
mexPrintf("piv_v[%d]=%f\n",j,piv_v[j]);
mexPrintf("N_max=%d\n",N_max);
mexPrintf("NR[%d]=%d\n",j,NR[j]);*/
markovitz=exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double(NR[j])/double(N_max)));
//mexPrintf("markovitz=%f\n",markovitz);
if (markovitz>markovitz_max)
{
piv=piv_v[j];
pivj=pivj_v[j]; //Line number
pivk=pivk_v[j]; //positi
markovitz_max=markovitz;
}
}
}
else
{
for (j=0;j<l;j++)
{
markovitz=exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(NR[j]/N_max));
if (markovitz>markovitz_max && NR[j]==1)
{
piv=piv_v[j];
pivj=pivj_v[j]; //Line number
pivk=pivk_v[j]; //positi
markovitz_max=markovitz;
}
}
}
#endif
//mexPrintf("OK i=%d\n", i);
pivot[i]=pivj;
//mexPrintf("pivot[%d]=%d\n",i,pivot[i]);
pivotk[i]=pivk;
pivotv[i]=piv;
line_done[pivj]=true;
if (piv_abs<eps)
{
mexPrintf("Error: singular system in Simulate_NG\n");
mexEvalString("drawnow;");
mexEvalString("st=fclose('all');clear all;");
filename+=" stopped";
mexErrMsgTxt(filename.c_str());
}
/*divide all the non zeros elements of the line pivj by the max_pivot*/
int nb_var=At_Row(pivj,&first);
for (j=0;j<nb_var;j++)
{
u[first->u_index]/=piv;
first=first->NZE_R_N;
}
u[b[pivj]]/=piv;
/*substract the elements on the non treated lines*/
nb_eq=At_Col(i,&first);
NonZeroElem *first_piva;
int nb_var_piva=At_Row(pivj,&first_piva);
for (j=0;j<Size;j++)
{
int row=first->r_index;
if (!line_done[row])
{
first_elem=u[first->u_index];
int nb_var_piv=nb_var_piva;
first_piv=first_piva;
int nb_var_sub=At_Row(row,&first_sub);
int l_sub=0, l_piv=0;
int sub_c_index=first_sub->c_index, piv_c_index=first_piv->c_index;
//int tmp_lag=first_sub->lag_index;
while (l_sub<nb_var_sub || l_piv<nb_var_piv)
{
if (l_sub<nb_var_sub && (sub_c_index<piv_c_index || l_piv>=nb_var_piv))
{
first_sub=first_sub->NZE_R_N;
if (first_sub)
sub_c_index=first_sub->c_index;
else
sub_c_index=Size/* *periods*/;
l_sub++;
}
else if (sub_c_index>piv_c_index || l_sub>=nb_var_sub)
{
int tmp_u_count=Get_u();
//int lag=first_piv->c_index/Size-row/Size;
Insert(row,first_piv->c_index,tmp_u_count,/*lag*/0);
u[tmp_u_count]=-u[first_piv->u_index]*first_elem;
first_piv=first_piv->NZE_R_N;
if (first_piv)
piv_c_index=first_piv->c_index;
else
piv_c_index=Size/* *periods*/;
l_piv++;
}
else /*first_sub->c_index==first_piv->c_index*/
{
if (i==sub_c_index)
{
firsta=first;
first_suba=first_sub->NZE_R_N;
Delete(first_sub->r_index,first_sub->c_index, Size, b);
first=firsta->NZE_C_N;
first_sub=first_suba;
if (first_sub)
sub_c_index=first_sub->c_index;
else
sub_c_index=Size/* *periods*/;
l_sub++;
first_piv=first_piv->NZE_R_N;
if (first_piv)
piv_c_index=first_piv->c_index;
else
piv_c_index=Size;
l_piv++;
}
else
{
u[first_sub->u_index]-=u[first_piv->u_index]*first_elem;
first_sub=first_sub->NZE_R_N;
if (first_sub)
sub_c_index=first_sub->c_index;
else
sub_c_index=Size /* *periods*/;
l_sub++;
first_piv=first_piv->NZE_R_N;
if (first_piv)
piv_c_index=first_piv->c_index;
else
piv_c_index=Size/* *periods*/;
l_piv++;
}
}
}
u[b[row]]-=u[b[pivj]]*first_elem;
}
else
first=first->NZE_C_N;
}
}
double slowc_lbx=slowc, res1bx;
//mexPrintf("before bksub it_=%d\n",it_);
for (i=0;i<y_size;i++)
ya[i+it_*y_size]=y[i+it_*y_size];
slowc_save=slowc;
//res1bx=bksub( NULL, 1, y_size, slowc_lbx);
res1bx=simple_bksub(it_,Size,slowc_lbx);
//mexPrintf("End of simulate_NG\n");
return(0);
}
@ -1976,6 +2384,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
double piv_abs, first_elem;
//SparseMatrix sparse_matrix;
//mexPrintf("->u_count=%d &u_count=%x\n",u_count,&u_count);
mexPrintf("GNU version=%d\n",GNUVER);
#ifdef RECORD_ALL
int save_u_count=u_count;
#endif
@ -2208,8 +2617,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
bool one=false;
piv_abs=0;
#endif
//#pragma omp parallel for
for (j=0;j<nb_eq;j++)
{
//mexPrintf("j=%d\n",j);
#ifdef PRINT_OUT
mexPrintf("first=%x\n",first);
mexPrintf("examine col %d row %d with lag 0 line_done=%d \n",i,first->r_index,line_done[first->r_index]);
@ -2275,6 +2686,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
double markovitz=0, markovitz_max=-9e70;
if (!one)
{
/*#pragma omp parallel for private(j) shared(markovitz_max, piv, pivj, pivk)
for (j=0;j<l;j++)
{
markovitz=exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double(NR[j])/double(N_max)));
@ -2285,19 +2697,45 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
pivk=pivk_v[j]; //positi
markovitz_max=markovitz;
}
}*/
//#pragma omp parallel private(j, markovitz)
{
//#pragma omp for schedule(dynamic, 10) nowait
//#pragma omp parallel for private(j, markovitz)
#pragma omp parallel for private(j,markovitz)
for (j=0;j<l;j++)
{
markovitz=exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double(NR[j])/double(N_max)));
//printf("j=%d, markovitz=%f, markovitz_max=%f\n",j,markovitz, markovitz_max);
if (markovitz>markovitz_max)
{
piv=piv_v[j];
pivj=pivj_v[j]; //Line number
pivk=pivk_v[j]; //positi
markovitz_max=markovitz;
//printf("j=%d, markovitz_max=%f\n",j, markovitz);
}
}
}
}
else
{
for (j=0;j<l;j++)
//#pragma omp parallel for private(j) shared(markovitz_max, piv, pivj, pivk)
//#pragma omp parallel private(j, markovitz)
{
markovitz=exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(NR[j]/N_max));
if (markovitz>markovitz_max && NR[j]==1)
//#pragma omp for schedule(dynamic, 10) nowait
//#pragma omp parallel for private(j, markovitz)
#pragma omp parallel for private(j,markovitz)
for (j=0;j<l;j++)
{
piv=piv_v[j];
pivj=pivj_v[j]; //Line number
pivk=pivk_v[j]; //positi
markovitz_max=markovitz;
markovitz=exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(NR[j]/N_max));
if (markovitz>markovitz_max && NR[j]==1)
{
piv=piv_v[j];
pivj=pivj_v[j]; //Line number
pivk=pivk_v[j]; //positi
markovitz_max=markovitz;
}
}
}
}
@ -2373,7 +2811,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
#endif
if (piv_abs<eps)
{
mexPrintf("Error: singular system\n");
mexPrintf("Error: singular system in Simulate_NG1\n");
mexEvalString("drawnow;");
mexEvalString("st=fclose('all');clear all;");
filename+=" stopped";
@ -2388,6 +2826,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
#ifdef PRINT_OUT
mexPrintf("nb_var=%d\n",nb_var);
#endif
//#pragma omp parallel for
for (j=0;j<nb_var;j++)
{
#ifdef PRINT_OUT
@ -2520,6 +2959,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
mexEvalString("drawnow;");
}
#endif
//#pragma omp parallel for
for (j=0;j<nb_eq;j++)
{
int row=first->r_index;
@ -3099,8 +3539,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
void
SparseMatrix::fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1)
{
//mexPrintf("u_count_int=%d\n",u_count_int);
u_count=u_count_int * periods;
u_count_alloc = 2*u_count;
//mexPrintf("mxMalloc(%d*sizeof(double)=%d)\n",u_count_alloc,u_count_alloc*sizeof(double));
(*u)=(double*)mxMalloc(u_count_alloc*sizeof(double));
memset((*u), 0, u_count_alloc*sizeof(double));
u_count_init=max_lag_plus_max_lead_plus_1;

10
mex/sources/simulate/SparseMatrix.hh
View File

@ -26,7 +26,12 @@
#include <map>
#include <ctime>
#include "Mem_Mngr.hh"
//! Openmp is available in GCC since version 4.3.2
//! Test if GCC version is greater then 4.3.2 order to avoid a copilation error
#define GNUVER 100*__GNUC__+10*__GNUC_MINOR__+__GNUC_PATCHLEVEL__
#if GNUVER >= 432
#include <omp.h>
#endif
#define NEW_ALLOC
#define MARKOVITZ
//#define MEMORY_LEAKS
@ -71,6 +76,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 simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter);
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);
//void initialize(int periods_arg, int nb_endo_arg, int y_kmin_arg, int y_kmax_arg, int y_size_arg, int u_count_arg, int u_count_init_arg, double *u_arg, double *y_arg, double *ya_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg, double res1_arg, double res2_arg, double max_res_arg);
@ -80,6 +86,7 @@ class SparseMatrix
private:
void Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int> ,int>, int> IM);
void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int> ,int>, int> IM);
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(int Size);
bool compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
#ifdef PROFILER
@ -106,6 +113,7 @@ class SparseMatrix
, long int *nmul
#endif
);
double simple_bksub(int it_, int Size, double slowc_l);
void run_triangular(int nop_all,int *op_all);
void run_it(int nop_all,int *op_all);
void run_u_period1(int periods);

17
preprocessor/ExprNode.cc
View File

@ -147,10 +147,11 @@ NumConstNode::eval(const eval_context_type &eval_context) const throw (EvalExcep
}
void
NumConstNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const
NumConstNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
{
CompileCode.write(&FLDC, sizeof(FLDC));
double vard=atof(datatree.num_constants.get(id).c_str());
//double vard=id;
#ifdef DEBUGC
cout << "FLDC " << vard << "\n";
#endif
@ -427,7 +428,7 @@ VariableNode::eval(const eval_context_type &eval_context) const throw (EvalExcep
}
void
VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const
VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
{
int i, lagl;
#ifdef DEBUGC
@ -892,7 +893,7 @@ UnaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExcept
}
void
UnaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const
UnaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
{
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode *>(this));
if (it != temporary_terms.end())
@ -1164,8 +1165,6 @@ BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
reference_count[this2]++;
if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C)
{
if(this2->idx==2280)
cout << "==>Curr_block= " << Curr_block << " equation= " << equation << " first_occurence[this2].first=" << first_occurence[this2].first << " first_occurence[this2].second=" << first_occurence[this2].second << "\n";
temporary_terms.insert(this2);
ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2);
}
@ -1226,7 +1225,7 @@ BinaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExcep
}
void
BinaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const
BinaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
{
// If current node is a temporary term
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
@ -1580,8 +1579,6 @@ TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
reference_count[this2]++;
if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C)
{
if(this2->idx==2280)
cout << "==>Curr_block= " << Curr_block << " equation= " << equation << " first_occurence[this2].first=" << first_occurence[this2].first << " first_occurence[this2].second=" << first_occurence[this2].second << "\n";
temporary_terms.insert(this2);
ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2);
}
@ -1613,7 +1610,7 @@ TrinaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExce
void
TrinaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms, map_idx_type map_idx) const
const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
{
// If current node is a temporary term
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<TrinaryOpNode *>(this));
@ -1788,7 +1785,7 @@ UnknownFunctionNode::eval(const eval_context_type &eval_context) const throw (Ev
}
void
UnknownFunctionNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const
UnknownFunctionNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
{
cerr << "UnknownFunctionNode::compile: operation impossible!" << endl;
exit(EXIT_FAILURE);

1
preprocessor/Makefile
View File

@ -33,7 +33,6 @@ OBJS = \
IncidenceMatrix.o \
BlockTriangular.o \
Model_Graph.o \
SymbolGaussElim.o \
DynareMain.o \
DynareMain2.o

14
preprocessor/ModFile.cc
View File

@ -125,7 +125,7 @@ ModFile::evalAllExpressions()
cout << "error in evaluation of variable\n";
}
}
if(mod_file_struct.load_params_and_steady_state_present && init_values.size()<symbol_table.endo_nbr+symbol_table.exo_nbr+symbol_table.exo_det_nbr)
if(mod_file_struct.load_params_and_steady_state_present && int(init_values.size())<symbol_table.endo_nbr+symbol_table.exo_nbr+symbol_table.exo_det_nbr)
{
for(j=0;j <symbol_table.endo_nbr; j++)
{
@ -200,7 +200,7 @@ ModFile::evalAllExpressions()
}
}
}
if(init_values.size()<symbol_table.endo_nbr+symbol_table.exo_nbr+symbol_table.exo_det_nbr)
if(int(init_values.size())<symbol_table.endo_nbr+symbol_table.exo_nbr+symbol_table.exo_det_nbr)
{
cout << "\nWarning: Uninitialized variable: \n";
cout << "Endogenous\n";
@ -246,7 +246,7 @@ ModFile::evalAllExpressions()
cout << "error in evaluation of pound\n";
}
}
if(model_tree.local_variables_table.size()!=symbol_table.model_local_variable_nbr+symbol_table.modfile_local_variable_nbr)
if(int(model_tree.local_variables_table.size())!=symbol_table.model_local_variable_nbr+symbol_table.modfile_local_variable_nbr)
{
cout << "Warning: Unitilialized pound: \n";
cout << "Local variable in a model\n";
@ -398,7 +398,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
mOutputFile << "logname_ = '" << basename << ".log';" << endl;
mOutputFile << "diary " << basename << ".log" << endl;
mOutputFile << "options_.model_mode = " << model_tree.mode << ";\n";
if (model_tree.mode == eSparseMode)
if (model_tree.mode == eSparseMode || model_tree.mode == eSparseDLLMode)
{
mOutputFile << "addpath " << basename << ";\n";
mOutputFile << "delete('" << basename << "_static.m');\n";
@ -408,7 +408,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
if (model_tree.equation_number() > 0)
{
if (model_tree.mode == eDLLMode)
if (model_tree.mode == eDLLMode || model_tree.mode == eSparseDLLMode)
{
mOutputFile << "if exist('" << basename << "_static.c')" << endl;
mOutputFile << " clear " << basename << "_static" << endl;
@ -425,7 +425,6 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
mOutputFile << "erase_compiled_function('" + basename +"_dynamic');" << endl;
}
}
cout << "Processing outputs ...";
symbol_table.writeOutput(mOutputFile);
@ -448,9 +447,10 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
mOutputFile << "save('" << basename << "_results.mat', 'oo_', 'M_', 'options_');" << endl;
mOutputFile << "diary off" << endl;
if (model_tree.mode == eSparseMode)
if (model_tree.mode == eSparseMode || model_tree.mode == eSparseDLLMode)
mOutputFile << "rmpath " << basename << ";\n";
mOutputFile << endl << "disp(['Total computing time : ' dynsec2hms(toc) ]);" << endl;
mOutputFile.close();
cout << "done\n";
}

60
preprocessor/ModelTree.cc
View File

@ -73,7 +73,7 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
}
void
ModelTree::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type map_idx) const
ModelTree::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const
{
first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(eEndogenous, symb_id, lag)));
if (it != first_derivatives.end())
@ -1082,7 +1082,7 @@ end:
void
ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type) const
ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const
{
struct Uff_l
{
@ -1107,7 +1107,8 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
map<NodeID, int> reference_count;
map<int,int> ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number;
int prev_Simulation_Type=-1;
SymbolicGaussElimination SGE;
//SymbolicGaussElimination SGE;
bool file_open=false;
temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
//----------------------------------------------------------------------
string main_name=file_name;
@ -1180,14 +1181,15 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
v=block_triangular.ModelBlock->Block_List[j].Max_Lead;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
{
//if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
//{
int u_count_int=0;
Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,SGE.file_open);
Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open,
ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE);
v=u_count_int;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
SGE.file_is_open();
}
file_open=true;
//}
}
for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++)
{
@ -1821,7 +1823,7 @@ ModelTree::reform(const string name1) const
void
ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, const int &num,
int &u_count_int, bool &file_open) const
int &u_count_int, bool &file_open, bool is_two_boundaries) const
{
int j;
std::ofstream SaveCode;
@ -1850,19 +1852,22 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b
u_count_int++;
}
}
for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
if(is_two_boundaries)
{
int eqr1=j;
int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods
+block_triangular.incidencematrix.Model_Max_Lead_Endo);
int k1=0;
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
u_count_int++;
for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
{
int eqr1=j;
int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods
+block_triangular.incidencematrix.Model_Max_Lead_Endo);
int k1=0;
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
u_count_int++;
}
}
//cout << "u_count_int=" << u_count_int << "\n";
for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
{
int varr=block_triangular.ModelBlock->Block_List[num].Variable[j];
@ -2103,7 +2108,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
ofstream mDynamicModelFile;
ostringstream tmp, tmp1, tmp_eq;
int prev_Simulation_Type, tmp_i;
SymbolicGaussElimination SGE;
//SymbolicGaussElimination SGE;
bool OK;
chdir(basename.c_str());
string filename = dynamic_basename + ".m";
@ -2121,7 +2126,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
mDynamicModelFile << "%/\n";
int i, k, Nb_SGE=0;
bool printed = false, skip_head, open_par=false;
bool skip_head, open_par=false;
if (computeJacobian || computeJacobianExo || computeHessian)
{
mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n";
@ -2413,10 +2418,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
if (open_par)
mDynamicModelFile << " end\n";
open_par=false;
if (!printed)
{
printed = true;
}
Nb_SGE++;
int nze, m;
for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++)
@ -2453,8 +2454,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
writeModelEquationsOrdered_M( block_triangular.ModelBlock, dynamic_basename);
chdir("..");
if (printed)
cout << "done\n";
}
void
@ -3074,9 +3073,10 @@ ModelTree::writeStaticFile(const string &basename) const
switch (mode)
{
case eStandardMode:
case eSparseDLLMode:
/*case eSparseDLLMode:*/
writeStaticMFile(basename + "_static");
break;
case eSparseDLLMode:
case eSparseMode:
// create a directory to store all files
#ifdef _WIN32
@ -3117,7 +3117,7 @@ ModelTree::writeDynamicFile(const string &basename) const
#else
mkdir(basename.c_str(), 0777);
#endif
writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL);
writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL, map_idx);
block_triangular.Free_Block(block_triangular.ModelBlock);
block_triangular.incidencematrix.Free_IM();
//block_triangular.Free_IM_X(block_triangular.First_IM_X);

2
preprocessor/ParsingDriver.cc
View File

@ -596,7 +596,7 @@ ParsingDriver::add_to_row(NodeID v)
void
ParsingDriver::steady()
{
if (mod_file->model_tree.mode == eSparseMode)
if (mod_file->model_tree.mode == eSparseMode || mod_file->model_tree.mode == eSparseDLLMode)
mod_file->addStatement(new SteadySparseStatement(options_list));
else
mod_file->addStatement(new SteadyStatement(options_list));

1827
preprocessor/SymbolGaussElim.cc
View File

File diff suppressed because it is too large Load Diff

14
preprocessor/include/ExprNode.hh
View File

@ -159,7 +159,7 @@ public:
};
virtual double eval(const eval_context_type &eval_context) const throw (EvalException) = 0;
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const = 0;
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const = 0;
};
//! Object used to compare two nodes (using their indexes)
@ -186,7 +186,7 @@ public:
virtual void collectExogenous(set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
};
//! Symbol or variable node
@ -214,7 +214,7 @@ public:
map_idx_type &map_idx) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
};
//! Unary operator node
@ -243,7 +243,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
};
//! Binary operator node
@ -273,7 +273,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
virtual NodeID get_arg1() { return(arg1);};
virtual NodeID get_arg2() { return(arg2);};
};
@ -310,7 +310,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
};
//! Unknown function node
@ -337,7 +337,7 @@ public:
virtual void collectExogenous(set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
};
//! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model

2
preprocessor/include/ModFile.hh
View File

@ -23,6 +23,8 @@
using namespace std;
#include <ostream>
#include <ctime>
#include "SymbolTable.hh"
#include "NumericalConstants.hh"

7
preprocessor/include/ModelTree.hh
View File

@ -32,7 +32,6 @@ using namespace std;
#include "NumericalConstants.hh"
#include "DataTree.hh"
#include "BlockTriangular.hh"
#include "SymbolGaussElim.hh"
//! The three in which ModelTree can work
enum ModelTreeMode
@ -85,7 +84,7 @@ private:
//! Write derivative of an equation w.r. to a variable
void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, SymbolType type) const;
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type map_idx) const;
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const;
//! Computes temporary terms
void computeTemporaryTerms(int order);
void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock);
@ -109,7 +108,7 @@ private:
//! Writes the Block reordred structure of the static model in M output
void writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const string &static_basename) const;
//! Writes the code of the Block reordred structure of the model in virtual machine bytecode
void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type) const;
void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const;
//! Writes static model file (Matlab version)
void writeStaticMFile(const string &static_basename) const;
//! Writes static model file (C version)
@ -173,7 +172,7 @@ public:
BlockTriangular block_triangular;
//! Adds informations for simulation in a binary file
void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename,
const int &num, int &u_count_int, bool &file_open) const;
const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
//! Returns the number of equations in the model
int equation_number() const;
};

8
tests/ferhat/ls2003.mod
View File

@ -17,8 +17,8 @@ rho_ys = 0.9;
rho_pies = 0.7;
//model(sparse_dll,cutoff=1e-17);
model(sparse);
model(sparse_dll);
//model(sparse);
//model;
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
@ -84,7 +84,7 @@ options_.maxit_=100;
steady;
model_info;
check;
//check;
shocks;
var e_q;
@ -99,5 +99,5 @@ rplot ww;
rplot A;
rplot pie;
stoch_simul(periods=200,order=1);
//stoch_simul(periods=200,order=1);

4
tests/ferhat/ramst.mod
View File

@ -10,8 +10,8 @@ aa=0.5;
//model(sparse);
//model(sparse_dll);
model;
model(sparse_dll);
//model;
c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
end;