diff --git a/preprocessor/BlockTriangular.cc b/preprocessor/BlockTriangular.cc index a87ce54bf..336334625 100644 --- a/preprocessor/BlockTriangular.cc +++ b/preprocessor/BlockTriangular.cc @@ -147,59 +147,23 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc } else { - k = -k; if(Cur_IM[i_1 + Index_Var_IM[*count_Equ].index]) { - tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; nb_lead_lag_endo++; - if(k > Lag) - Lag = k; + if(-k > Lag) + Lag = -k; } } } } - - /*Cur_IM = incidencematrix.Get_First(eEndogenous); - while(Cur_IM) - { - k = Cur_IM->lead_lag; - i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr; - if(k > 0) - { - if(Cur_IM->IM[i_1 + Index_Var_IM[*count_Equ].index]) - { - nb_lead_lag_endo++; - tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; - if(k > Lead) - Lead = k; - } - } - else - { - k = -k; - if(Cur_IM->IM[i_1 + Index_Var_IM[*count_Equ].index]) - { - tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++; - nb_lead_lag_endo++; - if(k > Lag) - Lag = k; - } - } - Cur_IM = Cur_IM->pNext; - }*/ - Lag_Endo = Lag; Lead_Endo = Lead; tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); tmp_nb_exo = 0; - - /*Cur_IM = incidencematrix.Get_First(eExogenous); - k = Cur_IM->lead_lag; - while(Cur_IM) - {*/ for(k = -incidencematrix.Model_Max_Lag_Exo; k<=incidencematrix.Model_Max_Lead_Exo; k++) { Cur_IM = incidencematrix.Get_IM(k, eExogenous); @@ -257,23 +221,22 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc else ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE; - //Cur_IM = incidencematrix.Get_First(eExogenous); tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); tmp_nb_exo = 0; - /*while(Cur_IM) - {*/ for(k = -incidencematrix.Model_Max_Lag_Exo; k <=incidencematrix.Model_Max_Lead_Exo; k++) { Cur_IM = incidencematrix.Get_IM(k, eExogenous); - i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr; - for(j=0;jpNext; + if(Cur_IM) + { + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr; + for(j=0;jBlock_List[*count_Block].nb_exo = tmp_nb_exo; ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int)); @@ -295,7 +258,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc { if(incidencematrix.Model_Max_Lag_Endo - Lag + li>=0) { - //cout << "ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag[" << li << "].size=" << tmp_size[Model_Max_Lag_Endo - Lag + li] << "\n"; ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li]; ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_endo = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li]; ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); @@ -339,7 +301,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc if(incidencematrix.Model_Max_Lag_Exo - Lag + li>=0) { ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_exo = tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li]; - //cout << "ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag[" << li << "].Exogenous= malloc(" << tmp_size_exo[Model_Max_Lag_Exo - Lag + li] << ")\n"; ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); @@ -384,7 +345,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].Simulation_Type = UNKNOWN; ModelBlock->Block_List[*count_Block].Equation = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); ModelBlock->Block_List[*count_Block].Variable = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); - //ModelBlock->Block_List[*count_Block].Variable_Sorted = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); ModelBlock->Block_List[*count_Block].Own_Derivative = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); Lead = Lag = 0; first_count_equ = *count_Equ; @@ -402,12 +362,9 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc { ModelBlock->Block_List[*count_Block].Equation[i] = Index_Equ_IM[*count_Equ].index; ModelBlock->Block_List[*count_Block].Variable[i] = Index_Var_IM[*count_Equ].index; - //Cur_IM = incidencematrix.Get_First(eEndogenous); i_1 = Index_Var_IM[*count_Equ].index; - //while(Cur_IM) for(k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++) { - //k = Cur_IM->lead_lag; Cur_IM = incidencematrix.Get_IM(k, eEndogenous); if(Cur_IM) { @@ -432,20 +389,19 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc } else { - k = -k; for(j = 0;j < size;j++) { if(Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr]) { - tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; if (!OK) { - tmp_endo[incidencematrix.Model_Max_Lag - k]++; + tmp_endo[incidencematrix.Model_Max_Lag + k]++; nb_lead_lag_endo++; OK = true; } - if(k > Lag) - Lag = k; + if(-k > Lag) + Lag = -k; } } } @@ -453,6 +409,8 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc } (*count_Equ)++; } + + if ((Lag > 0) && (Lead > 0)) ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE; else if(size > 1) @@ -478,9 +436,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc tmp_nb_exo = 0; for(i = 0;i < size;i++) { - //Cur_IM = incidencematrix.Get_First(eExogenous); - //k = Cur_IM->lead_lag; - //while(Cur_IM) for(k = -incidencematrix.Model_Max_Lag_Exo; k<=incidencematrix.Model_Max_Lead_Exo; k++) { Cur_IM = incidencematrix.Get_IM(k, eExogenous); @@ -860,13 +815,16 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_ //First create a static model incidence matrix size = symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(*SIM); SIM = (bool*)malloc(size); - memset(SIM, size, 0); + for(i = 0; i< symbol_table.endo_nbr * symbol_table.endo_nbr; i++) SIM[i] = 0; for(k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++) { Cur_IM = incidencematrix.Get_IM(k, eEndogenous); - for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) + if(Cur_IM) { - SIM[i] = (SIM[i]) || (Cur_IM[i]); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) + { + SIM[i] = (SIM[i]) || (Cur_IM[i]); + } } } if(bt_verbose) diff --git a/preprocessor/InidenceMatrix.cc b/preprocessor/IncidenceMatrix.cc similarity index 97% rename from preprocessor/InidenceMatrix.cc rename to preprocessor/IncidenceMatrix.cc index 65ac87255..286eb1a26 100644 --- a/preprocessor/InidenceMatrix.cc +++ b/preprocessor/IncidenceMatrix.cc @@ -38,7 +38,7 @@ IncidenceMatrix::Build_IM(int lead_lag, SymbolType type) { size = symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(IM[0]); List_IM[lead_lag] = IM = (bool*)malloc(size); - memset(IM, size, NULL); + for(int i = 0; i< symbol_table.endo_nbr * symbol_table.endo_nbr; i++) IM[i] = 0; if(lead_lag > 0) { if(lead_lag > Model_Max_Lead_Endo) @@ -62,7 +62,7 @@ IncidenceMatrix::Build_IM(int lead_lag, SymbolType type) { //eExogenous size = symbol_table.endo_nbr * symbol_table.exo_nbr * sizeof(IM[0]); List_IM_X[lead_lag] = IM = (bool*)malloc(size); - memset(IM, size, NULL); + for(int i = 0; i< symbol_table.endo_nbr * symbol_table.exo_nbr; i++) IM[i] = 0; if(lead_lag > 0) { if(lead_lag > Model_Max_Lead_Exo) diff --git a/preprocessor/Makefile b/preprocessor/Makefile index d1a85afe8..1712c2ac9 100644 --- a/preprocessor/Makefile +++ b/preprocessor/Makefile @@ -29,6 +29,7 @@ OBJS = \ ExprNode.o \ ModelNormalization.o \ ModelBlocks.o \ + IncidenceMatrix.o \ BlockTriangular.o \ Model_Graph.o \ SymbolGaussElim.o \ diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index fdc79d687..ee7fd9db0 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -396,7 +396,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms); - /*tmp_output << "[" << block_triangular.periods + variable_table.max_lag+variable_table.max_lead << "]";*/ } global_output << " global " << tmp_output.str() << " M_ ;\n"; //For each block @@ -588,8 +587,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock for(i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++) { int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; - //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; - //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; output << " g1(" << eq+1 << ", " << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr/*ModelBlock->Block_List[j].nb_exo*/ << ") = "; @@ -629,12 +626,10 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock { int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; - //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; output << " g1(" << eq+1 << ", " << var+1+(m+variable_table.max_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ") = "; writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << /*variable_table.getLag(variable_table.getSymbolID(var))*/k + << "(" << k << ") " << var+1 << ", equation=" << eq+1 << endl; } @@ -646,8 +641,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock for(i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++) { int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; - //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; - //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; output << " g1(" << eq+1 << ", " << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = "; writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous); @@ -1883,7 +1876,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b string filename; ofstream mStaticModelFile; int i, k, prev_Simulation_Type, ga_index = 1; - bool /*printed = false,*/ skip_head, open_par=false; + bool skip_head, open_par=false; filename = static_basename + ".m"; mStaticModelFile.open(filename.c_str(), ios::out | ios::binary); if (!mStaticModelFile.is_open()) @@ -2132,10 +2125,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " y_size=M_.endo_nbr;\n"; mDynamicModelFile << " if(length(varargin)>0)\n"; mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n"; - //mDynamicModelFile << " global it_;\n"; mDynamicModelFile << " it_=varargin{3};\n"; - //mDynamicModelFile << " g=zeros(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n"; - //mDynamicModelFile << " g1=spalloc(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1),y_size*y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n"; i = symbol_table.endo_nbr*(variable_table.max_endo_lag+variable_table.max_endo_lead+1)+ symbol_table.exo_nbr*(variable_table.max_exo_lag+variable_table.max_exo_lead+1); mDynamicModelFile << " g1=spalloc(" << symbol_table.endo_nbr << ", " << i << ", " << i*symbol_table.endo_nbr << ");\n"; @@ -2357,26 +2347,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " end\n"; mDynamicModelFile << " end\n"; } - /*else if ((k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (open_par) - mDynamicModelFile << " end\n"; - open_par=false; - if (!printed) - { - printed = true; - } - SGE.SGE_compute(block_triangular.ModelBlock, i, true, basename, symbol_table.endo_nbr); - Nb_SGE++; -#ifdef PRINT_OUT - cout << "end of Gaussian elimination\n"; -#endif - mDynamicModelFile << " Read_file(\"" << reform(basename) << "\",periods," << - block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ", " << symbol_table.endo_nbr << - ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << ");\n"; - cerr << "Not implemented block SOLVE_TWO_BOUNDARIES_COMPLETE" << endl; - exit(EXIT_FAILURE); - }*/ else if ((k == SOLVE_FORWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size)) { if (open_par) @@ -2406,8 +2376,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " end\n"; mDynamicModelFile << " end\n"; - /*cerr << "Not implemented block SOLVE_FORWARD_COMPLETE" << endl; - exit(EXIT_FAILURE);*/ } else if ((k == SOLVE_BACKWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size)) { @@ -2437,8 +2405,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " return;\n"; mDynamicModelFile << " end\n"; mDynamicModelFile << " end\n"; - /*cerr << "Not implemented block SOLVE_BACKWARD_COMPLETE" << endl; - exit(EXIT_FAILURE);*/ } else if ((k == SOLVE_TWO_BOUNDARIES_COMPLETE || k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) { @@ -2875,7 +2841,6 @@ ModelTree::writeOutput(ostream &output) const k++; count_lead_lag_incidence = 0; int Block_size=block_triangular.ModelBlock->Block_List[j].Size; - //int Block_exo_nbr=block_triangular.ModelBlock->Block_List[j].exo_nbr; max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ; max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead; max_lag_endo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Endo ; @@ -2883,7 +2848,7 @@ ModelTree::writeOutput(ostream &output) const max_lag_exo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Exo ; max_lead_exo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Exo; bool evaluate=false; - vector exogenous;//, tmp_exogenous(block_triangular.exo_nbr,-1); + vector exogenous; vector::iterator it_exogenous; exogenous.clear(); ostringstream tmp_s, tmp_s_eq; @@ -2935,7 +2900,6 @@ ModelTree::writeOutput(ostream &output) const for(int i=0;iBlock_List[j+Block_size].nb_exo;i++) { int ii=block_triangular.ModelBlock->Block_List[j+Block_size].Exogenous[i]; - //for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) cout << "*it_exogenous=" << *it_exogenous << "\n"; if(it_exogenous==exogenous.end()) exogenous.push_back(ii); } @@ -3038,8 +3002,7 @@ ModelTree::writeOutput(ostream &output) const } ii++; } - //if(done_some_where) - output << tmp_s.str() << "\n"; + output << tmp_s.str() << "\n"; tmp_s.str(""); } } @@ -3055,6 +3018,7 @@ ModelTree::writeOutput(ostream &output) const if(IM) { bool new_entry=true; + output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").lead_lag = " << j << ";\n"; output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").sparse_IM = ["; for(int i=0;i