diff --git a/BlockTriangular.cc b/BlockTriangular.cc
index a87ce54bf8a63c4c65371c5a64f7ca72e032ed35..33633462535978d8ce03adf6ebfe1799e5d71063 100644
--- a/BlockTriangular.cc
+++ b/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;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)
diff --git a/InidenceMatrix.cc b/IncidenceMatrix.cc
similarity index 97%
rename from InidenceMatrix.cc
rename to IncidenceMatrix.cc
index 65ac872552550f7bfdccfecd7647e293c9a6cfcc..286eb1a26cd277a912489378f9bd37791401b04a 100644
--- a/InidenceMatrix.cc
+++ b/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/Makefile b/Makefile
index d1a85afe8d34454d550cad9599f858eaac87e3e8..1712c2ac9ee8f19e3b8d4dee735e6451279c78ef 100644
--- a/Makefile
+++ b/Makefile
@@ -29,6 +29,7 @@ OBJS = \
ExprNode.o \
ModelNormalization.o \
ModelBlocks.o \
+ IncidenceMatrix.o \
BlockTriangular.o \
Model_Graph.o \
SymbolGaussElim.o \
diff --git a/ModelTree.cc b/ModelTree.cc
index fdc79d687f21109e1c7ca2e95624caacffc47f16..ee7fd9db06cb518bfb86172b55ec3d88d0632b68 100644
--- a/ModelTree.cc
+++ b/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;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";