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

Correction of few bugs in IncidenceMatrix, BlockTriangular and ModelTree 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2258 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2008-11-16 18:02:52 +00:00
parent af27f63392
commit 0b0f0704c4
4 changed files with 35 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
preprocessor/BlockTriangular.cc
View File

@ -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;j<symbol_table.exo_nbr;j++)
if(Cur_IM[i_1 + j] && (!tmp_exo[j]))
{
tmp_exo[j] = 1;
tmp_nb_exo++;
}
//Cur_IM = Cur_IM->pNext;
if(Cur_IM)
{
i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr;
for(j=0;j<symbol_table.exo_nbr;j++)
if(Cur_IM[i_1 + j] && (!tmp_exo[j]))
{
tmp_exo[j] = 1;
tmp_nb_exo++;
}
}
}
ModelBlock->Block_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)

4
preprocessor/InidenceMatrix.cc → preprocessor/IncidenceMatrix.cc
View File

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

1
preprocessor/Makefile
View File

@ -29,6 +29,7 @@ OBJS = \
ExprNode.o \
ModelNormalization.o \
ModelBlocks.o \
IncidenceMatrix.o \
BlockTriangular.o \
Model_Graph.o \
SymbolGaussElim.o \

49
preprocessor/ModelTree.cc
View File

@ -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;i<ModelBlock->Block_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;i<ModelBlock->Block_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<int> exogenous;//, tmp_exogenous(block_triangular.exo_nbr,-1);
vector<int> exogenous;
vector<int>::iterator it_exogenous;
exogenous.clear();
ostringstream tmp_s, tmp_s_eq;
@ -2935,7 +2900,6 @@ ModelTree::writeOutput(ostream &output) const
for(int i=0;i<block_triangular.ModelBlock->Block_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<symbol_table.endo_nbr*symbol_table.endo_nbr;i++)
{
@ -3119,7 +3083,6 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
int j=0;
bool *IM=NULL;
int a_variable_lag=-9999;
//block_triangular.Print_IM(2);
for(first_derivatives_type::iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
@ -3254,10 +3217,12 @@ ModelTree::computingPass(const eval_context_type &eval_context, bool no_tmp_term
if (mode == eSparseDLLMode || mode == eSparseMode)
{
BuildIncidenceMatrix();
jacob_map j_m;
evaluateJacobian(eval_context, &j_m);
if (block_triangular.bt_verbose)
{
cout << "The gross incidence matrix \n";