diff --git a/preprocessor/BlockTriangular.cc b/preprocessor/BlockTriangular.cc index bce845494..80631141a 100644 --- a/preprocessor/BlockTriangular.cc +++ b/preprocessor/BlockTriangular.cc @@ -88,17 +88,22 @@ BlockTriangular::init_incidence_matrix(int nb_endo) void -BlockTriangular::Free_IM(List_IM* First_IM) +BlockTriangular::Free_IM(List_IM* First_IM) const { +#ifdef DEBUG + cout << "Free_IM\n"; +#endif List_IM *Cur_IM, *SFirst_IM; Cur_IM = SFirst_IM = First_IM; while(Cur_IM) { First_IM = Cur_IM->pNext; free(Cur_IM->IM); + delete Cur_IM; Cur_IM = First_IM; } - free(SFirst_IM); + //free(SFirst_IM); + //delete SFirst_IM; } //------------------------------------------------------------------------------ @@ -223,7 +228,7 @@ BlockTriangular::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Inde SIM[pos2*n + j] = tmp_b; } } - /* ...and variables (colomn)*/ + /* ...and variables (column)*/ if(pos1 != pos3) { tmp_i = Index_Var_IM[pos1].index; @@ -241,7 +246,7 @@ BlockTriangular::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Inde //------------------------------------------------------------------------------ // Find the prologue and the epilogue of the model void -BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM) +BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0) { bool modifie = 1; int i, j, k, l = 0; @@ -261,7 +266,7 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n l = j; } } - if ((k == 1) /* && (l==i)*/) + if ((k == 1) && IM0[Index_Equ_IM[i].index*n + Index_Var_IM[l].index]) { modifie = 1; swap_IM_c(IM, *prologue, i, l, Index_Var_IM, Index_Equ_IM, n); @@ -273,6 +278,8 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n } *epilogue = 0; modifie = 1; + /*print_SIM(IM,n); + print_SIM(IM*/ while(modifie) { modifie = 0; @@ -287,7 +294,7 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n l = j; } } - if ((k == 1) /* && (l==i)*/) + if ((k == 1) && IM0[Index_Equ_IM[l].index*n + Index_Var_IM[i].index]) { modifie = 1; swap_IM_c(IM, n - (1 + *epilogue), l, i, Index_Var_IM, Index_Equ_IM, n); @@ -330,10 +337,10 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int while(Cur_IM) { k = Cur_IM->lead_lag; - i_1 = Index_Var_IM[*count_Equ].index * endo_nbr; + i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr; if(k > 0) { - if(Cur_IM->IM[i_1 + Index_Equ_IM[ /*j*/*count_Equ].index]) + if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index]) { nb_lead_lag_endo++; tmp_size[Model_Max_Lag + k]++; @@ -348,7 +355,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int else { k = -k; - if(Cur_IM->IM[i_1 + Index_Equ_IM[ /*j*/*count_Equ].index]) + if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index]) { tmp_size[Model_Max_Lag - k]++; nb_lead_lag_endo++; @@ -374,39 +381,45 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int ModelBlock->in_Block_Var[Index_Var_IM[*count_Equ].index] = *count_Block; ModelBlock->in_Equ_of_Block[Index_Equ_IM[*count_Equ].index] = ModelBlock->in_Var_of_Block[Index_Var_IM[*count_Equ].index] = 0; Index_Equ_IM[*count_Equ].block = *count_Block; - cout << "Lead=" << Lead << " Lag=" << Lag << "\n"; if ((Lead > 0) && (Lag > 0)) + ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE; + else if((Lead > 0) && (Lag == 0)) + ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_BACKWARD_SIMPLE; + else + ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FOREWARD_SIMPLE; + ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact)); + ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo; + if(nb_lead_lag_endo) { - ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE; - cout << "alloc ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag = (" << (Lead + Lag + 1) * sizeof(IM_compact) << ")\n"; - ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact)); - ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo; ModelBlock->Block_List[*count_Block].variable_dyn_index = (int*)malloc(nb_lead_lag_endo * sizeof(int)); ModelBlock->Block_List[*count_Block].variable_dyn_leadlag = (int*)malloc(nb_lead_lag_endo * sizeof(int)); - ls = l = 1; - i1 = 0; - for(i = 0;i < Lead + Lag + 1;i++) + } + ls = l = 1; + i1 = 0; + for(int li = 0;li < Lead + Lag + 1;li++) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = tmp_size[Model_Max_Lag - Lag + li]; + if(tmp_size[Model_Max_Lag - Lag + li]) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = tmp_size[Model_Max_Lag - Lag + i]; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_endo = tmp_size[Model_Max_Lag - Lag + i]; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l; - IM = bGet_IM(i - Lag); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_endo = tmp_size[Model_Max_Lag - Lag + li]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_dyn_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_init = l; + IM = bGet_IM(li - Lag); if(IM == NULL) { - cout << "Error IM(" << i - Lag << ") doesn't exist\n"; + cout << "Error IM(" << li - Lag << ") doesn't exist\n"; exit( -1); } - if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*endo_nbr]) + if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*endo_nbr] && nb_lead_lag_endo) { ModelBlock->Block_List[*count_Block].variable_dyn_index[i1] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = i - Lag; + ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = li - Lag; tmp_var[0] = i1; i1++; } @@ -414,31 +427,28 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr; if(IM[Index_Var_IM[*count_Equ].index + i_1]) { - if(i == Lag) + if(li == Lag) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us[m] = ls; ls++; } - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = 0; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = 0; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[0]]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u[m] = l; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ[m] = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var[m] = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index[m] = Index_Equ_IM[*count_Equ].index; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index[m] = Index_Var_IM[*count_Equ].index; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[0]]; l++; m++; } - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1; - } + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1; + } } - else if((Lead > 0) && (Lag == 0)) - ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_BACKWARD_SIMPLE; - else - ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FOREWARD_SIMPLE; (*count_Equ)++; (*count_Block)++; free(tmp_size); free(tmp_endo); + free(tmp_var); } } else @@ -535,7 +545,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int } ModelBlock->Block_List[*count_Block].Max_Lag = Lag; ModelBlock->Block_List[*count_Block].Max_Lead = Lead; - cout << "alloc ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag = (" << (Lead + Lag + 1) * sizeof(IM_compact) << ")\n"; ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact)); ls = l = size; i1 = 0; @@ -621,7 +630,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, int void -BlockTriangular::Free_Block(Model_Block* ModelBlock) +BlockTriangular::Free_Block(Model_Block* ModelBlock) const { int blk, i; #ifdef DEBUG @@ -637,12 +646,33 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) free(ModelBlock->Block_List[blk].Equation); free(ModelBlock->Block_List[blk].Variable); free(ModelBlock->Block_List[blk].Variable_Sorted); + free(ModelBlock->Block_List[blk].Own_Derivative); + if(ModelBlock->Block_List[blk].Nb_Lead_Lag_Endo) + { + free(ModelBlock->Block_List[blk].variable_dyn_index); + free(ModelBlock->Block_List[blk].variable_dyn_leadlag); + } + for(i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++) + { + free(ModelBlock->Block_List[blk].IM_lead_lag[i].u); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].us); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_dyn_Index); + } + free(ModelBlock->Block_List[blk].IM_lead_lag); + delete(ModelBlock->Block_List[blk].Temporary_terms); } else { free(ModelBlock->Block_List[blk].Equation); free(ModelBlock->Block_List[blk].Variable); free(ModelBlock->Block_List[blk].Variable_Sorted); + free(ModelBlock->Block_List[blk].Own_Derivative); + free(ModelBlock->Block_List[blk].variable_dyn_index); + free(ModelBlock->Block_List[blk].variable_dyn_leadlag); for(i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++) { free(ModelBlock->Block_List[blk].IM_lead_lag[i].u); @@ -651,8 +681,10 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var); free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index); free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_dyn_Index); } free(ModelBlock->Block_List[blk].IM_lead_lag); + delete(ModelBlock->Block_List[blk].Temporary_terms); } } free(ModelBlock->in_Block_Equ); @@ -661,6 +693,8 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) free(ModelBlock->in_Var_of_Block); free(ModelBlock->Block_List); free(ModelBlock); + free(Index_Equ_IM); + free(Index_Var_IM); } //------------------------------------------------------------------------------ @@ -672,10 +706,10 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int i, j, Nb_TotalBlocks, Nb_RecursBlocks; int count_Block, count_Equ; block_result_t* res; - List_IM * p_First_IM, *p_Cur_IM, *Cur_IM; + //List_IM * p_First_IM, *p_Cur_IM, *Cur_IM; Equation_set* Equation_gr = (Equation_set*) malloc(sizeof(Equation_set)); bool* SIM0, *SIM00; - p_First_IM = (List_IM*)malloc(sizeof(*p_First_IM)); + /*p_First_IM = (List_IM*)malloc(sizeof(*p_First_IM)); p_Cur_IM = p_First_IM; Cur_IM = First_IM; i = endo_nbr * endo_nbr; @@ -692,82 +726,110 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, } else p_Cur_IM->pNext = NULL; - } - Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM); + }*/ + SIM0 = (bool*)malloc(n * n * sizeof(bool)); + //cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n"; + memcpy(SIM0,IM_0,n*n*sizeof(bool)); + Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0); + //cout << "free SIM0=" << SIM0 << "\n"; + free(SIM0); if(bt_verbose) { cout << "prologue : " << *prologue << " epilogue : " << *epilogue << "\n"; + cout << "IM_0\n"; + Print_SIM(IM_0, n); + cout << "IM\n"; Print_SIM(IM, n); for(i = 0;i < n;i++) cout << "Index_Var_IM[" << i << "]=" << Index_Var_IM[i].index << " Index_Equ_IM[" << i << "]=" << Index_Equ_IM[i].index << "\n"; } - if(Do_Normalization) + if(*prologue+*epilogue, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) + cout << "Normalizing the model ...\n"; + if(mixing) { - if(fabs(iter->second)>max_val[iter->first.first]) - max_val[iter->first.first]=fabs(iter->second); - } - for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) - iter->second/=max_val[iter->first.first]; - free(max_val); - bool OK=false; - double bi=0.99999999; - int suppressed=0; - while(!OK && bi>1e-14) - { - int suppress=0; - SIM0 = (bool*)malloc(n * n * sizeof(bool)); - memset(SIM0,0,n*n*sizeof(bool)); - SIM00 = (bool*)malloc(n * n * sizeof(bool)); - memset(SIM00,0,n*n*sizeof(bool)); + double* max_val=(double*)malloc(n*sizeof(double)); + //cout << "n=" << n << "\n"; + memset(max_val,0,n*sizeof(double)); for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) - { - if(fabs(iter->second)>bi) - { - SIM0[iter->first.first*n+iter->first.second]=1; - if(!IM_0[iter->first.first*n+iter->first.second]) + { + //cout << "iter->first.first=" << iter->first.first << "\n"; + if(fabs(iter->second)>max_val[iter->first.first]) + max_val[iter->first.first]=fabs(iter->second); + } + for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) + iter->second/=max_val[iter->first.first]; + free(max_val); + bool OK=false; + double bi=0.99999999; + int suppressed=0; + while(!OK && bi>1e-14) + { + int suppress=0; + SIM0 = (bool*)malloc(n * n * sizeof(bool)); + //cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n"; + memset(SIM0,0,n*n*sizeof(bool)); + SIM00 = (bool*)malloc(n * n * sizeof(bool)); + //cout << "Allocate SIM00=" << SIM00 << " size=" << n * n * sizeof(bool) << "\n"; + memset(SIM00,0,n*n*sizeof(bool)); + //cout << "n*n=" << n*n << "\n"; + for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) + { + if(fabs(iter->second)>bi) + { + //cout << "iter->first.first*n+iter->first.second=" << iter->first.first*n+iter->first.second << "\n"; + SIM0[iter->first.first*n+iter->first.second]=1; + if(!IM_0[iter->first.first*n+iter->first.second]) + { + cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << "\n"; + } + } + else + suppress++; + } + //cout << "n*n=" << n*n << "\n"; + for(i = 0;i < n;i++) + for(j = 0;j < n;j++) { - cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << "\n"; + //cout << "Index_Equ_IM[i].index * n + Index_Var_IM[j].index=" << Index_Equ_IM[i].index * n + Index_Var_IM[j].index << "\n"; + SIM00[i*n + j] = SIM0[Index_Equ_IM[i].index * n + Index_Var_IM[j].index]; } - } - else - suppress++; - } - for(i = 0;i < n;i++) - for(j = 0;j < n;j++) - SIM00[i*n + j] = SIM0[Index_Equ_IM[i].index * n + Index_Var_IM[j].index]; - free(SIM0); - if(suppress!=suppressed) - OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM); - suppressed=suppress; + //cout << "free SIM0=" << SIM0 << "\n"; + free(SIM0); + if(suppress!=suppressed) + { + OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM); + if(!OK) + normalization.Free_Equation(n-*prologue-*epilogue,Equation_gr); + } + suppressed=suppress; + if(!OK) + bi/=1.07; + if(bi>1e-14) + free(SIM00); + } if(!OK) - bi/=1.07; - if(bi>1e-14) - free(SIM00); - } - if(!OK) - { - normalization.Set_fp_verbose(true); - OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM); - cout << "Error\n"; - exit(-1); + { + normalization.Set_fp_verbose(true); + OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM); + cout << "Error\n"; + exit(-1); + } } + else + normalization.Normalize(n, *prologue, *epilogue, IM, Index_Equ_IM, Equation_gr, 0, 0); } else - normalization.Normalize(n, *prologue, *epilogue, IM, Index_Equ_IM, Equation_gr, 0, 0); + normalization.Gr_to_IM_basic(n, *prologue, *epilogue, IM, Equation_gr, false); } - else - normalization.Gr_to_IM_basic(n, *prologue, *epilogue, IM, Equation_gr, false); cout << "Finding the optimal block decomposition of the model ...\n"; if(bt_verbose) blocks.Print_Equation_gr(Equation_gr); res = blocks.sc(Equation_gr); + normalization.Free_Equation(n-*prologue-*epilogue,Equation_gr); + free(Equation_gr); if(bt_verbose) blocks.block_result_print(res); if ((*prologue) || (*epilogue)) @@ -775,8 +837,10 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, else j = 0; for(i = 0;i < res->n_sets;i++) - if ((res->sets_f[i] - res->sets_s[i] + 1) > j) - j = res->sets_f[i] - res->sets_s[i] + 1; + { + if ((res->sets_f[i] - res->sets_s[i] + 1) > j) + j = res->sets_f[i] - res->sets_s[i] + 1; + } Nb_RecursBlocks = *prologue + *epilogue; Nb_TotalBlocks = res->n_sets + Nb_RecursBlocks; cout << Nb_TotalBlocks << " block(s) found:\n"; @@ -790,8 +854,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, ModelBlock->in_Var_of_Block = (int*)malloc(n * sizeof(int)); ModelBlock->Block_List = (Block*)malloc(sizeof(ModelBlock->Block_List[0]) * Nb_TotalBlocks); blocks.block_result_to_IM(res, IM, *prologue, endo_nbr, Index_Equ_IM, Index_Var_IM); - Free_IM(p_First_IM); + //Free_IM(p_First_IM); count_Equ = count_Block = 0; + //Print_IM(endo_nbr); if (*prologue) Allocate_Block(*prologue, &count_Equ, &count_Block, PROLOGUE, ModelBlock); for(j = 0;j < res->n_sets;j++) @@ -803,6 +868,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, } if (*epilogue) Allocate_Block(*epilogue, &count_Equ, &count_Block, EPILOGUE, ModelBlock); + blocks.block_result_free(res); return 0; } @@ -852,6 +918,8 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_ for(i = 0;i < endo_nbr*endo_nbr;i++) SIM_0[i] = Cur_IM->IM[i]; Normalize_and_BlockDecompose(SIM, ModelBlock, endo_nbr, &prologue, &epilogue, Index_Var_IM, Index_Equ_IM, 1, 1, SIM_0, j_m); + free(SIM_0); + free(SIM); if(bt_verbose) for(i = 0;i < endo_nbr;i++) cout << "Block=" << Index_Equ_IM[i].block << " Equ=" << Index_Equ_IM[i].index << " Var= " << Index_Var_IM[i].index << " " << symbol_table.getNameByID(eEndogenous, Index_Var_IM[i].index) << "\n"; diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index f139a0192..26f8fff99 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -68,6 +68,23 @@ CheckStatement::checkPass(ModFileStructure &mod_file_struct) mod_file_struct.check_present = true; } +Model_InfoStatement::Model_InfoStatement(const OptionsList &options_list_arg) : + options_list(options_list_arg) +{ +} + +void Model_InfoStatement::checkPass(ModFileStructure &mod_file_struct) +{ + //mod_file_struct.model_info_present = true; +} + +void Model_InfoStatement::writeOutput(ostream &output, const string &basename) const +{ + options_list.writeOutput(output); + output << "model_info();\n"; +} + + SimulStatement::SimulStatement(const OptionsList &options_list_arg) : options_list(options_list_arg) { diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 85c0c0865..2fee306cc 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -183,6 +183,7 @@ statement : declaration | dynatype | dynasave | model_comparison + | model_info | planner_objective | ramsey_policy | bvar_density @@ -612,6 +613,10 @@ check_options_list : check_options_list COMMA check_options check_options : o_solve_algo; +model_info : MODEL_INFO ';' + { driver.model_info(); } + ; + simul : SIMUL ';' { driver.simulate(); } | SIMUL '(' simul_options_list ')' ';' diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index e222fe5df..6a2ebd90f 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -104,6 +104,7 @@ int sigma_e = 0; periods {BEGIN DYNARE_STATEMENT; return token::PERIODS;} cutoff {BEGIN DYNARE_STATEMENT; return token::CUTOFF;} markowitz {BEGIN DYNARE_STATEMENT; return token::MARKOWITZ;} +model_info {BEGIN DYNARE_STATEMENT; return token::MODEL_INFO;} estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;} prior_analysis {BEGIN DYNARE_STATEMENT; return token::PRIOR_ANALYSIS;} posterior_analysis {BEGIN DYNARE_STATEMENT; return token::POSTERIOR_ANALYSIS;} diff --git a/preprocessor/DynareMain.cc b/preprocessor/DynareMain.cc index 61dc9bbb8..7a32a3a0a 100644 --- a/preprocessor/DynareMain.cc +++ b/preprocessor/DynareMain.cc @@ -86,6 +86,7 @@ main(int argc, char** argv) } // Do the rest + main2(macro_output, basename, debug, clear_all); return 0; diff --git a/preprocessor/DynareMain2.cc b/preprocessor/DynareMain2.cc index 53a5584d2..7061e71e2 100644 --- a/preprocessor/DynareMain2.cc +++ b/preprocessor/DynareMain2.cc @@ -28,10 +28,10 @@ void main2(stringstream &in, string &basename, bool debug, bool clear_all) { ParsingDriver p; - + //cout << "OK\n"; // Do parsing and construct internal representation of mod file ModFile *mod_file = p.parse(in, debug); - + //cout << "OK1\n"; // Run checking pass mod_file->checkPass(); diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index b09b789a2..5a9699e01 100644 --- a/preprocessor/ModFile.cc +++ b/preprocessor/ModFile.cc @@ -150,6 +150,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const << "warning warning_old_state" << endl; mOutputFile << "logname_ = '" << basename << ".log';" << endl; mOutputFile << "diary " << basename << ".log" << endl; + mOutputFile << "options_.model_mode = " << model_tree.mode << ";\n"; if (model_tree.equation_number() > 0) diff --git a/preprocessor/ModelBlocks.cc b/preprocessor/ModelBlocks.cc index f37117a2a..a7a6ea12b 100644 --- a/preprocessor/ModelBlocks.cc +++ b/preprocessor/ModelBlocks.cc @@ -207,6 +207,9 @@ Blocks::sc(Equation_set *g) void Blocks::block_result_free(block_result_t *r) { +#ifdef DEBUG + cout << "block_result_free\n"; +#endif free(r->vertices); free(r->sets_s); free(r->sets_f); diff --git a/preprocessor/ModelNormalization.cc b/preprocessor/ModelNormalization.cc index fbd8cdb05..75dfe62b5 100644 --- a/preprocessor/ModelNormalization.cc +++ b/preprocessor/ModelNormalization.cc @@ -283,9 +283,23 @@ void Normalization::Gr_to_IM_basic(int n0, int prologue, int epilogue, bool* IM, Equation_set *Equation, bool transpose) { int i, j, edges, n; - Edge *e1; + Edge *e1, *e2; n = n0 - prologue - epilogue; Equation->size = n; + if(Equation->Number) + { + for(i = 0;i < n;i++) + { + e1 = Equation->Number[i].First_Edge; + while(e1 != NULL) + { + e2 = e1->next; + free(e1); + e1 = e2; + } + } + free(Equation->Number); + } Equation->Number = (Equation_vertex*)malloc(n * sizeof(Equation_vertex)); edges = 0; if(transpose) @@ -382,8 +396,12 @@ Normalization::Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* In } free(SIM); free(Index_Equ_IM_tmp); + //cout << "mixing=" << mixing << "\n"; if(mixing) - Gr_to_IM_basic(n0, prologue, epilogue, IM, Equation, true); + { + //Free_Equation(n,Equation); + Gr_to_IM_basic(n0, prologue, epilogue, IM, Equation, true); + } else { // In this step we : @@ -453,17 +471,21 @@ Normalization::Free_Equation(int n, Equation_set* Equation) while(e1 != NULL) { e2 = e1->next; + free(e1); e1 = e2; } } free(Equation->Number); - free(Equation); + //free(Equation); } void Normalization::Free_Other(Variable_set* Variable) { //free unused space +#ifdef DEBUG + cout << "Free_Other\n"; +#endif free(Local_Heap); free(Variable->Number); free(Variable); @@ -526,6 +548,7 @@ Normalization::Normalize(int n, int prologue, int epilogue, bool* IM, simple* In } } Free_Other(Variable); + //Free_All(n,Equation,Variable); #ifdef DEBUG cout << "end of Normalize\n"; #endif diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index a50e0932d..efd3661e2 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -335,7 +335,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) if (order == 2) for(second_derivatives_type::iterator it = second_derivatives.begin(); it != second_derivatives.end(); it++) - it->second->computeTemporaryTerms(reference_count, temporary_terms, false); + it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, 0, ModelBlock, map_idx); /*New*/ j=0; for(temporary_terms_type::const_iterator it = temporary_terms.begin(); @@ -391,7 +391,7 @@ ModelTree::writeModelEquationsOrdered_C(ostream &output, Model_Block *ModelBlock } else lhs_rhs_done=false; - if (prev_Simulation_Type==ModelBlock->Block_List[j].Simulation_Type + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) && (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R @@ -593,7 +593,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock bool OK, lhs_rhs_done, skip_the_head; ostringstream Uf[symbol_table.endo_nbr]; map reference_count; - int prev_Simulation_Type=-1; + int prev_Simulation_Type=-1, count_derivates=0; temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); //---------------------------------------------------------------------- //Temporary variables declaration @@ -610,12 +610,9 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock /*tmp_output << "[" << block_triangular.periods + variable_table.max_lag+variable_table.max_lead << "]";*/ } - - if (tmp_output.str().length()>0) - { - global_output << " global " << tmp_output.str() << " M_ ;\n"; - } + global_output << " global " << tmp_output.str() << " M_ ;\n"; //For each block + int gen_blocks=0; for(j = 0;j < ModelBlock->Size;j++) { //For a block composed of a single equation determines wether we have to evaluate or to solve the equation @@ -630,7 +627,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock } else lhs_rhs_done=false; - if (prev_Simulation_Type==ModelBlock->Block_List[j].Simulation_Type + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) && (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R @@ -640,6 +637,8 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock skip_the_head=false; if (!skip_the_head) { + count_derivates=0; + gen_blocks++; if (j>0) { output << "return;\n\n\n"; @@ -650,12 +649,12 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R) - output << "function [y] = " << dynamic_basename << "_" << j+1 << "(y, x, it_)\n"; + output << "function [y, g1, g2, g3] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, jacobian_eval, g1, g2, g3)\n"; else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE ||ModelBlock->Block_List[j].Simulation_Type==SOLVE_FOREWARD_SIMPLE) - output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_)\n"; + output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, g1, g2, g3, y_index, jacobian_eval)\n"; else - output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, y_kmin, y_size, periods)\n"; + output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, y_kmin, y_size, periods, g1, g2, g3)\n"; output << " % ////////////////////////////////////////////////////////////////////////" << endl << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << " //" << endl @@ -676,8 +675,8 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock int nze; for(nze=0,m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) nze+=ModelBlock->Block_List[j].IM_lead_lag[m].size; - output << " Jacobian_Size=" << ModelBlock->Block_List[j].Size << "*(y_kmin+" << ModelBlock->Block_List[j].Max_Lead << " +periods);\n"; - output << " g1=spalloc( y_size*periods, Jacobian_Size, " << nze << "*periods" << ");\n"; + //output << " Jacobian_Size=" << ModelBlock->Block_List[j].Size << "*(y_kmin+" << ModelBlock->Block_List[j].Max_Lead << " +periods);\n"; + //output << " g1=spalloc( y_size*periods, Jacobian_Size, " << nze << "*periods" << ");\n"; output << " for it_ = y_kmin+1:(periods+y_kmin)\n"; output << " Per_y_=it_*y_size;\n"; output << " Per_J_=(it_-y_kmin-1)*y_size;\n"; @@ -706,8 +705,8 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock { ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ; - output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i] << " variable : " << sModel - << " (" << ModelBlock->Block_List[j].Variable[i] << ")" << endl; + output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel + << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl; if (!lhs_rhs_done) { eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; @@ -762,101 +761,142 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R - && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD_R) - { + && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD_R + && ModelBlock->Block_List[j].Simulation_Type!=SOLVE_FOREWARD_SIMPLE + && ModelBlock->Block_List[j].Simulation_Type!=SOLVE_BACKWARD_SIMPLE) output << " " << sps << "% Jacobian " << endl; - switch(ModelBlock->Block_List[j].Simulation_Type) + else + { + output << " " << sps << "% Jacobian " << endl << " if jacobian_eval" << endl; + } + switch(ModelBlock->Block_List[j].Simulation_Type) + { + case SOLVE_BACKWARD_SIMPLE: + case SOLVE_FOREWARD_SIMPLE: + case EVALUATE_BACKWARD: + case EVALUATE_FOREWARD: + case EVALUATE_BACKWARD_R: + case EVALUATE_FOREWARD_R: + count_derivates++; + for(m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++) { - case SOLVE_BACKWARD_SIMPLE: - case SOLVE_FOREWARD_SIMPLE: - output << " g1(1)="; + k=m-ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) + { + if(ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]==ModelBlock->Block_List[j].Variable[0]) + { + output << " g1(M_.block_structure.block(" << gen_blocks << ").equation(" << count_derivates << "), M_.block_structure.block(" << gen_blocks << ").variable(" << count_derivates << ")+" << (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")="; + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], k, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) + << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) + << ") " << ModelBlock->Block_List[j].Variable[0]+1 + << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; + } + } + } + if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE + || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FOREWARD_SIMPLE) + { + output << " else\n"; + output << " g1="; writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabDynamicModelSparse, temporary_terms); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) - << ") " << ModelBlock->Block_List[j].Variable[0] - << ", equation=" << ModelBlock->Block_List[j].Equation[0] << endl; - break; - case SOLVE_BACKWARD_COMPLETE: - case SOLVE_FOREWARD_COMPLETE: - m=ModelBlock->Block_List[j].Max_Lag; + << ") " << ModelBlock->Block_List[j].Variable[0]+1 + << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; + } + if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD + || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD + || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R + || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R + || ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE + || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FOREWARD_SIMPLE) + output << " end;" << endl; + break; + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FOREWARD_COMPLETE: + m=ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) + { + 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 u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u(" << u << ")*y(Per_y_+" << var << ")"; + output << " u(" << u+1 << ") = "; + writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) + << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; + break; + case SOLVE_TWO_BOUNDARIES_COMPLETE: + output << " g2=0;g3=0;\n"; + for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) { 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 u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i]; + //int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u(" << u << ")*y(Per_y_+" << var << ")"; - output << " u(" << u+1 << ") = "; - writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms); + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; + if (k==0) + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_)*y(it_, " << var+1 << ")"; + else if (k==1) + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_)*y(it_+1, " << var+1 << ")"; + else if (k>0) + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")"; + else if (k<0) + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")"; + //output << " u(" << u+1 << "+Per_u_) = "; + if(k==0) + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; + else if(k==1) + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = "; + else if(k>0) + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = "; + else if(k<0) + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var - << ", equation=" << eq << endl; - } - for(i = 0;i < ModelBlock->Block_List[j].Size;i++) - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; - break; - case SOLVE_TWO_BOUNDARIES_COMPLETE: - output << " g2=0;g3=0;\n"; - for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=m-ModelBlock->Block_List[j].Max_Lag; - for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) - { - 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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; - if (k==0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_)*y(it_, " << var+1 << ")"; - else if (k>0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k << "-1))*y(it_+" << k << ", " << var+1 << ")"; - else if (k<0) - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k << "-1))*y(it_" << k << ", " << var+1 << ")"; - //output << " u(" << u+1 << "+Per_u_) = "; - if(k==0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; - else if(k>0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k << "-1)) = "; - else if(k<0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k << "-1)) = "; - writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << k << ") " << var - << ", equation=" << eq << endl; + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; #ifdef CONDITION - output << " if (fabs(condition[" << eqr << "])Block_List[j].Size;i++) - { - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; -#ifdef CONDITION - output << " if (fabs(condition(" << i+1 << "))Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=m-ModelBlock->Block_List[j].Max_Lag; - for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) - { - 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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n"; - } - } - for(i = 0;i < ModelBlock->Block_List[j].Size;i++) - output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; -#endif - output << " end;\n"; - break; } + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; +#ifdef CONDITION + output << " if (fabs(condition(" << i+1 << "))Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) + { + 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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n"; + } + } + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; +#endif + output << " end;\n"; + break; } prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; } @@ -890,10 +930,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode tmp_output << " "; (*it)->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms); } - if (tmp_output.str().length()>0) - { - global_output << " global " << tmp_output.str() << " M_ ;\n"; - } + global_output << " global " << tmp_output.str() << " M_ ;\n"; //For each block for(j = 0;j < ModelBlock->Size;j++) { @@ -909,7 +946,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode } else lhs_rhs_done=false; - if (prev_Simulation_Type==ModelBlock->Block_List[j].Simulation_Type + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) && (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R @@ -925,7 +962,13 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode } else output << "\n\n"; - output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x)\n"; + if(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD + ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD + ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R + ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R ) + output << "function [y] = " << static_basename << "_" << j+1 << "(y, x)\n"; + else + output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x)\n"; output << " % ////////////////////////////////////////////////////////////////////////" << endl << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << " //" << endl @@ -939,43 +982,36 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode output << " end\n"; } - temporary_terms_type tt2; int n=ModelBlock->Block_List[j].Size; int n1=symbol_table.endo_nbr; - //cout << "n1=" << n1 << "\n"; IM=(bool*)malloc(n*n*sizeof(bool)); memset(IM, 0, n*n*sizeof(bool)); - //cout << "ModelBlock->Block_List[j].Max_Lead" << ModelBlock->Block_List[j].Max_Lag << " ModelBlock->Block_List[j].Max_Lag=" << ModelBlock->Block_List[j].Max_Lead << "\n"; for(m=-ModelBlock->Block_List[j].Max_Lag;m<=ModelBlock->Block_List[j].Max_Lead;m++) { - //cout << "bGet_IM(" << m << ")\n"; IMl=block_triangular.bGet_IM(m); - //cout <<"OK\n"; for(i=0;iBlock_List[j].Equation[i]; for(k=0;kBlock_List[j].Variable[k]; - //cout << "eq=" << eq << " var=" << var << "\n"; IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var]; - /*if(IM[i*n+k]) - cout << " ->i=" << i << " j=" << j << "\n";*/ } } - //cout << "done\n"; } for(nze=0, i=0;iBlock_List[j].Simulation_Type!=EVALUATE_BACKWARD && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD_R) - output << " g1=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n"; + { + output << " g1=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n"; + output << " residual=zeros(" << ModelBlock->Block_List[j].Size << ",1);\n"; + } sps=""; if (ModelBlock->Block_List[j].Temporary_terms->size()) output << " " << sps << "% //Temporary variables" << endl; @@ -997,8 +1033,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode { ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ; - output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i] << " variable : " - << sModel << " (" << ModelBlock->Block_List[j].Variable[i] << ")" << endl; + output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " + << sModel << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl; if (!lhs_rhs_done) { eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; @@ -1057,8 +1093,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) - << ") " << ModelBlock->Block_List[j].Variable[0] - << ", equation=" << ModelBlock->Block_List[j].Equation[0] << endl; + << ") " << ModelBlock->Block_List[j].Variable[0]+1 + << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; break; case SOLVE_BACKWARD_COMPLETE: case SOLVE_FOREWARD_COMPLETE: @@ -1073,8 +1109,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode output << " u(" << u+1 << ") = "; writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var - << ", equation=" << eq << endl; + << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 + << ", equation=" << eq+1 << endl; } for(i = 0;i < ModelBlock->Block_List[j].Size;i++) output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; @@ -1091,8 +1127,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode //int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; - output << "% i=" << i << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << "\n"; - cout << "% i=" << i << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << "\n"; + //output << "% i=" << i << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << "\n"; if(!IM[eqr*ModelBlock->Block_List[j].Size+varr]) { Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 @@ -1102,8 +1137,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode output << " g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + "; writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << k << ") " << var - << ", equation=" << eq << endl; + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; #ifdef CONDITION output << " if (fabs(condition[" << eqr << "])Block_List[j].Simulation_Type; + free(IM); } output << "return;\n\n\n"; - free(IM); + //free(IM); } @@ -1190,7 +1226,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl ModelBlock_Aggregated_Count=-1; for(j = 0;j < ModelBlock->Size;j++) { - if (prev_Simulation_Type==ModelBlock->Block_List[j].Simulation_Type + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) && (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R @@ -1266,10 +1302,10 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl } else lhs_rhs_done=false; - if (ModelBlock->Block_List[j].Size==1) + /*if (ModelBlock->Block_List[j].Size==1) lhs_rhs_done=true; else - lhs_rhs_done=false; + lhs_rhs_done=false;*/ //The Temporary terms temporary_terms_type tt2; i=0; @@ -1994,7 +2030,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); SaveCode.write(reinterpret_cast(&u), sizeof(u)); - cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " u=" << u << "\n"; + //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " u=" << u << "\n"; u_count_int++; } } @@ -2008,7 +2044,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); - cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " eqr1=" << eqr1 << "\n"; + //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " eqr1=" << eqr1 << "\n"; u_count_int++; } for(j=0;jBlock_List[num].Size;j++) @@ -2052,9 +2088,10 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b mStaticModelFile << " if(length(varargin)>0)\n"; mStaticModelFile << " %it is a simple evaluation of the dynamic model for time _it\n"; mStaticModelFile << " global it_;\n"; - mStaticModelFile << " y=varargin{1}(y_kmin,:);\n"; + mStaticModelFile << " y=varargin{1}(:);\n"; mStaticModelFile << " ys=y;\n"; - mStaticModelFile << " x=varargin{2}(y_kmin,:);\n"; + mStaticModelFile << " g1=[];\n"; + mStaticModelFile << " x=varargin{2}(:);\n"; mStaticModelFile << " residual=zeros(1, " << symbol_table.endo_nbr << ");\n"; prev_Simulation_Type=-1; for(i=0;iSize;i++) @@ -2066,7 +2103,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b } mStaticModelFile << " ];\n"; k=block_triangular.ModelBlock->Block_List[i].Simulation_Type; - if (prev_Simulation_Type==k && + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)) skip_head=true; else @@ -2078,7 +2115,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b case EVALUATE_FOREWARD_R: case EVALUATE_BACKWARD_R: if(!skip_head) - mStaticModelFile << " " << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n"; mStaticModelFile << " residual(y_index)=ys(y_index)-y(y_index);\n"; break; case SOLVE_FOREWARD_COMPLETE: @@ -2092,168 +2129,152 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b } prev_Simulation_Type=k; } - mStaticModelFile << " varargout{1}=residual;\n"; - mStaticModelFile << " varargout{2}=g1;\n"; - mStaticModelFile << " return;\n"; - mStaticModelFile << " end;\n"; - mStaticModelFile << " %it is the deterministic simulation of the block decomposed static model\n"; - mStaticModelFile << " periods=options_.periods;\n"; - mStaticModelFile << " maxit_=options_.maxit_;\n"; - mStaticModelFile << " solve_tolf=options_.solve_tolf;\n"; - mStaticModelFile << " y=oo_.steady_state;\n"; - mStaticModelFile << " x=oo_.exo_steady_state;\n"; - prev_Simulation_Type=-1; - for(i = 0;i < block_triangular.ModelBlock->Size;i++) - { - k = block_triangular.ModelBlock->Block_List[i].Simulation_Type; - if (prev_Simulation_Type==k && - (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)) - skip_head=true; - else - skip_head=false; - if ((k == EVALUATE_FOREWARD || k == EVALUATE_FOREWARD_R || k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (!skip_head) - { - if (open_par) - { - mStaticModelFile << " end\n"; - } - mStaticModelFile << " " << static_basename << "_" << i + 1 << "(y, x);\n"; - } - open_par=false; - } + mStaticModelFile << " varargout{1}=residual;\n"; + mStaticModelFile << " varargout{2}=g1;\n"; + mStaticModelFile << " return;\n"; + mStaticModelFile << " end;\n"; + mStaticModelFile << " %it is the deterministic simulation of the block decomposed static model\n"; + mStaticModelFile << " periods=options_.periods;\n"; + mStaticModelFile << " maxit_=options_.maxit_;\n"; + mStaticModelFile << " solve_tolf=options_.solve_tolf;\n"; + mStaticModelFile << " y=oo_.steady_state;\n"; + mStaticModelFile << " x=oo_.exo_steady_state;\n"; + mStaticModelFile << " varargout{2}=1;\n"; + prev_Simulation_Type=-1; + for(i = 0;i < block_triangular.ModelBlock->Size;i++) + { + k = block_triangular.ModelBlock->Block_List[i].Simulation_Type; + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && + (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)) + skip_head=true; + else + skip_head=false; + if ((k == EVALUATE_FOREWARD || k == EVALUATE_FOREWARD_R || k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size)) + { + if (!skip_head) + { + if (open_par) + { + mStaticModelFile << " end\n"; + } + mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n"; + } + open_par=false; + } else if ((k == SOLVE_FOREWARD_SIMPLE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (open_par) - { + { + if (open_par) + { mStaticModelFile << " end\n"; - } - open_par=false; - mStaticModelFile << " g1=0;\n"; - mStaticModelFile << " r=0;\n"; - /*mStaticModelFile << " for it_=y_kmin+1:periods+y_kmin\n"; - mStaticModelFile << " cvg=0;\n"; - mStaticModelFile << " iter=0;\n"; - mStaticModelFile << " Per_y_=it_*y_size;\n"; - mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n"; - mStaticModelFile << " y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] = y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]-r/g1;\n"; - mStaticModelFile << " cvg=((r[0]*r[0])maxit_),\n"; - mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x);\n"; - mStaticModelFile << " y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; - mStaticModelFile << " cvg=((r*r)Block_List[i].Size)) - { - if (open_par) - { - mStaticModelFile << "end\n"; - } - open_par=false; - mStaticModelFile << " y_index=["; - for(int ik=0;ikBlock_List[i].Size;ik++) - { - mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; - } - mStaticModelFile << " ];\n"; - mStaticModelFile << " g1=0;g2=0;g3=0;\n"; - mStaticModelFile << " r=0;\n"; - /*mStaticModelFile << " for it_=y_kmin+1:periods+y_kmin\n"; - mStaticModelFile << " cvg=0;\n"; - mStaticModelFile << " iter=0;\n"; - mStaticModelFile << " Per_y_=it_*y_size;\n"; - mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n"; - mStaticModelFile << " [L, U] = LU(g1);\n"; - mStaticModelFile << " y(it_, y_index) = U\\(L\\b);\n"; - mStaticModelFile << " cvg=((r'*r)maxit_),\n"; - mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x);\n"; - mStaticModelFile << " max_res=max(abs(r));\n"; - mStaticModelFile << " if(iter>0)\n"; - mStaticModelFile << " if(~isreal(max_res) | max_resa1e-6)\n"; - mStaticModelFile << " lambda=lambda/2;\n"; - mStaticModelFile << " y(y_index)=y_save+lambda*dx;\n"; - mStaticModelFile << " continue;\n"; - mStaticModelFile << " else\n"; - mStaticModelFile << " disp(['No convergence after ' num2str(iter,'%d') ' iterations']);\n"; - mStaticModelFile << " return;\n"; - mStaticModelFile << " end;\n"; - mStaticModelFile << " else\n"; - mStaticModelFile << " if(lambda<1)\n"; - mStaticModelFile << " lambda=max(lambda*2, 1);\n"; - mStaticModelFile << " end;\n"; - mStaticModelFile << " end;\n"; - mStaticModelFile << " end;\n"; - mStaticModelFile << " max_resa=max_res;\n"; - mStaticModelFile << " cvg=(max_resmaxit_),\n"; + mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n"; + mStaticModelFile << " y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] = y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]-r/g1;\n"; + mStaticModelFile << " cvg=((r[0]*r[0])maxit_),\n"; + mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; + mStaticModelFile << " cvg=((r*r)Block_List[i].Size)) + { + if (open_par) + { + mStaticModelFile << "end\n"; + } + open_par=false; + mStaticModelFile << " y_index=["; + for(int ik=0;ikBlock_List[i].Size;ik++) + { + mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + } + mStaticModelFile << " ];\n"; + mStaticModelFile << " g1=0;g2=0;g3=0;\n"; + mStaticModelFile << " r=0;\n"; + /*mStaticModelFile << " for it_=y_kmin+1:periods+y_kmin\n"; + mStaticModelFile << " cvg=0;\n"; + mStaticModelFile << " iter=0;\n"; + mStaticModelFile << " Per_y_=it_*y_size;\n"; + mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n"; + mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n"; + mStaticModelFile << " [L, U] = LU(g1);\n"; + mStaticModelFile << " y(it_, y_index) = U\\(L\\b);\n"; + mStaticModelFile << " cvg=((r'*r)maxit_),\n"; + mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " max_res=max(abs(r));\n"; + mStaticModelFile << " cvg=(max_res0)\n"; mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n"; - mDynamicModelFile << " global 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"; mDynamicModelFile << " Per_u_=0;\n"; mDynamicModelFile << " Per_y_=it_*y_size;\n"; - mDynamicModelFile << " y1=varargin{1};\n"; + mDynamicModelFile << " y=varargin{1};\n"; + mDynamicModelFile << " ys=y(it_,:);\n"; + /*mDynamicModelFile << " y1=varargin{1};\n"; mDynamicModelFile << " cnb_nz_elem=1;\n"; mDynamicModelFile << " for i = -y_kmin:y_kmax\n"; mDynamicModelFile << " nz_elem=find(M_.lead_lag_incidence(:,1+i+y_kmin));\n"; @@ -2397,19 +2425,47 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons mDynamicModelFile << " nz_elem_s=nz_elem;\n"; mDynamicModelFile << " end;\n"; mDynamicModelFile << " cnb_nz_elem=cnb_nz_elem+nb_nz_elem;\n"; - mDynamicModelFile << " end;\n"; + mDynamicModelFile << " end;\n";*/ mDynamicModelFile << " x=varargin{2};\n"; prev_Simulation_Type=-1; + tmp.str(""); + tmp_eq.str(""); for(i = 0;i < block_triangular.ModelBlock->Size;i++) { - mDynamicModelFile << " y_index=["; + + k=block_triangular.ModelBlock->Block_List[i].Simulation_Type; + if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) && + ((prev_Simulation_Type==EVALUATE_FOREWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FOREWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R) + || (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))) + { + mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n"; + tmp_eq.str(""); + mDynamicModelFile << " y_index=[" << tmp.str() << "];\n"; + tmp.str(""); + mDynamicModelFile << tmp1.str(); + tmp1.str(""); + } + //mDynamicModelFile << " y_index=["; for(int ik=0;ikBlock_List[i].Size;ik++) { - mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; + tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; } - mDynamicModelFile << " ];\n"; - k=block_triangular.ModelBlock->Block_List[i].Simulation_Type; - if (prev_Simulation_Type==k && + //mDynamicModelFile << " ];\n"; + if(k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R + ) + { + if(i==block_triangular.ModelBlock->Size-1) + { + mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n"; + tmp_eq.str(""); + mDynamicModelFile << " y_index=[" << tmp.str() << "];\n"; + tmp.str(""); + mDynamicModelFile << tmp1.str(); + tmp1.str(""); + } + } + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)) skip_head=true; else @@ -2421,19 +2477,49 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons case EVALUATE_FOREWARD_R: case EVALUATE_BACKWARD_R: if(!skip_head) - mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_kmin, Per_u_, Per_y_, y_size);\n"; - mDynamicModelFile << " residual(y_index)=ys(y_index)-y(it_, y_index);\n"; + { + tmp1 << " [y, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n"; + tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n"; + } + break; + case SOLVE_FOREWARD_SIMPLE: + case SOLVE_BACKWARD_SIMPLE: + mDynamicModelFile << " y_index=" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ";\n"; + mDynamicModelFile << " [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, y_index, 1);\n"; + mDynamicModelFile << " residual(y_index_eq)=r;\n"; break; case SOLVE_FOREWARD_COMPLETE: case SOLVE_BACKWARD_COMPLETE: case SOLVE_TWO_BOUNDARIES_COMPLETE: - mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_size, it_);\n"; - mDynamicModelFile << " residual(y_index)=r;\n"; + //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_size, 1);\n"; + mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; + mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; + tmp.str(""); + tmp_eq.str(""); + int tmp_i=variable_table.max_lag+variable_table.max_lead+1; + mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n"; + mDynamicModelFile << " y_index_c = y_index;\n"; + tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1; + mDynamicModelFile << " for i=1:" << tmp_i-1 << ",\n"; + mDynamicModelFile << " y_index_c = [y_index_c (y_index+i*y_size)];\n"; + mDynamicModelFile << " end;\n"; + mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n"; + if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead) + mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga;\n"; + else + mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n"; + mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n"; break; } prev_Simulation_Type=k; } - mDynamicModelFile << " varagout{1}=residual;\n"; + if(tmp1.str().length()) + { + mDynamicModelFile << tmp1.str(); + tmp1.str(""); + } + mDynamicModelFile << " varargout{1}=residual;\n"; + mDynamicModelFile << " varargout{2}=g1;\n"; mDynamicModelFile << " return;\n"; mDynamicModelFile << " end;\n"; mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n"; @@ -2459,7 +2545,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons for(i = 0;i < block_triangular.ModelBlock->Size;i++) { k = block_triangular.ModelBlock->Block_List[i].Simulation_Type; - if (prev_Simulation_Type==k && + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)) skip_head=true; else @@ -2490,7 +2576,8 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons mDynamicModelFile << " Per_u_=0;\n"; mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n"; mDynamicModelFile << " Per_y_=it_*y_size;\n"; - mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, it_);\n"; + mDynamicModelFile << " g1=[];g2=[];g3=[];\n"; + mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n"; } } if(mode==eSparseDLLMode) @@ -2524,7 +2611,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons mDynamicModelFile << " Per_u_=0;\n"; mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n"; mDynamicModelFile << " Per_y_=it_*y_size;\n"; - mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, it_);\n"; + mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3);\n"; } } if(mode==eSparseDLLMode) @@ -2585,7 +2672,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons mDynamicModelFile << " iter=0;\n"; mDynamicModelFile << " Per_y_=it_*y_size;\n"; mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_);\n"; + mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n"; mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; mDynamicModelFile << " cvg=((r*r)maxit_),\n"; - mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_);\n"; - mDynamicModelFile << " y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] = y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]-r[it_]/g1;\n"; - mDynamicModelFile << " cvg=((r[it_]*r[it_])Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; + mDynamicModelFile << " cvg=((r*r)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 << ";\n"; mDynamicModelFile << " Per_y_=it_*y_size;\n"; - mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, r, g1, g2);\n"; + mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, r, g1, g2, g3);\n"; mDynamicModelFile << "#ifdef PRINT_OUT\n"; mDynamicModelFile << " for(j=0;j<" << 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 << ";j++)\n"; mDynamicModelFile << " {\n"; @@ -2800,7 +2887,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons } else { - mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n"; + mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2, g3);\n"; mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", 0, false);\n"; } mDynamicModelFile << " }\n"; @@ -2997,6 +3084,19 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons mDynamicModelFile << " y_kmax_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lead << ";\n"; mDynamicModelFile << " lambda=options_.slowc;\n"; mDynamicModelFile << " correcting_factor=0.01;\n"; + mDynamicModelFile << " luinc_tol=1e-10;\n"; + mDynamicModelFile << " max_resa=1e100;\n"; + 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++) + nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; + mDynamicModelFile << " Jacobian_Size=" << block_triangular.ModelBlock->Block_List[i].Size << "*(y_kmin+" << block_triangular.ModelBlock->Block_List[i].Max_Lead << " +periods);\n"; + mDynamicModelFile << " g1=spalloc( length(y_index)*periods, Jacobian_Size, " << nze << "*periods" << ");\n"; + mDynamicModelFile << " cpath=path;\n"; + mDynamicModelFile << " addpath(fullfile(matlabroot,'toolbox','matlab','sparfun'));\n"; + mDynamicModelFile << " bicgstabh=@bicgstab;\n"; + mDynamicModelFile << " path(cpath);\n"; + mDynamicModelFile << sp << " reduced = 0;\n"; + //mDynamicModelFile << " functions(bicgstabh)\n"; if (!block_triangular.ModelBlock->Block_List[i].is_linear) { sp=" "; @@ -3006,7 +3106,7 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons { sp=""; } - mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, y_kmin, Blck_size, periods);\n"; + mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, y_kmin, Blck_size, periods, g1, g2, g3);\n"; mDynamicModelFile << sp << " g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);\n"; mDynamicModelFile << sp << " b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';\n"; mDynamicModelFile << sp << " if(~isreal(r))\n"; @@ -3014,9 +3114,8 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons mDynamicModelFile << sp << " else\n"; mDynamicModelFile << sp << " max_res=max(max(abs(r)));\n"; mDynamicModelFile << sp << " end;\n"; - mDynamicModelFile << sp << " if(iter>0)\n"; - mDynamicModelFile << sp << " if(~isreal(max_res) | isnan(max_res) | max_resa1))\n"; mDynamicModelFile << sp << " if(isnan(max_res))\n"; mDynamicModelFile << sp << " detJ=det(g1aa);\n"; mDynamicModelFile << sp << " if(abs(detJ)<1e-7)\n"; @@ -3037,8 +3136,9 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons mDynamicModelFile << sp << " return;\n"; mDynamicModelFile << sp << " end;\n"; mDynamicModelFile << sp << " end;\n"; - mDynamicModelFile << sp << " elseif(lambda>1e-6)\n"; + mDynamicModelFile << sp << " elseif(lambda>1e-8)\n"; mDynamicModelFile << sp << " lambda=lambda/2;\n"; + mDynamicModelFile << sp << " reduced = 1;\n"; mDynamicModelFile << sp << " disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);\n"; mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';\n"; if (!block_triangular.ModelBlock->Block_List[i].is_linear) @@ -3059,41 +3159,50 @@ ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, cons mDynamicModelFile << sp << " g1aa=g1a;\n"; mDynamicModelFile << sp << " ba=b;\n"; mDynamicModelFile << sp << " max_resa=max_res;\n"; - - mDynamicModelFile << sp << " if(options_.simulation_method==0),\n"; mDynamicModelFile << sp << " dx = g1a\\b- ya;\n"; mDynamicModelFile << sp << " ya = ya + lambda*dx;\n"; mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n"; mDynamicModelFile << sp << " elseif(options_.simulation_method==2),\n"; - mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,1e-6);\n"; - mDynamicModelFile << sp << " [za,flag1] = gmres(g1a,b," << block_triangular.ModelBlock->Block_List[i].Size << ",1e-6," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n"; - mDynamicModelFile << sp << " dx = za - ya;\n"; - mDynamicModelFile << sp << " ya = ya + lambda*dx;\n"; - mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n"; - mDynamicModelFile << sp << " if (flag1>0)\n"; - mDynamicModelFile << sp << " if(flag1==1)\n"; - mDynamicModelFile << sp << " disp(['No convergence inside GMRES after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n"; - mDynamicModelFile << sp << " elseif(flag1==2)\n"; - mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n"; - mDynamicModelFile << sp << " elseif(flag1==3)\n"; - mDynamicModelFile << sp << " disp(['GMRES stagnated. (Two consecutive iterates were the same.)']);\n"; + mDynamicModelFile << sp << " flag1=1;\n"; + mDynamicModelFile << sp << " while(flag1>0)\n"; + mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n"; + mDynamicModelFile << sp << " [za,flag1] = gmres(g1a,b," << block_triangular.ModelBlock->Block_List[i].Size << ",1e-6," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n"; + mDynamicModelFile << sp << " if (flag1>0 | reduced)\n"; + mDynamicModelFile << sp << " if(flag1==1)\n"; + mDynamicModelFile << sp << " disp(['No convergence inside GMRES after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n"; + mDynamicModelFile << sp << " elseif(flag1==2)\n"; + mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n"; + mDynamicModelFile << sp << " elseif(flag1==3)\n"; + mDynamicModelFile << sp << " disp(['GMRES stagnated. (Two consecutive iterates were the same.)']);\n"; + mDynamicModelFile << sp << " end;\n"; + mDynamicModelFile << sp << " luinc_tol = luinc_tol/10;\n"; + mDynamicModelFile << sp << " reduced = 0;\n"; + mDynamicModelFile << sp << " else\n"; + mDynamicModelFile << sp << " dx = za - ya;\n"; + mDynamicModelFile << sp << " ya = ya + lambda*dx;\n"; + mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n"; mDynamicModelFile << sp << " end;\n"; mDynamicModelFile << sp << " end;\n"; mDynamicModelFile << sp << " elseif(options_.simulation_method==3),\n"; - mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,1e-7);\n"; - mDynamicModelFile << sp << " [za,flag1] = bicgstab(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n"; - mDynamicModelFile << sp << " dx = za - ya;\n"; - mDynamicModelFile << sp << " ya = ya + lambda*dx;\n"; - //mDynamicModelFile << sp << " [ya,flag1] = eval(strcat(fullfile(matlabroot,'toolbox','matlab','sparfun','bicgstab'),'(g1a,b,1e-6," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1)'));\n"; - mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n"; - mDynamicModelFile << sp << " if (flag1>0)\n"; - mDynamicModelFile << sp << " if(flag1==1)\n"; - mDynamicModelFile << sp << " disp(['No convergence inside BICGSTAB after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n"; - mDynamicModelFile << sp << " elseif(flag1==2)\n"; - mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n"; - mDynamicModelFile << sp << " elseif(flag1==3)\n"; - mDynamicModelFile << sp << " disp(['BICGSTAB stagnated. (Two consecutive iterates were the same.)']);\n"; + mDynamicModelFile << sp << " flag1=1;\n"; + mDynamicModelFile << sp << " while(flag1>0)\n"; + mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n"; + mDynamicModelFile << sp << " [za,flag1] = bicgstabh(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n"; + mDynamicModelFile << sp << " if (flag1>0 | reduced)\n"; + mDynamicModelFile << sp << " if(flag1==1)\n"; + mDynamicModelFile << sp << " disp(['No convergence inside BICGSTAB after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n"; + mDynamicModelFile << sp << " elseif(flag1==2)\n"; + mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n"; + mDynamicModelFile << sp << " elseif(flag1==3)\n"; + mDynamicModelFile << sp << " disp(['BICGSTAB stagnated. (Two consecutive iterates were the same.)']);\n"; + mDynamicModelFile << sp << " end;\n"; + mDynamicModelFile << sp << " luinc_tol = luinc_tol/10;\n"; + mDynamicModelFile << sp << " reduced = 0;\n"; + mDynamicModelFile << sp << " else\n"; + mDynamicModelFile << sp << " dx = za - ya;\n"; + mDynamicModelFile << sp << " ya = ya + lambda*dx;\n"; + mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n"; mDynamicModelFile << sp << " end;\n"; mDynamicModelFile << sp << " end;\n"; mDynamicModelFile << sp << " end;\n"; @@ -3509,7 +3618,171 @@ ModelTree::writeOutput(ostream &output) const output << ";"; } output << "]';\n"; + //In case of sparse model, writes the block structure of the model + if(mode==eSparseMode || mode==eSparseDLLMode) + { + int prev_Simulation_Type=-1; + bool skip_the_head; + int k=0; + int count_lead_lag_incidence = 0; + int max_lead, max_lag; + for(int j = 0;j < block_triangular.ModelBlock->Size;j++) + { + //For a block composed of a single equation determines wether we have to evaluate or to solve the equation + if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(block_triangular.ModelBlock->Block_List[j].Simulation_Type) + && (block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD + ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD + ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R + ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R )) + skip_the_head=true; + else + { + skip_the_head=false; + k++; + count_lead_lag_incidence = 0; + int Block_size=block_triangular.ModelBlock->Block_List[j].Size; + max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ; + max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead; + bool evaluate=false; + ostringstream tmp_s, tmp_s_eq; + tmp_s.str(""); + tmp_s_eq.str(""); + for(int i=0;iBlock_List[j].Size;i++) + { + tmp_s << " " << block_triangular.ModelBlock->Block_List[j].Variable[i]+1; + tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j].Equation[i]+1; + } + if (block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD + ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD + ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R + ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R + && j+Block_size<(block_triangular.ModelBlock->Size)) + { + bool OK=true; + evaluate=true; + while(j+Block_size<(block_triangular.ModelBlock->Size) && OK) + { + if(BlockTriangular::BlockSim(block_triangular.ModelBlock->Block_List[j].Simulation_Type)!=BlockTriangular::BlockSim(block_triangular.ModelBlock->Block_List[j+Block_size].Simulation_Type)) + OK=false; + else + { + if(max_lag Block_List[j+Block_size].Max_Lag ) + max_lag =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag ; + if(max_leadBlock_List[j+Block_size].Max_Lead) + max_lead=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead; + //cout << "block_triangular.ModelBlock->Block_List[" << j+Block_size << "].Size=" << block_triangular.ModelBlock->Block_List[j+Block_size].Size << "\n"; + for(int i=0;iBlock_List[j+Block_size].Size;i++) + { + tmp_s << " " << block_triangular.ModelBlock->Block_List[j+Block_size].Variable[i]+1; + tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j+Block_size].Equation[i]+1; + } + Block_size+=block_triangular.ModelBlock->Block_List[j+Block_size].Size; + } + //cout << "i=" << i << " max_lag=" << max_lag << " max_lead=" << max_lead << "\n"; + } + } + output << "M_.block_structure.block(" << k << ").num = " << j+1 << ";\n"; + //output << "M_.block_structure.block(" << k << ").size = " << block_triangular.ModelBlock->Block_List[j].Size << ";\n"; + output << "M_.block_structure.block(" << k << ").Simulation_Type = " << block_triangular.ModelBlock->Block_List[j].Simulation_Type << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead << ";\n"; + output << "M_.block_structure.block(" << k << ").endo_nbr = " << Block_size << ";\n"; + output << "M_.block_structure.block(" << k << ").equation = [" << tmp_s_eq.str() << "];\n"; + output << "M_.block_structure.block(" << k << ").variable = [" << tmp_s.str() << "];\n"; + tmp_s.str(""); + bool done_IM=false; + if(!evaluate) + { + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [];\n"; + for(int l=-max_lag;lBlock_List[j].Size;l_var++) + { + for(int l_equ=0;l_equBlock_List[j].Size;l_equ++) + if(tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[j].Variable[l_var]]) + { + count_lead_lag_incidence++; + if(tmp_s.str().length()) + tmp_s << " "; + tmp_s << count_lead_lag_incidence; + done_IM=true; + break; + } + if(!done_IM) + tmp_s << " 0"; + done_IM=false; + } + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n"; + tmp_s.str(""); + } + } + else + { + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n"; + for(int l=-max_lag;lBlock_List[ii].Size;l_var++) + { + for(int l_equ=0;l_equBlock_List[ii].Size;l_equ++) + if(tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]]) + { + //if(not_increm && l==-max_lag) + count_lead_lag_incidence++; + not_increm=false; + if(tmp_s.str().length()) + tmp_s << " "; + //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size; + tmp_s << count_lead_lag_incidence; + done_IM=true; + break; + } + if(!done_IM) + tmp_s << " 0"; + done_IM=false; + } + ii++; + } + output << tmp_s.str() << "\n"; + tmp_s.str(""); + } + output << "];\n"; + } + } + prev_Simulation_Type=block_triangular.ModelBlock->Block_List[j].Simulation_Type; + + } + for(int j=-block_triangular.Model_Max_Lag;jaddStatement(new Model_InfoStatement(options_list)); + options_list.clear(); +} + + void ParsingDriver::check() { diff --git a/preprocessor/include/BlockTriangular.hh b/preprocessor/include/BlockTriangular.hh index ceb96dda5..3d31c7485 100644 --- a/preprocessor/include/BlockTriangular.hh +++ b/preprocessor/include/BlockTriangular.hh @@ -57,14 +57,14 @@ public: void unfill_IM(int equation, int variable_endo, int lead_lag); void init_incidence_matrix(int nb_endo); void Print_IM(int n) const; - void Free_IM(List_IM* First_IM); + void Free_IM(List_IM* First_IM) const; void Print_SIM(bool* IM, int n) const; void Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m); bool Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int* prologue, int* epilogue, simple* Index_Var_IM, simple* Index_Equ_IM, bool Do_Normalization, bool mixing, bool* IM_0 , jacob_map j_m); - void Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM); + void Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0); void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n); void Allocate_Block(int size, int *count_Equ, int *count_Block, int type, Model_Block * ModelBlock); - void Free_Block(Model_Block* ModelBlock); + void Free_Block(Model_Block* ModelBlock) const; List_IM *First_IM ; List_IM *Last_IM ; simple *Index_Equ_IM; diff --git a/preprocessor/include/ComputingTasks.hh b/preprocessor/include/ComputingTasks.hh index aae31cf6d..f497a45bf 100644 --- a/preprocessor/include/ComputingTasks.hh +++ b/preprocessor/include/ComputingTasks.hh @@ -76,6 +76,16 @@ public: virtual void writeOutput(ostream &output, const string &basename) const; }; +class Model_InfoStatement : public Statement +{ +private: + const OptionsList options_list; +public: + Model_InfoStatement(const OptionsList &options_list_arg); + virtual void checkPass(ModFileStructure &mod_file_struct); + virtual void writeOutput(ostream &output, const string &basename) const; +}; + class StochSimulStatement : public Statement { private: diff --git a/preprocessor/include/ModelNormalization.hh b/preprocessor/include/ModelNormalization.hh index dc99b87c3..47460b94a 100644 --- a/preprocessor/include/ModelNormalization.hh +++ b/preprocessor/include/ModelNormalization.hh @@ -73,6 +73,7 @@ public: void Gr_to_IM_basic(int n0, int prologue, int epilogue, bool* IM, Equation_set *Equation,bool transpose); const SymbolTable &symbol_table; void Set_fp_verbose(bool ok); + void Free_Equation(int n, Equation_set* Equation); private: void IM_to_Gr(int n0, int prologue, int epilogue, bool* IM, Equation_set *Equation, Variable_set *Variable ); void Inits(Equation_set *Equation); @@ -83,7 +84,6 @@ private: int MeasureMatching(Equation_set *Equation); void OutputMatching(Equation_set* Equation); void Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* Index_Var_IM, Equation_set *Equation,bool mixing, bool* IM_s); - void Free_Equation(int n, Equation_set* Equation); void Free_Other(Variable_set* Variable); void Free_All(int n, Equation_set* Equation, Variable_set* Variable); int eq, eex; diff --git a/preprocessor/include/ParsingDriver.hh b/preprocessor/include/ParsingDriver.hh index d13dd5046..c9da8d05b 100644 --- a/preprocessor/include/ParsingDriver.hh +++ b/preprocessor/include/ParsingDriver.hh @@ -46,6 +46,7 @@ using namespace std; # undef yyFlexLexer #endif + //! The lexer class /*! Actually it was necessary to subclass the DynareFlexLexer class generated by Flex, since the prototype for DynareFlexLexer::yylex() was not convenient. @@ -53,7 +54,7 @@ using namespace std; class DynareFlex : public DynareFlexLexer { public: - DynareFlex(istream* in = 0, ostream* out = 0); + DynareFlex(std::istream* in = 0, ostream* out = 0); //! The main lexing function Dynare::parser::token_type lex(Dynare::parser::semantic_type *yylval, @@ -145,7 +146,7 @@ public: //! Starts parsing, and constructs the MOD file representation /*! The returned pointer should be deleted after use */ - ModFile *parse(istream &in, bool debug); + ModFile *parse(std::istream &in, bool debug); //! Reference to the lexer class DynareFlex *lexer; @@ -282,6 +283,8 @@ public: void simul(); //! Writes check command void check(); + //! Writes model_info command + void model_info(); //! Writes estimated params command void estimated_params(); //! Writes estimated params init command diff --git a/preprocessor/macro/MacroDriver.hh b/preprocessor/macro/MacroDriver.hh index cb8a29e25..63bded30b 100644 --- a/preprocessor/macro/MacroDriver.hh +++ b/preprocessor/macro/MacroDriver.hh @@ -53,12 +53,12 @@ private: class ScanContext { public: - istream *input; + std::istream *input; struct yy_buffer_state *buffer; const Macro::parser::location_type yylloc; const string for_body; const Macro::parser::location_type for_body_loc; - ScanContext(istream *input_arg, struct yy_buffer_state *buffer_arg, + ScanContext(std::istream *input_arg, struct yy_buffer_state *buffer_arg, Macro::parser::location_type &yylloc_arg, const string &for_body_arg, Macro::parser::location_type &for_body_loc_arg) : input(input_arg), buffer(buffer_arg), yylloc(yylloc_arg), for_body(for_body_arg), @@ -70,7 +70,7 @@ private: //! Input stream used for initialization of current scanning context /*! Kept for deletion at end of current scanning buffer */ - istream *input; + std::istream *input; //! If current context is the body of a loop, contains the string of the loop body. Empty otherwise. string for_body; @@ -125,7 +125,7 @@ private: and initialise a new scanning context with the loop body */ bool iter_loop(MacroDriver &driver, Macro::parser::location_type *yylloc); public: - MacroFlex(istream* in = 0, ostream* out = 0); + MacroFlex(std::istream* in = 0, ostream* out = 0); //! The main lexing function Macro::parser::token_type lex(Macro::parser::semantic_type *yylval,