diff --git a/matlab/model_info.m b/matlab/model_info.m index 3886735d5..985297344 100644 --- a/matlab/model_info.m +++ b/matlab/model_info.m @@ -72,17 +72,20 @@ function model_info; function ret=Sym_type(type); - EVALUATE_FORWARD=0; - EVALUATE_BACKWARD=1; - SOLVE_FORWARD_SIMPLE=2; - SOLVE_BACKWARD_SIMPLE=3; - SOLVE_TWO_BOUNDARIES_SIMPLE=4; - SOLVE_FORWARD_COMPLETE=5; - SOLVE_BACKWARD_COMPLETE=6; - SOLVE_TWO_BOUNDARIES_COMPLETE=7; - EVALUATE_FORWARD_R=8; - EVALUATE_BACKWARD_R=9; + UNKNOWN=0; + EVALUATE_FORWARD=1; + EVALUATE_BACKWARD=2; + SOLVE_FORWARD_SIMPLE=3; + SOLVE_BACKWARD_SIMPLE=4; + SOLVE_TWO_BOUNDARIES_SIMPLE=5; + SOLVE_FORWARD_COMPLETE=6; + SOLVE_BACKWARD_COMPLETE=7; + SOLVE_TWO_BOUNDARIES_COMPLETE=8; + EVALUATE_FORWARD_R=9; + EVALUATE_BACKWARD_R=10; switch (type) + case (UNKNOWN), + ret='UNKNOWN '; case {EVALUATE_FORWARD,EVALUATE_FORWARD_R}, ret='EVALUATE FORWARD '; case {EVALUATE_BACKWARD,EVALUATE_BACKWARD_R}, diff --git a/mex/sources/simulate/Interpreter.cc b/mex/sources/simulate/Interpreter.cc index 36f63c042..68292dd9c 100644 --- a/mex/sources/simulate/Interpreter.cc +++ b/mex/sources/simulate/Interpreter.cc @@ -62,7 +62,7 @@ Interpreter::pow1(double a, double b) void -Interpreter::compute_block_time() /*throw(EvalException)*/ +Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ { int var, lag, op; double v1, v2; @@ -586,7 +586,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba #endif //SparseMatrix sparse_matrix; - int nb_endo, u_count_init; + //int nb_endo, u_count_init; //mexPrintf("simulate_a_block\n"); @@ -608,12 +608,20 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba { set_code_pointer(begining); Per_y_=it_*y_size; - compute_block_time(); + compute_block_time(0); #ifdef PRINT_OUT for (j = 0; j in SOLVE_BACKWARD_SIMPLE : OK\n"); - mexEvalString("drawnow;"); + compute_block_time(0); + /*mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_SIMPLE : OK\n"); + mexEvalString("drawnow;");*/ y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0]; double rr; rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]); @@ -781,19 +789,39 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba } break;*/ case SOLVE_FORWARD_COMPLETE : -//#ifdef DEBUGC - mexPrintf("SOLVE_FORWARD_COMPLETE\n"); -//#endif +#ifdef DEBUGC + mexPrintf("SOLVE FORWARD_COMPLETE\n"); + mexPrintf("confirmation!\n"); + mexEvalString("drawnow;"); +#endif is_linear=get_code_bool; + //mexPrintf("is_linear=%d\n",is_linear); + //mexEvalString("drawnow;"); max_lag_plus_max_lead_plus_1=get_code_int; + //mexPrintf("max_lag_plus_max_lead_plus_1=%d\n",max_lag_plus_max_lead_plus_1); + //mexEvalString("drawnow;"); symbol_table_endo_nbr=get_code_int; + //mexPrintf("symbol_table_endo_nbr=%d\n",symbol_table_endo_nbr); + //mexEvalString("drawnow;"); Block_List_Max_Lag=get_code_int; + //mexPrintf("Block_List_Max_Lag=%d\n",Block_List_Max_Lag); + //mexEvalString("drawnow;"); Block_List_Max_Lead=get_code_int; + //mexPrintf("Block_List_Max_Lead=%d\n",Block_List_Max_Lead); + //mexEvalString("drawnow;"); + u_count_int=get_code_int + //mexPrintf("u_count_int=%d\n",u_count_int); + //mexEvalString("drawnow;"); + fixe_u(&u, u_count_int, u_count_int); //Read_file(file_name, periods, 0, symbol_table_endo_nbr, Block_List_Max_Lag, Block_List_Max_Lead, nb_endo, u_count, u_count_init, u); //sparse_matrix.initialize(periods, nb_endo, y_kmin, y_kmax, y_size, u_count, u_count_init, u, y, ya, slowc, y_decal, markowitz_c, res1, res2, max_res); + Read_SparseMatrix(bin_basename, size, 1, 0, 0); + //mexPrintf("size=%d\n",size); g1=(double*)mxMalloc(size*size*sizeof(double)); r=(double*)mxMalloc(size*sizeof(double)); begining=get_code_pointer; + Per_u_ = 0; + if (!is_linear) { for (it_=y_kmin;it_maxit_))) { set_code_pointer(begining); - compute_block_time(); - mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK\n"); - mexEvalString("drawnow;"); - Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + compute_block_time(0); + /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK\n"); + mexEvalString("drawnow;");*/ + //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter); res2=0; res1=0; max_res=0; @@ -841,15 +870,26 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba res1=res2=max_res=0; /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE before compute_block_time OK\n"); mexEvalString("drawnow;");*/ - compute_block_time(); + compute_block_time(0); + //mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK\n"); /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : %d OK\n",it_); mexEvalString("drawnow;");*/ //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); - simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, /*true*/false, cvg, iter); + //simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, /*true*/false, cvg, iter); + simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter); //simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter); } + /*mexPrintf("solve forward complete simulation\n"); + for (it_=y_kmin;it_maxit_))) { set_code_pointer(begining); - compute_block_time(); - mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_COMPLETE : OK\n"); - mexEvalString("drawnow;"); - Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + compute_block_time(0); + /*mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_COMPLETE : OK\n"); + mexEvalString("drawnow;");*/ + //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter); res2=0; res1=0; max_res=0; @@ -911,16 +957,19 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba { set_code_pointer(begining); Per_y_=it_*y_size; - compute_block_time(); - Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 1, false, iter); + compute_block_time(0); + //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 1, false, iter); + simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter); } } + memset(direction,0,size_of_direction); mxFree(g1); mxFree(r); mxFree(u); break; case SOLVE_TWO_BOUNDARIES_SIMPLE : case SOLVE_TWO_BOUNDARIES_COMPLETE: + mexPrintf("omp_get_max_threads=%d\n",omp_get_max_threads()); #ifdef DEBUGC mexPrintf("SOLVE_TWO_BOUNDARIES_COMPLETE\n"); mexEvalString("drawnow;"); @@ -982,7 +1031,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba Per_u_=0; Per_y_=it_*y_size; set_code_pointer(begining); - compute_block_time(); + compute_block_time(0); linbcg.Initialize(filename, res1, res2, max_res, slowc, ya, direction, iter); linbcg.Preconditioner(periods, y_kmin, y_kmax, size, IM_i, index_vara, index_equa, y_size, y, true, 0, a, indx); #endif @@ -1012,7 +1061,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba //mexPrintf("ok\n"); //mexPrintf("compute_block_time\n"); set_code_pointer(begining); - compute_block_time(); + compute_block_time(Per_u_); //mexPrintf("end of compute_block_time\n"); /*if(Gaussian_Elimination) initialize(periods, nb_endo, y_kmin, y_kmax, y_size, u_count, u_count_init, u, y, ya, slowc, y_decal, markowitz_c, res1, res2, max_res);*/ @@ -1084,7 +1133,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba Per_u_=(it_-y_kmin)*max_lag_plus_max_lead_plus_1; Per_y_=it_*y_size; set_code_pointer(begining); - compute_block_time(); + compute_block_time(Per_u_); #ifdef PRINT_OUT for (j=0; j Block_Contain; diff --git a/mex/sources/simulate/Mem_Mngr.cc b/mex/sources/simulate/Mem_Mngr.cc index a3c2a2e3d..f48ffd889 100644 --- a/mex/sources/simulate/Mem_Mngr.cc +++ b/mex/sources/simulate/Mem_Mngr.cc @@ -277,7 +277,7 @@ void Mem_Mngr::Free_All() { int i; - for (int i=0;i(&var), sizeof(var)); SaveCode.read(reinterpret_cast(&lag), sizeof(lag)); SaveCode.read(reinterpret_cast(&j), sizeof(j)); - /*mexPrintf("0i=%d eq=%d var=%d lag=%d j=%d\n",i,eq,var,lag,j ); - mexEvalString("drawnow;");*/ + //mexPrintf("i=%d eq=%d var=%d lag=%d j=%d\n",i,eq,var,lag,j ); + //mexEvalString("drawnow;"); IM_i[std::make_pair(std::make_pair(eq, var), lag)] = j; /*mexPrintf("1i=%d\n",i); mexEvalString("drawnow;");*/ @@ -419,8 +419,6 @@ void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, in #ifdef MEM_ALLOC_CHK mexPrintf("index_vara=(int*)mxMalloc(%d*sizeof(int))\n",Size*(periods+y_kmin+y_kmax)); #endif - /*mexPrintf("Size=%d periods=%d y_kmin=%d y_kmax=%d\n",Size,periods, y_kmin, y_kmax); - mexEvalString("drawnow;");*/ index_vara=(int*)mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int)); #ifdef MEM_ALLOC_CHK mexPrintf("ok\n"); @@ -430,18 +428,21 @@ void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, in SaveCode.read(reinterpret_cast(&index_vara[j]), sizeof(*index_vara)); //mexPrintf("index_vara[%d]=%d\n",j,index_vara[j]); } - for (i=1;i1) { - for (j=0;j ,int>, int> IM) +{ + int i, eq, var, lag; + //double tmp_b=0.0; + std::map ,int>, int>::iterator it4; + NonZeroElem* first; + //mexPrintf("periods=%d, y_kmin=%d, y_kmax=%d, SizeInit=%d, IM.size()=%d\n",periods, y_kmin, y_kmax, Size, IM.size()); + pivot=(int*)mxMalloc(Size*sizeof(int)); + pivotk=(int*)mxMalloc(Size*sizeof(int)); + pivotv=(double*)mxMalloc(Size*sizeof(double)); + pivotva=(double*)mxMalloc(Size*sizeof(double)); + b=(int*)mxMalloc(Size*sizeof(int)); + line_done=(bool*)mxMalloc(Size*sizeof(bool)); + memset(line_done, 0, Size*sizeof(*line_done)); + mem_mngr.init_CHUNK_BLCK_SIZE(u_count); + g_save_op=NULL; + g_nop_all=0; + i=Size*sizeof(NonZeroElem*); + FNZE_R=(NonZeroElem**)mxMalloc(i); + FNZE_C=(NonZeroElem**)mxMalloc(i); + memset(FNZE_R, 0, i); + memset(FNZE_C, 0, i); + NonZeroElem** temp_NZE_R=(NonZeroElem**)mxMalloc(i); + NonZeroElem** temp_NZE_C=(NonZeroElem**)mxMalloc(i); + memset(temp_NZE_R, 0, i); + memset(temp_NZE_C, 0, i); + i=Size*sizeof(int); + NbNZRow=(int*)mxMalloc(i); + NbNZCol=(int*)mxMalloc(i); + memset(NbNZRow, 0, i); + memset(NbNZCol, 0, i); + i=Size*sizeof(*b); + memset(b,0,i); + it4=IM.begin(); + eq=-1; + double tmp_b[Size]; + for(i=0; i< Size;i++); + tmp_b[i]=0;//u[i]; + int u_count1=Size; + while (it4!=IM.end()) + { + var=it4->first.first.second; + /*if (eq!=it4->first.first.first) + tmp_b=0;*/ + eq=it4->first.first.first; + lag=it4->first.second; + if (lag==0) /*Build the index for sparse matrix containing the jacobian : u*/ + { + //mexPrintf("Add eq=%d, var=%d, lag=%d at it_=%d u=%f\n",eq,var,lag, it_, u[u_count1]); + //mexPrintf(" u_index=%d\n",/*it4->second+u_count_init*it_*/u_count1); + NbNZRow[eq]++; + NbNZCol[var]++; +#ifdef NEW_ALLOC + first=mem_mngr.mxMalloc_NZE(); +#else + first=(NonZeroElem*)mxMalloc(sizeof(NonZeroElem)); +#endif + first->NZE_C_N=NULL; + first->NZE_R_N=NULL; + first->u_index=u_count1/*it4->second+u_count_init*it_*/; + first->r_index=eq; + first->c_index=var; + first->lag_index=lag; + //mexPrintf(" u[%d](%f)*y[%d](%f)=%f\n",u_count1, u[u_count1], index_vara[var]+it_*y_size, y[index_vara[var]+it_*y_size], u[u_count1]*y[index_vara[var]+it_*y_size]); + tmp_b[eq] += u[u_count1]*y[index_vara[var]+it_*y_size]; + if (FNZE_R[eq]==NULL) + { + FNZE_R[eq]=first; + } + if (FNZE_C[var]==NULL) + FNZE_C[var]=first; + if (temp_NZE_R[eq]!=NULL) + temp_NZE_R[eq]->NZE_R_N=first; + if (temp_NZE_C[var]!=NULL) + temp_NZE_C[var]->NZE_C_N=first; + temp_NZE_R[eq]=first; + temp_NZE_C[var]=first; + u_count1++; + } + it4++; + } + for(i=0;i ,int>, int> IM) { int t,i, eq, var, lag; double tmp_b=0.0; std::map ,int>, int>::iterator it4; NonZeroElem* first; - //mexPrintf("periods=%d, y_kmin=%d, y_kmax=%d, Size=%d, IM.size()=%d\n",periods, y_kmin, y_kmax, Size, IM.size()); + //mexPrintf("periods=%d, y_kmin=%d, y_kmax=%d, SizeInit=%d, IM.size()=%d\n",periods, y_kmin, y_kmax, Size, IM.size()); #ifdef MEM_ALLOC_CHK mexPrintf("pivot=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods); #endif @@ -1319,73 +1416,6 @@ SparseMatrix::complete(int beg_t, int Size, int periods, int *b) return(beg_t); } -double -SparseMatrix::bksub( int tbreak, int last_period, int Size, double slowc_l -#ifdef PROFILER -, /*NonZeroElem *first,*/ long int *nmul -#endif -) -{ - NonZeroElem *first; - int i, j, k; - double yy; - res1 = res2 = max_res = 0; - for (i=0;i=y_kmin;t--) - { - int ti=(t-y_kmin)*Size; -#ifdef PRINT_OUT - mexPrintf("t=%d ti=%d\n",t,ti); -#endif - int cal=y_kmin*Size; - int cal_y=y_size*y_kmin; - for (i=ti-1;i>=ti-Size;i--) - { - j=i+cal; -#ifdef PRINT_OUT_y - mexPrintf("t=%d, ti=%d j=%d i+Size=%d\n",t,ti,j,i+Size); -#endif - int pos=pivot[/*j*/i+Size]; -#ifdef PRINT_OUT_y - mexPrintf("i-ti+Size=%d pos=%d j=%d\n",i-ti+Size,pos,j); -#endif - int nb_var=At_Row(pos,&first); - first=first->NZE_R_N; - nb_var--; - int eq=index_vara[j]+y_size; -#ifdef PRINT_OUT_y1 - mexPrintf("y[index_vara[%d]=%d]=",j,index_vara[j]+y_size); -#endif - yy=0; - for (k=0;kc_index]+cal_y]*u[first->u_index]; -#ifdef PROFILER - (*nmul)++; -#endif - first=first->NZE_R_N; - } -#ifdef PRINT_OUT_y1 - mexPrintf("|u[%d](%f)|",b[pos],double(u[b[pos]])); -#endif - yy=-(yy+y[eq]+u[b[pos]]); - direction[eq]=yy; - //mexPrintf("direction[%d] = %f\n",eq,yy); - y[eq] += slowc_l*yy; -#ifdef PRINT_OUT_y1 - mexPrintf("=%f (%f)\n",double(yy),double(y[eq])); -#endif - } - } - return res1; -} - - void SparseMatrix::close_swp_file() { @@ -1955,6 +1985,384 @@ SparseMatrix::Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_k +double +SparseMatrix::bksub( int tbreak, int last_period, int Size, double slowc_l +#ifdef PROFILER +, /*NonZeroElem *first,*/ long int *nmul +#endif +) +{ + NonZeroElem *first; + int i, j, k; + double yy; + res1 = res2 = max_res = 0; + for (i=0;i=y_kmin;t--) + { + int ti=(t-y_kmin)*Size; +#ifdef PRINT_OUT + mexPrintf("t=%d ti=%d\n",t,ti); +#endif + int cal=y_kmin*Size; + int cal_y=y_size*y_kmin; + for (i=ti-1;i>=ti-Size;i--) + { + j=i+cal; +#ifdef PRINT_OUT_y + mexPrintf("t=%d, ti=%d j=%d i+Size=%d\n",t,ti,j,i+Size); +#endif + int pos=pivot[/*j*/i+Size]; +#ifdef PRINT_OUT_y + mexPrintf("i-ti+Size=%d pos=%d j=%d\n",i-ti+Size,pos,j); +#endif + int nb_var=At_Row(pos,&first); + first=first->NZE_R_N; + nb_var--; + int eq=index_vara[j]+y_size; +#ifdef PRINT_OUT_y1 + mexPrintf("y[index_vara[%d]=%d]=",j,index_vara[j]+y_size); +#endif + yy=0; + for (k=0;kc_index]+cal_y]*u[first->u_index]; +#ifdef PROFILER + (*nmul)++; +#endif + first=first->NZE_R_N; + } +#ifdef PRINT_OUT_y1 + mexPrintf("|u[%d](%f)|",b[pos],double(u[b[pos]])); +#endif + yy=-(yy+y[eq]+u[b[pos]]); + direction[eq]=yy; + //mexPrintf("direction[%d] = %f\n",eq,yy); + y[eq] += slowc_l*yy; +#ifdef PRINT_OUT_y1 + mexPrintf("=%f (%f)\n",double(yy),double(y[eq])); +#endif + } + } + return res1; +} + + +double +SparseMatrix::simple_bksub(int it_, int Size, double slowc_l) +{ + int i,k; + double yy; + NonZeroElem *first; + res1 = res2 = max_res = 0; + //mexPrintf("simple_bksub\n"); + for (i=0;i=0;i--) + { + //mexPrintf("i=%d\n",i); + int pos=pivot[i]; + //mexPrintf("pos=%d\n",pos); + int nb_var=At_Row(pos,&first); + //mexPrintf("nb_var=%d\n",nb_var); + first=first->NZE_R_N; + nb_var--; + //mexPrintf("i=%d\n",i); + int eq=index_vara[i]; + //mexPrintf("eq=%d\n",eq); + yy = 0; + for (k=0;kc_index, index_vara[first->c_index], y[index_vara[first->c_index]+it_*y_size]); + //mexPrintf("u[first->u_index=%d]=%f\n",first->u_index, u[first->u_index]); + yy+=y[index_vara[first->c_index]+it_*y_size]*u[first->u_index]; + first=first->NZE_R_N; + } + yy=-(yy+y[eq+it_*y_size]+u[b[pos]]); + //mexPrintf("yy=%f\n",yy); + direction[eq+it_*y_size]=yy; + y[eq+it_*y_size] += slowc_l*yy; + } + return res1; +} + + +int +SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter) +{ + int i, j, k; + int pivj=0, pivk=0; + double piv_abs, first_elem; + NonZeroElem *first, *firsta, *first_sub, *first_piv, *first_suba; +#ifdef MARKOVITZ + double piv_v[Size]; + int pivj_v[Size], pivk_v[Size], NR[Size], l, N_max; + bool one; +#endif + + /*mexPrintf("In simulate_NG it_=%d\n",it_); + mexEvalString("drawnow;");*/ + Simple_Init(it_, y_kmin, y_kmax, Size, IM_i); + if (isnan(res1) || isinf(res1)) + { + if (slowc_save<1e-8) + { + mexPrintf("Dynare cannot improve the simulation\n"); + mexEvalString("drawnow;"); + mexEvalString("st=fclose('all');clear all;"); + filename+=" stopped"; + mexErrMsgTxt(filename.c_str()); + } + slowc_save/=2; + mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n",slowc_save); + for (i=0;ir_index=%d\n", first->r_index); + if (!line_done[first->r_index]) + { + k=first->u_index; + //mexPrintf("k=%d\n",k); + int jj=first->r_index; + //mexPrintf("jj=%d\n",jj); + int NRow_jj=NRow(jj); + //mexPrintf("NRow_jj=%d\n",NRow_jj); +#ifdef PROFILER + nbpivot++; + nbpivot_it++; +#endif + +#ifdef MARKOVITZ + //mexPrintf("l=%d\n",l); + piv_v[l]=u[k]; + double piv_fabs=fabs(u[k]); + pivj_v[l]=jj; + pivk_v[l]=k; + NR[l]=NRow_jj; + //mexPrintf("one=%d\n",one); + if (NRow_jj==1 && !one) + { + one=true; + piv_abs=piv_fabs; + N_max=NRow_jj; + } + if (!one) + { + if (piv_fabs>piv_abs) + piv_abs=piv_fabs; + if (NRow_jj>N_max) + N_max=NRow_jj; + } + else + { + if (NRow_jj==1) + { + if (piv_fabs>piv_abs) + piv_abs=piv_fabs; + if (NRow_jj>N_max) + N_max=NRow_jj; + } + } + l++; +#else + if (piv_absNZE_C_N=%x\n",j,first->NZE_C_N); + first=first->NZE_C_N; + } + //mexPrintf("Before Markovitz\n"); +#ifdef MARKOVITZ + double markovitz=0, markovitz_max=-9e70; + //mexPrintf("l=%d one=%d\n",l, one); + if (!one) + { + for (j=0;jmarkovitz_max) + { + piv=piv_v[j]; + pivj=pivj_v[j]; //Line number + pivk=pivk_v[j]; //positi + markovitz_max=markovitz; + } + } + } + else + { + for (j=0;jmarkovitz_max && NR[j]==1) + { + piv=piv_v[j]; + pivj=pivj_v[j]; //Line number + pivk=pivk_v[j]; //positi + markovitz_max=markovitz; + } + } + } +#endif + //mexPrintf("OK i=%d\n", i); + pivot[i]=pivj; + //mexPrintf("pivot[%d]=%d\n",i,pivot[i]); + pivotk[i]=pivk; + pivotv[i]=piv; + line_done[pivj]=true; + if (piv_absu_index]/=piv; + first=first->NZE_R_N; + } + u[b[pivj]]/=piv; + /*substract the elements on the non treated lines*/ + nb_eq=At_Col(i,&first); + NonZeroElem *first_piva; + int nb_var_piva=At_Row(pivj,&first_piva); + for (j=0;jr_index; + if (!line_done[row]) + { + first_elem=u[first->u_index]; + int nb_var_piv=nb_var_piva; + first_piv=first_piva; + int nb_var_sub=At_Row(row,&first_sub); + int l_sub=0, l_piv=0; + int sub_c_index=first_sub->c_index, piv_c_index=first_piv->c_index; + //int tmp_lag=first_sub->lag_index; + while (l_sub=nb_var_piv)) + { + first_sub=first_sub->NZE_R_N; + if (first_sub) + sub_c_index=first_sub->c_index; + else + sub_c_index=Size/* *periods*/; + l_sub++; + } + else if (sub_c_index>piv_c_index || l_sub>=nb_var_sub) + { + int tmp_u_count=Get_u(); + //int lag=first_piv->c_index/Size-row/Size; + Insert(row,first_piv->c_index,tmp_u_count,/*lag*/0); + u[tmp_u_count]=-u[first_piv->u_index]*first_elem; + first_piv=first_piv->NZE_R_N; + if (first_piv) + piv_c_index=first_piv->c_index; + else + piv_c_index=Size/* *periods*/; + l_piv++; + } + else /*first_sub->c_index==first_piv->c_index*/ + { + if (i==sub_c_index) + { + firsta=first; + first_suba=first_sub->NZE_R_N; + Delete(first_sub->r_index,first_sub->c_index, Size, b); + first=firsta->NZE_C_N; + first_sub=first_suba; + if (first_sub) + sub_c_index=first_sub->c_index; + else + sub_c_index=Size/* *periods*/; + l_sub++; + first_piv=first_piv->NZE_R_N; + if (first_piv) + piv_c_index=first_piv->c_index; + else + piv_c_index=Size; + l_piv++; + } + else + { + u[first_sub->u_index]-=u[first_piv->u_index]*first_elem; + first_sub=first_sub->NZE_R_N; + if (first_sub) + sub_c_index=first_sub->c_index; + else + sub_c_index=Size /* *periods*/; + l_sub++; + first_piv=first_piv->NZE_R_N; + if (first_piv) + piv_c_index=first_piv->c_index; + else + piv_c_index=Size/* *periods*/; + l_piv++; + } + } + } + u[b[row]]-=u[b[pivj]]*first_elem; + } + else + first=first->NZE_C_N; + } + } + double slowc_lbx=slowc, res1bx; + //mexPrintf("before bksub it_=%d\n",it_); + for (i=0;iu_count=%d &u_count=%x\n",u_count,&u_count); + mexPrintf("GNU version=%d\n",GNUVER); #ifdef RECORD_ALL int save_u_count=u_count; #endif @@ -2208,8 +2617,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax bool one=false; piv_abs=0; #endif + //#pragma omp parallel for for (j=0;jr_index,line_done[first->r_index]); @@ -2275,6 +2686,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax double markovitz=0, markovitz_max=-9e70; if (!one) { + /*#pragma omp parallel for private(j) shared(markovitz_max, piv, pivj, pivk) for (j=0;jmarkovitz_max) + { + piv=piv_v[j]; + pivj=pivj_v[j]; //Line number + pivk=pivk_v[j]; //positi + markovitz_max=markovitz; + //printf("j=%d, markovitz_max=%f\n",j, markovitz); + } + } } } else { - for (j=0;jmarkovitz_max && NR[j]==1) + //#pragma omp for schedule(dynamic, 10) nowait + //#pragma omp parallel for private(j, markovitz) + #pragma omp parallel for private(j,markovitz) + for (j=0;jmarkovitz_max && NR[j]==1) + { + piv=piv_v[j]; + pivj=pivj_v[j]; //Line number + pivk=pivk_v[j]; //positi + markovitz_max=markovitz; + } } } } @@ -2373,7 +2811,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax #endif if (piv_absr_index; @@ -3099,8 +3539,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax void SparseMatrix::fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1) { + //mexPrintf("u_count_int=%d\n",u_count_int); u_count=u_count_int * periods; u_count_alloc = 2*u_count; + //mexPrintf("mxMalloc(%d*sizeof(double)=%d)\n",u_count_alloc,u_count_alloc*sizeof(double)); (*u)=(double*)mxMalloc(u_count_alloc*sizeof(double)); memset((*u), 0, u_count_alloc*sizeof(double)); u_count_init=max_lag_plus_max_lead_plus_1; diff --git a/mex/sources/simulate/SparseMatrix.hh b/mex/sources/simulate/SparseMatrix.hh index 1463d27b5..d8dc90a44 100644 --- a/mex/sources/simulate/SparseMatrix.hh +++ b/mex/sources/simulate/SparseMatrix.hh @@ -26,7 +26,12 @@ #include #include #include "Mem_Mngr.hh" - +//! Openmp is available in GCC since version 4.3.2 +//! Test if GCC version is greater then 4.3.2 order to avoid a copilation error +#define GNUVER 100*__GNUC__+10*__GNUC_MINOR__+__GNUC_PATCHLEVEL__ +#if GNUVER >= 432 + #include +#endif #define NEW_ALLOC #define MARKOVITZ //#define MEMORY_LEAKS @@ -71,6 +76,7 @@ class SparseMatrix public: SparseMatrix(); int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter); + int simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter); void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter); void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1); //void initialize(int periods_arg, int nb_endo_arg, int y_kmin_arg, int y_kmax_arg, int y_size_arg, int u_count_arg, int u_count_init_arg, double *u_arg, double *y_arg, double *ya_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg, double res1_arg, double res2_arg, double max_res_arg); @@ -80,6 +86,7 @@ class SparseMatrix private: void Init(int periods, int y_kmin, int y_kmax, int Size, map ,int>, int> IM); void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map ,int>, int> IM); + void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map ,int>, int> IM); void End(int Size); bool compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size #ifdef PROFILER @@ -106,6 +113,7 @@ class SparseMatrix , long int *nmul #endif ); + double simple_bksub(int it_, int Size, double slowc_l); void run_triangular(int nop_all,int *op_all); void run_it(int nop_all,int *op_all); void run_u_period1(int periods); diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc index 79b1b9ac5..879e6bd37 100644 --- a/preprocessor/ExprNode.cc +++ b/preprocessor/ExprNode.cc @@ -147,10 +147,11 @@ NumConstNode::eval(const eval_context_type &eval_context) const throw (EvalExcep } void -NumConstNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const +NumConstNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const { CompileCode.write(&FLDC, sizeof(FLDC)); double vard=atof(datatree.num_constants.get(id).c_str()); + //double vard=id; #ifdef DEBUGC cout << "FLDC " << vard << "\n"; #endif @@ -427,7 +428,7 @@ VariableNode::eval(const eval_context_type &eval_context) const throw (EvalExcep } void -VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const +VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const { int i, lagl; #ifdef DEBUGC @@ -892,7 +893,7 @@ UnaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExcept } void -UnaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const +UnaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const { temporary_terms_type::const_iterator it = temporary_terms.find(const_cast(this)); if (it != temporary_terms.end()) @@ -1164,8 +1165,6 @@ BinaryOpNode::computeTemporaryTerms(map &reference_count, reference_count[this2]++; if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C) { - if(this2->idx==2280) - cout << "==>Curr_block= " << Curr_block << " equation= " << equation << " first_occurence[this2].first=" << first_occurence[this2].first << " first_occurence[this2].second=" << first_occurence[this2].second << "\n"; temporary_terms.insert(this2); ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2); } @@ -1226,7 +1225,7 @@ BinaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExcep } void -BinaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const +BinaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const { // If current node is a temporary term temporary_terms_type::const_iterator it = temporary_terms.find(const_cast(this)); @@ -1580,8 +1579,6 @@ TrinaryOpNode::computeTemporaryTerms(map &reference_count, reference_count[this2]++; if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C) { - if(this2->idx==2280) - cout << "==>Curr_block= " << Curr_block << " equation= " << equation << " first_occurence[this2].first=" << first_occurence[this2].first << " first_occurence[this2].second=" << first_occurence[this2].second << "\n"; temporary_terms.insert(this2); ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2); } @@ -1613,7 +1610,7 @@ TrinaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExce void TrinaryOpNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, - const temporary_terms_type &temporary_terms, map_idx_type map_idx) const + const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const { // If current node is a temporary term temporary_terms_type::const_iterator it = temporary_terms.find(const_cast(this)); @@ -1788,7 +1785,7 @@ UnknownFunctionNode::eval(const eval_context_type &eval_context) const throw (Ev } void -UnknownFunctionNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const +UnknownFunctionNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const { cerr << "UnknownFunctionNode::compile: operation impossible!" << endl; exit(EXIT_FAILURE); diff --git a/preprocessor/Makefile b/preprocessor/Makefile index 4dd2f09c1..40d47a87e 100644 --- a/preprocessor/Makefile +++ b/preprocessor/Makefile @@ -33,7 +33,6 @@ OBJS = \ IncidenceMatrix.o \ BlockTriangular.o \ Model_Graph.o \ - SymbolGaussElim.o \ DynareMain.o \ DynareMain2.o diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index c2081049d..68a818f43 100644 --- a/preprocessor/ModFile.cc +++ b/preprocessor/ModFile.cc @@ -125,7 +125,7 @@ ModFile::evalAllExpressions() cout << "error in evaluation of variable\n"; } } - if(mod_file_struct.load_params_and_steady_state_present && init_values.size() 0) { - if (model_tree.mode == eDLLMode) + if (model_tree.mode == eDLLMode || model_tree.mode == eSparseDLLMode) { mOutputFile << "if exist('" << basename << "_static.c')" << endl; mOutputFile << " clear " << basename << "_static" << endl; @@ -425,7 +425,6 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const mOutputFile << "erase_compiled_function('" + basename +"_dynamic');" << endl; } } - cout << "Processing outputs ..."; symbol_table.writeOutput(mOutputFile); @@ -448,9 +447,10 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const mOutputFile << "save('" << basename << "_results.mat', 'oo_', 'M_', 'options_');" << endl; mOutputFile << "diary off" << endl; - if (model_tree.mode == eSparseMode) + if (model_tree.mode == eSparseMode || model_tree.mode == eSparseDLLMode) mOutputFile << "rmpath " << basename << ";\n"; mOutputFile << endl << "disp(['Total computing time : ' dynsec2hms(toc) ]);" << endl; mOutputFile.close(); + cout << "done\n"; } diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index 9c972fba2..4d369c29d 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -73,7 +73,7 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag, } void -ModelTree::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type map_idx) const +ModelTree::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const { first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(eEndogenous, symb_id, lag))); if (it != first_derivatives.end()) @@ -1082,7 +1082,7 @@ end: void -ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type) const +ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const { struct Uff_l { @@ -1107,7 +1107,8 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl map reference_count; map ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number; int prev_Simulation_Type=-1; - SymbolicGaussElimination SGE; + //SymbolicGaussElimination SGE; + bool file_open=false; temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); //---------------------------------------------------------------------- string main_name=file_name; @@ -1180,14 +1181,15 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl code_file.write(reinterpret_cast(&v),sizeof(v)); v=block_triangular.ModelBlock->Block_List[j].Max_Lead; code_file.write(reinterpret_cast(&v),sizeof(v)); - if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) - { + //if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) + //{ int u_count_int=0; - Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,SGE.file_open); + Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open, + ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE); v=u_count_int; code_file.write(reinterpret_cast(&v),sizeof(v)); - SGE.file_is_open(); - } + file_open=true; + //} } for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++) { @@ -1821,7 +1823,7 @@ ModelTree::reform(const string name1) const void ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, const int &num, - int &u_count_int, bool &file_open) const + int &u_count_int, bool &file_open, bool is_two_boundaries) const { int j; std::ofstream SaveCode; @@ -1850,19 +1852,22 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b u_count_int++; } } - for (j=0;jBlock_List[num].Size;j++) + if(is_two_boundaries) { - int eqr1=j; - int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods - +block_triangular.incidencematrix.Model_Max_Lead_Endo); - int k1=0; - SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); - SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); - SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); - SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); - u_count_int++; + for (j=0;jBlock_List[num].Size;j++) + { + int eqr1=j; + int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods + +block_triangular.incidencematrix.Model_Max_Lead_Endo); + int k1=0; + SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); + SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); + SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); + SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); + u_count_int++; + } } - + //cout << "u_count_int=" << u_count_int << "\n"; for (j=0;jBlock_List[num].Size;j++) { int varr=block_triangular.ModelBlock->Block_List[num].Variable[j]; @@ -2103,7 +2108,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string ofstream mDynamicModelFile; ostringstream tmp, tmp1, tmp_eq; int prev_Simulation_Type, tmp_i; - SymbolicGaussElimination SGE; + //SymbolicGaussElimination SGE; bool OK; chdir(basename.c_str()); string filename = dynamic_basename + ".m"; @@ -2121,7 +2126,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << "%/\n"; int i, k, Nb_SGE=0; - bool printed = false, skip_head, open_par=false; + bool skip_head, open_par=false; if (computeJacobian || computeJacobianExo || computeHessian) { mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n"; @@ -2413,10 +2418,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string if (open_par) mDynamicModelFile << " end\n"; open_par=false; - if (!printed) - { - printed = true; - } Nb_SGE++; int nze, m; for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++) @@ -2453,8 +2454,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string writeModelEquationsOrdered_M( block_triangular.ModelBlock, dynamic_basename); chdir(".."); - if (printed) - cout << "done\n"; } void @@ -3074,9 +3073,10 @@ ModelTree::writeStaticFile(const string &basename) const switch (mode) { case eStandardMode: - case eSparseDLLMode: + /*case eSparseDLLMode:*/ writeStaticMFile(basename + "_static"); break; + case eSparseDLLMode: case eSparseMode: // create a directory to store all files #ifdef _WIN32 @@ -3117,7 +3117,7 @@ ModelTree::writeDynamicFile(const string &basename) const #else mkdir(basename.c_str(), 0777); #endif - writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL); + writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL, map_idx); block_triangular.Free_Block(block_triangular.ModelBlock); block_triangular.incidencematrix.Free_IM(); //block_triangular.Free_IM_X(block_triangular.First_IM_X); diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc index 1c106e7be..3ff818bc8 100644 --- a/preprocessor/ParsingDriver.cc +++ b/preprocessor/ParsingDriver.cc @@ -596,7 +596,7 @@ ParsingDriver::add_to_row(NodeID v) void ParsingDriver::steady() { - if (mod_file->model_tree.mode == eSparseMode) + if (mod_file->model_tree.mode == eSparseMode || mod_file->model_tree.mode == eSparseDLLMode) mod_file->addStatement(new SteadySparseStatement(options_list)); else mod_file->addStatement(new SteadyStatement(options_list)); diff --git a/preprocessor/SymbolGaussElim.cc b/preprocessor/SymbolGaussElim.cc deleted file mode 100644 index f4d3452d9..000000000 --- a/preprocessor/SymbolGaussElim.cc +++ /dev/null @@ -1,1827 +0,0 @@ -/* - * Copyright (C) 2007-2008 Dynare Team - * - * This file is part of Dynare. - * - * Dynare is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * Dynare is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with Dynare. If not, see . - */ - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "ModelTree.hh" -//#include "pctimer_h.hh" -#include "Model_Graph.hh" -#include "SymbolGaussElim.hh" -//#define SIMPLIFY - -using namespace std; -int max_nb_table_y, max_nb_in_degree_edges=0, max_nb_out_degree_edges=0; - -SymbolicGaussElimination::SymbolicGaussElimination() -{ -#ifdef SAVE - file_open = false; -#endif -} - -#ifdef SIMPLIFY - -void -SymbolicGaussElimination::set_free_u_list(int index) -{ -#ifndef UNSORTED - int i = 0, j; - if(nb_free_u_list > 0) - { - while((free_u_list[i] > index) && (i < nb_free_u_list)) - i++; - for(j = nb_free_u_list;j > i;j--) - { - free_u_list[j] = free_u_list[j - 1]; - if(j >= MAX_FREE_U_LIST) - { - cout << "Error to much j=" << j << " Maximum allowed=" << MAX_FREE_U_LIST << "\n"; - system("pause"); - exit(EXIT_FAILURE); - } - } - } - if(i >= MAX_FREE_U_LIST) - { - cout << "Error to much i=" << i << " Maximum allowed=" << MAX_FREE_U_LIST << "\n"; - system("pause"); - exit(EXIT_FAILURE); - } - free_u_list[i] = index; - nb_free_u_list++; - if(nb_free_u_list >= MAX_FREE_U_LIST) - { - cout << "Error to much free_u_list=" << nb_free_u_list << " Maximum allowed=" << MAX_FREE_U_LIST << "\n"; - system("pause"); - exit(EXIT_FAILURE); - } - if(nb_free_u_list > max_nb_free_u_list) - max_nb_free_u_list = nb_free_u_list; -#else - free_u_list[nb_free_u_list++] = index; -#endif -} - -int -SymbolicGaussElimination::get_free_u_list(bool dynamic) -{ - int i; - if((nb_free_u_list > 0) && (simplification_allowed)) - { - i = free_u_list[--nb_free_u_list]; - return (i); - } - else - if(dynamic) - return (u_count++); - else - return (u_count++); -} -#endif /**SIMPLIFY**/ - -//================================================================================== -void -SymbolicGaussElimination::init_glb() -{ - tol = TOL; - vertex_count = 0; - nb_table_u = 0; - nb_prologue_save_table_y = 0; - nb_first_save_table_y = 0; - nb_middle_save_table_y = 0; - nb_last_save_table_y = 0; - nb_prologue_save_table_u = 0; - nb_first_save_table_u = 0; - nb_middle_save_table_u = 0; - nb_last_save_table_u = 0; - middle_count_loop = 0; - u_count = 0; - y_kmin = y_kmax = 0; - prologue_save_table_u = NULL; - first_save_table_u = NULL; - first_save_i_table_u = NULL; - middle_save_table_u = NULL; - middle_save_i_table_u = NULL; - last_save_table_u = NULL; - prologue_save_table_y = NULL; - first_save_table_y = NULL; - first_save_i_table_y = NULL; - middle_save_table_y = NULL; - middle_save_i_table_y = NULL; - last_save_table_y = NULL; - pos_nb_first_save_table_u = 0; -#ifdef SIMPLIFY - nb_free_u_list = 0; - max_nb_free_u_list1 = 0; - max_nb_free_u_list = 0; -#endif - save_direct=true; -} - - -void -SymbolicGaussElimination::write_to_file_table_y( t_table_y *save_table_y, t_table_y *save_i_table_y, int nb_save_table_y, int nb_save_i_table_y) -{ - int i, k; -#ifdef PRINT_OUT - cout << "nb_save_table_y=" << nb_save_table_y << "\n"; -#endif - SaveCode.write(reinterpret_cast (&nb_save_table_y), sizeof(nb_save_table_y)); - for(i = 0;i < nb_save_table_y;i++) - { - SaveCode.write(reinterpret_cast(&save_table_y[i].index), sizeof(save_table_y[i].index)); -#ifdef PRINT_OUT - cout << "+y[" << save_table_y[i].index << "]="; -#endif - SaveCode.write(reinterpret_cast(&save_table_y[i].nb), sizeof(save_table_y[i].nb)); - for(k = 0;k < save_table_y[i].nb;k++) - { - SaveCode.write(reinterpret_cast(&save_table_y[i].u_index[k]), sizeof(save_table_y[i].u_index[k])); - SaveCode.write(reinterpret_cast(&save_table_y[i].y_index[k]), sizeof(save_table_y[i].y_index[k])); -#ifdef PRINT_OUT - cout << "u[" << save_table_y[i].u_index[k] << "]*y[" << save_table_y[i].y_index[k] << "]"; - if(k + 1 < save_table_y[i].nb) - cout << "+"; - else - cout << "\n"; -#endif - } - } - for(i = 0;i < nb_save_i_table_y;i++) - { - SaveCode.write(reinterpret_cast(&save_i_table_y[i].index), sizeof(save_i_table_y[i].index)); -#ifdef PRINT_OUT - cout << "+i y[" << save_i_table_y[i].index << "]="; -#endif - SaveCode.write(reinterpret_cast(&save_i_table_y[i].nb), sizeof(save_i_table_y[i].nb)); - for(k = 0;k < save_i_table_y[i].nb;k++) - { - SaveCode.write(reinterpret_cast(&save_i_table_y[i].u_index[k]), sizeof(save_i_table_y[i].u_index[k])); - SaveCode.write(reinterpret_cast(&save_i_table_y[i].y_index[k]), sizeof(save_i_table_y[i].y_index[k])); -#ifdef PRINT_OUT - cout << "u[" << save_i_table_y[i].u_index[k] << "]*y[" << save_i_table_y[i].y_index[k] << "]"; - if(k + 1 < save_i_table_y[i].nb) - cout << "+"; - else - cout << "\n"; -#endif - } - } -} - -void -SymbolicGaussElimination::write_to_file_table_u_b(t_table_u *save_table_u, t_table_u *last_table_u, int *nb_save_table_u, bool chk) -{ - t_table_u *table_u; - bool OK=true; - while(/*save_table_u!=last_table_u*/OK) - { -#ifdef PRINT_OUT - cout << "**save_table_u->type=" << int(save_table_u->type) << "\n"; -#endif - switch (save_table_u->type) - { - case 3: - case 7: - (*nb_save_table_u)++; - //SaveCode << save_table_u->type << save_table_u->index << save_table_u->op1; - SaveCode.write(reinterpret_cast(&save_table_u->type), sizeof(save_table_u->type)); - SaveCode.write(reinterpret_cast(&save_table_u->index), sizeof(save_table_u->index)); - SaveCode.write(reinterpret_cast(&save_table_u->op1), sizeof(save_table_u->op1)); -#ifdef PRINT_OUT - if(save_table_u->type == 3) - cout << "+u[" << save_table_u->index << "]=1/(1-u[" << save_table_u->op1 << "])\n"; - else - cout << "+u[" << save_table_u->index << "]*=u[" << save_table_u->op1 << "]\n"; -#endif /**PRINT_OUT**/ - break; - case 1: - case 2: - case 6: - (*nb_save_table_u)++; - //SaveCode << save_table_u->type << save_table_u->index << save_table_u->op1 << save_table_u->op2; - SaveCode.write(reinterpret_cast(&save_table_u->type), sizeof(save_table_u->type)); - SaveCode.write(reinterpret_cast(&save_table_u->index), sizeof(save_table_u->index)); - SaveCode.write(reinterpret_cast(&save_table_u->op1), sizeof(save_table_u->op1)); - SaveCode.write(reinterpret_cast(&save_table_u->op2), sizeof(save_table_u->op2)); -#ifdef PRINT_OUT - if(save_table_u->type == 1) - cout << "+u[" << save_table_u->index << "]=" << "u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "]\n"; - else if(save_table_u->type == 2) - cout << "+u[" << save_table_u->index << "]+=u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "]\n"; - else - cout << "+u[" << save_table_u->index << "]=1/(1-u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "])\n"; -#endif /**PRINT_OUT**/ - break; - case 5: - (*nb_save_table_u)++; - //SaveCode << save_table_u->type << save_table_u->index; - SaveCode.write(reinterpret_cast(&save_table_u->type), sizeof(save_table_u->type)); - SaveCode.write(reinterpret_cast(&save_table_u->index), sizeof(save_table_u->index)); -#ifdef PRINT_OUT - cout << "+push(u[" << save_table_u->index << "])\n"; -#endif /**PRINT_OUT**/ - break; - } - if(chk) - { - OK=(save_table_u!=last_table_u); - if(OK) - { - table_u = save_table_u->pNext; - free(save_table_u); - save_table_u = table_u; - } - } - else - { - table_u = save_table_u->pNext; - free(save_table_u); - save_table_u = table_u; - OK=(save_table_u!=last_table_u); - } - } -} - - -void -store_code(t_table_u **save_table_u, t_table_u *First_table_u, t_table_u *stop_table_u, t_table_y *table_y, t_table_y **s_table_y, int *nb_save_table_u, int *nb_save_table_y, int vertex_count) -{ - int i, k; - t_table_u *s_table_u, *table_u; - table_u = First_table_u->pNext; - s_table_u = (t_table_u*)malloc(sizeof(t_table_u) - 3 * sizeof(int)); - (*save_table_u) = s_table_u; - (*nb_save_table_u) = 0; - while(table_u != stop_table_u) - { - if((table_u->type > 7) || (table_u->type < 1)) - { - cout << "Error : table_u->type=" << int(table_u->type) << " *nb_save_table_u=" << *nb_save_table_u << "\n"; - exit(EXIT_FAILURE); - } - else - { - (*nb_save_table_u)++; -#ifdef PRINT_OUT - cout << "--table_u->type=" << int(table_u->type) << "\n"; -#endif - switch (table_u->type) - { - case 3: - case 7: - s_table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - sizeof(int)); - s_table_u = s_table_u->pNext; - s_table_u->type = table_u->type; - s_table_u->index = table_u->index; - s_table_u->op1 = table_u->op1; -#ifdef PRINT_OUT - if(s_table_u->type == 3) - cout << "st u[" << s_table_u->index << "]=-1/u[" << s_table_u->op1 << "]\n"; - else - cout << "st u[" << s_table_u->index << "]*=u[" << s_table_u->op1 << "]\n"; -#endif - break; - case 1: - case 2: - case 6: - s_table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) ); - s_table_u = s_table_u->pNext; - s_table_u->type = table_u->type; - s_table_u->index = table_u->index; - s_table_u->op1 = table_u->op1; - s_table_u->op2 = table_u->op2; -#ifdef PRINT_OUT - if(s_table_u->type == 1) - cout << "st u[" << s_table_u->index << "]=" << "u[" << s_table_u->op1 << "]*u[" << s_table_u->op2 << "]\n"; - else if(s_table_u->type == 2) - cout << "st u[" << s_table_u->index << "]+=u[" << s_table_u->op1 << "]*u[" << s_table_u->op2 << "]\n"; - else - cout << "st u[" << s_table_u->index << "]=1/(1-u[" << s_table_u->op1 << "]*u[" << s_table_u->op2 << "])\n"; -#endif /**PRINT_OUT**/ - break; - case 5: - s_table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - 2 * sizeof(int)); - s_table_u = s_table_u->pNext; - s_table_u->type = table_u->type; - s_table_u->index = table_u->index; -#ifdef PRINT_OUT - cout << "st push(u[" << s_table_u->index << "])\n"; -#endif /**PRINT_OUT**/ - break; - } - } - table_u = table_u->pNext; - } - s_table_u->pNext = NULL; -#ifdef PRINT_OUT - cout << "*nb_save_table_u=" << *nb_save_table_u << "\n"; - cout << "vertex_count=" << vertex_count << " nb_save_table_y=" << nb_save_table_y << "\n"; -#endif - (*s_table_y) = (t_table_y*)malloc((vertex_count - *nb_save_table_y) * sizeof(t_table_y)); - for(i = 0;i < vertex_count - *nb_save_table_y;i++) - { - (*s_table_y)[i].u_index = (int*)malloc(table_y[i + *nb_save_table_y].nb * sizeof(int)); - (*s_table_y)[i].y_index = (int*)malloc(table_y[i + *nb_save_table_y].nb * sizeof(int)); - /*(*s_table_y)[i].y_lead_lag = (int*)malloc(table_y[i + *nb_save_table_y].nb * sizeof(int));*/ - (*s_table_y)[i].index = table_y[i + *nb_save_table_y].index; - (*s_table_y)[i].nb = table_y[i + *nb_save_table_y].nb; - for(k = 0;k < table_y[i + *nb_save_table_y].nb;k++) - { - (*s_table_y)[i].u_index[k] = table_y[i + *nb_save_table_y].u_index[k]; - (*s_table_y)[i].y_index[k] = table_y[i + *nb_save_table_y].y_index[k]; - } - } - *nb_save_table_y = vertex_count - *nb_save_table_y; -} - - - -bool -SymbolicGaussElimination::Check_Regularity(t_table_u *first_u_blck, t_table_u *second_u_blck, t_table_u *third_u_blck) -{ - t_table_u *c_first_table_u, *c_second_table_u, *c_third_table_u; - bool OK, Regular; - int i=0; - c_first_table_u = first_u_blck->pNext; - c_second_table_u = second_u_blck->pNext; - c_third_table_u = third_u_blck->pNext; - OK = true; - Regular = true; - while(OK) - { - i++; -#ifdef PRINT_OUT - cout << "c_first_table_u=" << c_first_table_u << " c_second_table_u=" << c_second_table_u << " c_third_table_u=" << c_third_table_u << "\n"; -#endif /**PRINT_OUT**/ - if((c_first_table_u->type != c_second_table_u->type) || - (c_second_table_u->index + (c_second_table_u->index - c_first_table_u->index) != c_third_table_u->index)) - { - Regular = false; -#ifdef PRINT_OUT - cout << "c_first_table_u->type=" << (int)c_first_table_u->type << "?=c_second_table_u->type=" << (int)c_second_table_u->type << "\n"; - cout << "c_second_table_u->index+(c_second_table_u->index-c_first_table_u->index)=" << c_second_table_u->index + (c_second_table_u->index - c_first_table_u->index) - << "?=c_third_table_u->index=" << c_third_table_u->index << "\n"; -#endif - break; - } - if(c_first_table_u->type != 5) - { - if(c_second_table_u->op1 + (c_second_table_u->op1 - c_first_table_u->op1) != c_third_table_u->op1) - { -#ifdef PRINT_OUT - cout << "c_second_table_u->op1+(c_second_table_u->op1-c_first_table_u->op1)=" << c_second_table_u->op1 + (c_second_table_u->op1 - c_first_table_u->op1) - << "?=c_third_table_u->op1=" << c_third_table_u->op1 << "\n"; -#endif - Regular = false; - break; - } - if((c_first_table_u->type != 3) && (c_first_table_u->type != 7)) - { - if(c_second_table_u->op2 + (c_second_table_u->op2 - c_first_table_u->op2) != c_third_table_u->op2) - { -#ifdef PRINT_OUT - cout << "c_second_table_u->op2+(c_second_table_u->op2-c_first_table_u->op2)=" << c_second_table_u->op2 + (c_second_table_u->op2 - c_first_table_u->op2) - << "?=c_third_table_u->op2=" << c_third_table_u->op2 << "\n"; -#endif - Regular = false; - break; - } - } - } - if(c_first_table_u->pNext != second_u_blck->pNext) - { - c_second_table_u = c_second_table_u->pNext; - c_first_table_u = c_first_table_u->pNext; - c_third_table_u = c_third_table_u->pNext; - } - else - OK = 0; - } -#ifdef PRINT_OUT - cout << "Regular=" << Regular << " OK=" << OK << " i=" << i << "\n"; - cout << "Check_Regularity Regular=" << Regular << "\n"; -#endif - return (Regular); -} - - -//================================================================================== -t_table_u* -SymbolicGaussElimination::interpolation(t_model_graph* model_graph, t_table_y* table_y, int to_add, bool middle, int per) -{ - t_table_u *tmp_table_u=NULL, *old_table_u=NULL; - t_table_u *c_first_table_u, *c_second_table_u, *c_third_table_u; - int c_first_y_blck, c_second_y_blck; - int i, j, k, up_to=0, op_count, k1, k2; - bool OK, modify_u_count; - int cur_pos, nb_table_u=0; -#ifdef PRINT_OUT - cout << "in interpolation\n"; - cout << "--- u_count=" << u_count << "\n"; - print_Graph(model_graph); -#endif /**PRINT_OUT**/ - cur_pos=SaveCode.tellp(); - SaveCode.write(reinterpret_cast(&up_to), sizeof(up_to)); - if(middle) - { - up_to = 2; - middle_count_loop = 1 + per; - //middle_count_loop = 4; - } - else - up_to = 2; -#ifdef SAVE - if(middle) - { - middle_save_table_u = NULL; -#ifdef PRINT_OUT - cout << "up_to=" << up_to << " middle_count_loop=" << middle_count_loop << "\n"; -#endif - } -#endif /**SAVE**/ - /*Adding the Gaussian Elimination coefficients for the uncomputed periods*/ - /*table_u = tmp_table_u = (t_table_u*)malloc(sizeof(t_table_u) - 3 * sizeof(int));*/ - c_third_table_u = third_u_blck; - c_first_y_blck = first_y_blck; - c_second_y_blck = second_y_blck; - modify_u_count = true; - for(i = 1;i < up_to;i++) - { - c_first_table_u = first_u_blck->pNext; - c_second_table_u = second_u_blck->pNext; -#ifdef PRINT_OUT - cout << "c_first_table_u=" << c_first_table_u << " c_second_table_u=" << c_second_table_u << "\n"; -#endif - old_table_u = tmp_table_u; - op_count = 0; - OK = true; - while(OK) - { - if(c_first_table_u->type != c_second_table_u->type) - { - cout << "Error : Interpolation the two lists are not synchronized (middle=" << middle << ", op_count=" << op_count << ")\n"; - cout << "c_first_table_u=" << c_first_table_u << " c_second_table_u=" << c_second_table_u << "\n"; - cout << "y_kmin=" << y_kmin << " y_kmax=" << y_kmax << "\n"; - cout << "first_u_blck=" << first_u_blck << " second_u_blck=" << second_u_blck << "\n"; - cout << "c_first_table_u->type=" << int(c_first_table_u->type) << " c_second_table_u->type=" << int(c_second_table_u->type) << "\n"; - system("pause"); - exit(EXIT_FAILURE); - } - switch (c_first_table_u->type) - { - case 3: - case 7: - nb_table_u++; -#ifdef SAVE - if(i == 1) - { - if(middle) - { - nb_middle_save_table_u++; - } - else - { -#ifdef PRINT_OUT - cout << "!!!!!nb_first_save_table_u=" << nb_first_save_table_u << "\n"; -#endif - nb_first_save_table_u++; - } - SaveCode.write(reinterpret_cast(&c_first_table_u->type), sizeof(c_first_table_u->type)); - SaveCode.write(reinterpret_cast(&c_first_table_u->index), sizeof(c_first_table_u->index)); - SaveCode.write(reinterpret_cast(&c_first_table_u->op1), sizeof(c_first_table_u->op1)); -#ifdef PRINT_OUT - if(c_first_table_u->type == 3) - cout << ">u[" << c_first_table_u->index << "]=1/(1-u[" << c_first_table_u->op1 << "])\n"; //"]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; - else - cout << ">u[" << c_first_table_u->index << "]*=u[" << c_first_table_u->op1 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; -#endif - k=c_second_table_u->index - c_first_table_u->index; - SaveCode.write(reinterpret_cast(&k), sizeof(k)); - k1=c_second_table_u->op1- c_first_table_u->op1; - SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); -#ifdef PRINT_OUT - if(c_first_table_u->type == 3) - cout << ">i u[" << k << "]=1/(1-u[" << k1 << "])\n"; //"]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; - else - cout << ">i u[" << k << "]*=u[" << k1 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; -#endif - } -#endif /**SAVE**/ - break; - case 1: - case 2: - case 6: - nb_table_u++; -#ifdef SAVE - if((i == 1) && (middle)) - { - nb_middle_save_table_u++; - SaveCode.write(reinterpret_cast(&c_first_table_u->type), sizeof(c_first_table_u->type)); - SaveCode.write(reinterpret_cast(&c_first_table_u->index), sizeof(c_first_table_u->index)); - SaveCode.write(reinterpret_cast(&c_first_table_u->op1), sizeof(c_first_table_u->op1)); - SaveCode.write(reinterpret_cast(&c_first_table_u->op2), sizeof(c_first_table_u->op2)); -#ifdef PRINT_OUT - if((c_first_table_u->type == 1) && (middle)) - cout << ">u[" << c_first_table_u->index << "]=u[" << c_first_table_u->op1 << "]*u[" << c_first_table_u->op2 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; - else if(c_first_table_u->type == 2) - cout << ">u[" << c_first_table_u->index << "]+=u[" << c_first_table_u->op1 << "]*u[" << c_first_table_u->op2 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; - else - cout << ">u[" << c_first_table_u->index << "]=1/(1-u[" << c_first_table_u->op1 << "]*u[" << c_first_table_u->op2 << "])\n"; //"]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; -#endif //PRINT_OUT - - k=c_second_table_u->index - c_first_table_u->index; - SaveCode.write(reinterpret_cast(&k), sizeof(c_first_table_u->index)); - k1=c_second_table_u->op1- c_first_table_u->op1; - SaveCode.write(reinterpret_cast(&k1), sizeof(c_first_table_u->op1)); - k2=c_second_table_u->op2- c_first_table_u->op2; - SaveCode.write(reinterpret_cast(&k2), sizeof(c_first_table_u->op2)); -#ifdef PRINT_OUT - if((c_first_table_u->type == 1) && (middle)) - cout << ">u[" << k << "]=u[" << k1 << "]*u[" << k2 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; - else if(c_first_table_u->type == 2) - cout << ">u[" << k << "]+=u[" << k1 << "]*u[" << k2 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; - else - cout << ">u[" << k << "]=1/(1-u[" << k1 << "]*u[" << k2 << "])\n"; //"]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; -#endif //PRINT_OUT - } -#endif /**SAVE**/ - break; - case 5: - nb_table_u++; -#ifdef SAVE - if((i == 1) && (middle)) - { - nb_middle_save_table_u++; - SaveCode.write(reinterpret_cast(&c_first_table_u->type), sizeof(c_first_table_u->type)); - SaveCode.write(reinterpret_cast(&c_first_table_u->index), sizeof(c_first_table_u->index)); -#ifdef PRINT_OUT - cout << ">push(u[" << c_first_table_u->index << "])\n"; // "]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; -#endif //PRINT_OUT - k=c_second_table_u->index - c_first_table_u->index; - SaveCode.write(reinterpret_cast(&k), sizeof(c_first_table_u->index)); -#ifdef PRINT_OUT - cout << ">push(u[" << k << "])\n"; // "]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n"; -#endif //PRINT_OUT - - } -#endif /**SAVE**/ - break; - } - if((c_first_table_u->pNext != second_u_blck->pNext) && (c_first_table_u->pNext != NULL) && (c_second_table_u->pNext != NULL)) - { - tmp_table_u=c_first_table_u ; - c_first_table_u = c_first_table_u->pNext; - free(tmp_table_u); - tmp_table_u=c_second_table_u ; - c_second_table_u = c_second_table_u->pNext; - free(tmp_table_u); - } - else - { - if(c_first_table_u->pNext != second_u_blck->pNext) /*||(second_u_blck->pNext!=NULL))*/ - { - cout << "c_first_table_u->pNext=" << c_first_table_u->pNext << " second_u_blck->pNext=" << second_u_blck->pNext << "\n"; - cout << "Error not synchronize graph interpolation\n"; - exit(EXIT_FAILURE); - } - OK = 0; - } - op_count++; - } - SaveCode.flush(); - k=SaveCode.tellp(); - SaveCode.seekp(cur_pos); - SaveCode.write(reinterpret_cast(&nb_table_u), sizeof(nb_table_u)); - SaveCode.seekp(k); - if(middle) - { - first_y_blck = c_first_y_blck; - second_y_blck = c_second_y_blck; -#ifdef SAVE - if((i == 1) && (middle)) - { - middle_save_table_y = (t_table_y*)malloc((nb_endo) * sizeof(t_table_y)); - middle_save_i_table_y = (t_table_y*)malloc((nb_endo) * sizeof(t_table_y)); - nb_middle_save_table_y = nb_endo; - } -#endif /**SAVE**/ - for(j = 0;j < nb_endo;j++) - { -#ifdef SAVE - if(i == 1) - { - if(middle) - { - middle_save_i_table_y[j].u_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int)); - middle_save_i_table_y[j].y_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int)); - middle_save_i_table_y[j].index = table_y[second_y_blck].index - table_y[first_y_blck].index; - middle_save_i_table_y[j].nb = table_y[first_y_blck].nb; - for(k = 0;k < table_y[first_y_blck].nb;k++) - { - middle_save_i_table_y[j].u_index[k] = table_y[second_y_blck].u_index[k] - table_y[first_y_blck].u_index[k]; - middle_save_i_table_y[j].y_index[k] = table_y[second_y_blck].y_index[k] - table_y[first_y_blck].y_index[k]; - } - middle_save_table_y[j].u_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int)); - middle_save_table_y[j].y_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int)); - /*middle_save_table_y[j].y_lead_lag = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));*/ - middle_save_table_y[j].index = table_y[first_y_blck].index; - middle_save_table_y[j].nb = table_y[first_y_blck].nb; - for(k = 0;k < table_y[first_y_blck].nb;k++) - { - middle_save_table_y[j].u_index[k] = table_y[first_y_blck].u_index[k]; - middle_save_table_y[j].y_index[k] = table_y[first_y_blck].y_index[k]; - } - } - else - { - first_save_i_table_y[j].u_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int)); - first_save_i_table_y[j].y_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int)); - first_save_i_table_y[j].index = table_y[second_y_blck].index - table_y[first_y_blck].index; - first_save_i_table_y[j].nb = table_y[first_y_blck].nb; - for(k = 0;k < table_y[first_y_blck].nb;k++) - { - first_save_i_table_y[j].u_index[k] = table_y[second_y_blck].u_index[k] - table_y[first_y_blck].u_index[k]; - first_save_i_table_y[j].y_index[k] = table_y[second_y_blck].y_index[k] - table_y[first_y_blck].y_index[k]; - } - first_save_table_y[j].u_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int)); - first_save_table_y[j].y_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int)); - first_save_table_y[j].index = table_y[first_y_blck].index; - first_save_table_y[j].nb = table_y[first_y_blck].nb; - for(k = 0;k < table_y[first_y_blck].nb;k++) - { - first_save_table_y[j].u_index[k] = table_y[first_y_blck].u_index[k]; - first_save_table_y[j].y_index[k] = table_y[first_y_blck].y_index[k]; - } - } - } -#endif /**SAVE**/ - table_y[vertex_count].nb = table_y[first_y_blck].nb; - table_y[vertex_count].index = table_y[first_y_blck].index + (i - 1) * (table_y[second_y_blck].index - table_y[first_y_blck].index); -#ifdef PRINT_OUT - cout << ">y[" << table_y[vertex_count].index << "]="; -#endif /**PRINT_OUT**/ - for(k = 0;k < table_y[vertex_count].nb;k++) - { - table_y[vertex_count].u_index[k] = table_y[first_y_blck].u_index[k] + (i - 1) * (table_y[second_y_blck].u_index[k] - table_y[first_y_blck].u_index[k]); - table_y[vertex_count].y_index[k] = table_y[first_y_blck].y_index[k] + (i - 1) * (table_y[second_y_blck].y_index[k] - table_y[first_y_blck].y_index[k]); -#ifdef PRINT_OUT - cout << "u[" << table_y[vertex_count].u_index[k] << "]*y[" << table_y[vertex_count].y_index[k] << "]"; - if(k + 1 == table_y[vertex_count].nb) - cout << "\n"; - else - cout << "+"; -#endif /**PRINT_OUT**/ - } - vertex_count++; - first_y_blck++; - second_y_blck++; - } - } - } -#ifdef PRINT_OUT - cout << "+++ u_count=" << u_count << "\n"; - cout << "out interpolation\n"; -#endif - - return (old_table_u); -} - - -//================================================================================== -bool -SymbolicGaussElimination::Loop_Elimination(t_model_graph* model_graph) -{ - int i, j, pos, i1, i_per; - bool there_is_a_loop, go_on = true; - /* destroy the loop*/ - int per = 0; - there_is_a_loop = 0; - i_per = 0; - for(i = 0;i < model_graph->nb_vertices;i++) - { - if(nstacked) - { - if(model_graph->vertex[i].nb_in_degree_edges > 0) - { - if(!(i_per % nb_endo)) - { - per++; - if(per == 1) - { - first_u_blck = table_u; -#ifdef PRINT_OUT - cout << "first_u_blck=" << first_u_blck << "\n"; -#endif - } - if(per == 2) - { - pos_nb_first_save_table_u = nb_table_u; - second_u_blck = table_u; -#ifdef PRINT_OUT - cout << "second_u_blck=" << second_u_blck << "\n"; -#endif - } - if(per == 3) - { - go_on = false; - third_u_blck = table_u; -#ifdef PRINT_OUT - cout << "third_u_blck\n"; -#endif - } - } - i_per++; - } - } - pos = -1; - for(j = 0;((j < model_graph->vertex[i].nb_in_degree_edges) && (pos < 0));j++) - if(model_graph->vertex[i].in_degree_edge[j].index == i) - { - // at least one loop - pos = j; - } - if(pos >= 0) - { - there_is_a_loop = 1; -#ifdef DEBUGR - cout << "--------------------------------------------------------------- \n"; - cout << "Loop elimination on vertex " << model_graph->vertex[i].index << "\n"; -#endif /**DEBUGR**/ - if(nstacked) - { - /*Store the u_count associated to the loop to re-use it in case of - loop creation during vertex supression step*/ - loop_table_u_count[i] = model_graph->vertex[i].in_degree_edge[pos].u_count; - loop_table_vertex_index[i] = model_graph->vertex[i].index; - } - if(go_on) - { - nb_table_u++; - table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - sizeof(int)); - table_u = table_u->pNext; - table_u->pNext = NULL; - table_u->type = 3; - table_u->index = model_graph->vertex[i].in_degree_edge[pos].u_count; - table_u->op1 = model_graph->vertex[i].in_degree_edge[pos].u_count; -#ifdef PRINT_OUT - cout << "u[" << table_u->op1 << "]=1/(-u[" << table_u->op1 << "]) " << table_u << "\n"; // "]) (nb_free_u_list : " << nb_free_u_list << " u_count=" << u_count << ")\n"; -#endif /**PRINT_OUT**/ - for(j = 0;j < model_graph->vertex[i].nb_in_degree_edges;j++) - { - if(j != pos) - { - i1 = model_graph->vertex[i].in_degree_edge[j].index; -#ifdef DIRECT_COMPUTE - nb_table_u++; - table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - sizeof(int)); - table_u = table_u->pNext; - table_u->pNext = NULL; - table_u->type = 7; - table_u->index = model_graph->vertex[i].in_degree_edge[j].u_count; - table_u->op1 = model_graph->vertex[i].in_degree_edge[pos].u_count; -#ifdef PRINT_OUT - cout << "u[" << model_graph->vertex[i].in_degree_edge[j].u_count << "]*=u[" << table_u->op1 << "] " << table_u << "\n"; // "] (nb_free_u_list : " << nb_free_u_list << " u_count=" << u_count << ")\n"; -#endif /**PRINT_OUT**/ -#endif /**DIRECT_COMPUTE**/ - } - } - } - /*elimination of the loop*/ - /*in the out_degree*/ - pos = -1; - for(j = 0;((j < model_graph->vertex[i].nb_out_degree_edges) && (pos < 0));j++) - if(model_graph->vertex[i].out_degree_edge[j].index == i) - pos = j; - if(pos < 0) - { - cout << "Error: not symetric on a loop on vertex " << model_graph->vertex[i].index << "\n"; - print_Graph(model_graph); - system("pause"); - exit(EXIT_FAILURE); - } - for(j = pos + 1;j < model_graph->vertex[i].nb_out_degree_edges;j++) - { - model_graph->vertex[i].out_degree_edge[j - 1].index = model_graph->vertex[i].out_degree_edge[j].index; - model_graph->vertex[i].out_degree_edge[j - 1].u_count = model_graph->vertex[i].out_degree_edge[j].u_count; - } - model_graph->vertex[i].nb_out_degree_edges--; - /*in the in_degree*/ - pos = -1; - for(j = 0;((j < model_graph->vertex[i].nb_in_degree_edges) && (pos < 0));j++) - if(model_graph->vertex[i].in_degree_edge[j].index == i) - pos = j; - if(pos < 0) - { - cout << "Error: not symetric on a loop on vertex " << model_graph->vertex[i].index << "\n"; - print_Graph(model_graph); - system("pause"); - exit(EXIT_FAILURE); - } - for(j = pos + 1;j < model_graph->vertex[i].nb_in_degree_edges;j++) - { - model_graph->vertex[i].in_degree_edge[j - 1].index = model_graph->vertex[i].in_degree_edge[j].index; - model_graph->vertex[i].in_degree_edge[j - 1].u_count = model_graph->vertex[i].in_degree_edge[j].u_count; - } - model_graph->vertex[i].nb_in_degree_edges--; - } - } - return (there_is_a_loop); -} - -//================================================================================== -bool -SymbolicGaussElimination::Vertex_Elimination(t_model_graph* model_graph, int pos, bool* interpolate, int length_markowitz, bool dynamic) -{ - int i, j, k, min_edge, vertex_to_eliminate = -1, to_add, to_add_index, curr_u_count; - int i1, i2, j1, j2, i3, i4, i5, size; - int i_max, j_max; - bool OK; - t_vertex* lvertex = model_graph->vertex; - t_vertex* ilvertex; - int nb_new = 0, nb_free = 0; - int a_loop = 0; - char usi; -#ifdef DIRECT_COMPUTE - int in_y = 0; -#endif /**DIRECT_COMPUTE**/ - size = model_graph->nb_vertices; - min_edge = size * size + 1; - vertex_to_eliminate = -1; - // minimize the product of in and out degree of the remaining vertices (Markowitz criteria) - for(i = starting_vertex;i < starting_vertex + length_markowitz;i++) - if(((lvertex[i].nb_out_degree_edges*lvertex[i].nb_in_degree_edges) < min_edge) && (lvertex[i].nb_in_degree_edges > 0)) - { - min_edge = lvertex[i].nb_out_degree_edges * lvertex[i].nb_in_degree_edges; - vertex_to_eliminate = i; - } - OK = 0; - //cout << "save_direct=" << save_direct << "\n"; - if(vertex_to_eliminate >= 0) - { -#ifdef PRINT_OUT_1 - cout << "--------------------------------------------------------------- \n"; - cout << "Elimination of vertex " << lvertex[vertex_to_eliminate].index << " interpolate=" << *interpolate << "\n"; - cout << "min_edge=" << min_edge << " length_markowitz=" << length_markowitz << "\n"; - print_Graph(model_graph); -#endif /**PRINT_OUT**/ -#ifdef DIRECT_COMPUTE - table_y[vertex_count].index = lvertex[vertex_to_eliminate].index; -#ifdef PRINT_OUT - cout << "vertex_count=" << vertex_count << " size=" << size << " vertex_to_eliminate=" << vertex_to_eliminate << "\n"; -#endif -#endif /**DIRECT_COMPUTE**/ - ilvertex = &(lvertex[vertex_to_eliminate]); - i_max = ilvertex->nb_in_degree_edges; - j_max = ilvertex->nb_out_degree_edges; - for(i = 0;i < i_max;i++) - { - i1 = ilvertex->in_degree_edge[i].index; - i2 = ilvertex->in_degree_edge[i].u_count; -#ifdef DIRECT_COMPUTE - table_y[vertex_count].u_index[in_y] = i2; - table_y[vertex_count].y_index[in_y] = lvertex[i1].index; - in_y++; - nb_table_u++; - if(save_direct) - { - usi=5; - //SaveCode << (&usi) << i2; - SaveCode.write(reinterpret_cast(&usi), sizeof(usi)); - SaveCode.write(reinterpret_cast(&i2), sizeof(i2)); -#ifdef PRINT_OUT_1 - cout << "push(u[" << i2 << "]) \n"; -#endif - } - else - { - table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - 2 * sizeof(int)); - table_u = table_u->pNext; - table_u->pNext = NULL; - table_u->type = 5; - table_u->index = i2; - } - //cout << "ok0\n"; -#ifdef PRINT_OUT_1 - cout << "push(u[" << i2 << "]) (" << table_u << ")\n"; -#endif /**PRINT_OUT**/ -#endif /**DIRECT_COMPUTE**/ - for(j = 0;j < j_max;j++) - { - to_add = -9999999; - j1 = ilvertex->out_degree_edge[j].index; - j2 = ilvertex->out_degree_edge[j].u_count; - //cout << "ok1\n"; - /* Is there already an edge from vertex i1 to vertex j1? */ - for(k = 0;((k < lvertex[i1].nb_out_degree_edges) && (to_add < 0));k++) - if(lvertex[i1].out_degree_edge[k].index == j1) - { - /*yes*/ - to_add = lvertex[i1].out_degree_edge[k].u_count; - to_add_index = k; - } - if(to_add != -9999999) - { -#ifdef DEBUGR - cout << " modification of edge between vertices " << lvertex[i1].index << " and " << lvertex[j1].index << "\n"; -#endif /**DEBUGR**/ -#ifdef DIRECT_COMPUTE - if(save_direct) - { - nb_table_u++; - usi=2; - SaveCode.write(reinterpret_cast(&usi), sizeof(usi)); - SaveCode.write(reinterpret_cast(&to_add), sizeof(to_add)); - SaveCode.write(reinterpret_cast(&j2), sizeof(j2)); - SaveCode.write(reinterpret_cast(&i2), sizeof(i2)); - //SaveCode << (unsigned short int)(2) << to_add << j2 << i2; -#ifdef PRINT_OUT_1 - cout << "u[" << to_add << "]+=u[" << j2 << "]*u[" << i2 << "]; \n"; -#endif - } - else - { - nb_table_u++; - table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u)); - table_u = table_u->pNext; - table_u->pNext = NULL; - table_u->type = 2; - table_u->index = to_add; - table_u->op1 = j2; - table_u->op2 = i2; - } -#ifdef PRINT_OUT_1 - cout << "u[" << to_add << "]+=u[" << j2 << "]*u[" << i2 << "]; (" << table_u << ")\n"; -#endif /**PRINT_OUT**/ -#endif /**DIRECT_COMPUTE**/ - } - else - { - /*now modifie the in_degree edge from vertex_to_elminate to j1 */ - if(j1 == i1) - { - /* its a loop */ - /*Store it to operate after*/ - if(a_loop >= size) - { - cout << "Error : a_loop (" << a_loop << ") >= " << size << "\n"; - exit(EXIT_FAILURE); - } - s_j2[a_loop] = j2; - s_i2[a_loop] = i2; - s_i1[a_loop] = i1; - a_loop++; - } - else - { -#ifdef DEBUGR - cout << " creation of edge between vertices " << lvertex[i1].index << " and " << lvertex[j1].index << "\n"; -#endif /**DEBUGR**/ -#ifdef SYMPLIFY /*no reuse of free numbers, to be sure that there is no cyclcial behaviour in numbering*/ - curr_u_count = get_free_u_list(dynamic); -#else - curr_u_count = u_count; - u_count++; -#endif -#ifdef DIRECT_COMPUTE - nb_table_u++; - nb_new++; - if(save_direct) - { - usi=1; - SaveCode.write(reinterpret_cast(&usi), sizeof(usi)); - SaveCode.write(reinterpret_cast(&curr_u_count), sizeof(curr_u_count)); - SaveCode.write(reinterpret_cast(&j2), sizeof(j2)); - SaveCode.write(reinterpret_cast(&i2), sizeof(i2)); -#ifdef PRINT_OUT_1 - cout << "u[" << curr_u_count << "]=u[" << j2 << "]*u[" << i2 << "]; \n"; -#endif - } - else - { - table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u)); - table_u = table_u->pNext; - table_u->pNext = NULL; - table_u->type = 1; - table_u->index = curr_u_count; - table_u->op1 = j2; - table_u->op2 = i2; - } -#ifdef PRINT_OUT_1 - cout << "u[" << curr_u_count << "]=u[" << j2 << "]*u[" << i2 << "]; (" << table_u << ")\n"; -#endif /**PRINT_OUT**/ -#endif /**DIRECT_COMPUTE**/ - /*its not a loop so adding a new adge*/ - /*modify the in_degree edge from vertex_to_elminate to j1*/ -#ifdef SORTED - k = lvertex[j1].nb_in_degree_edges; - while((i1 < lvertex[j1].in_degree_edge[k - 1].index) && (k > 0)) - { - lvertex[j1].in_degree_edge[k].index = lvertex[j1].in_degree_edge[k - 1].index; - lvertex[j1].in_degree_edge[k].u_count = lvertex[j1].in_degree_edge[k - 1].u_count; - k--; - } - lvertex[j1].in_degree_edge[k].index = i1; - lvertex[j1].in_degree_edge[k].u_count = curr_u_count; - (lvertex[j1].nb_in_degree_edges)++; - /*Tested*/ - if(lvertex[j1].max_nb_in_degree_edges<=lvertex[j1].nb_in_degree_edges) - { - lvertex[j1].max_nb_in_degree_edges=lvertex[j1].nb_in_degree_edges; - t_edge *tmp_edge; - tmp_edge=lvertex[j1].in_degree_edge; - int taille=lvertex[j1].nb_in_degree_edges*sizeof(t_edge); - tmp_edge = (t_edge*)realloc(tmp_edge, taille); - lvertex[j1].in_degree_edge=tmp_edge; - } - /*EndTested*/ -#else - k = lvertex[j1].nb_in_degree_edges; - lvertex[j1].in_degree_edge[k].index = i1; - lvertex[j1].in_degree_edge[k].u_count = curr_u_count; - (lvertex[j1].nb_in_degree_edges)++; -#endif - /*now modify the out_degree edge from i1 to vertex_to_elminate */ -#ifdef SORTED - k = lvertex[i1].nb_out_degree_edges; - while((j1 < lvertex[i1].out_degree_edge[k - 1].index) && (k > 0)) - { - lvertex[i1].out_degree_edge[k].index = lvertex[i1].out_degree_edge[k - 1].index; - lvertex[i1].out_degree_edge[k].u_count = lvertex[i1].out_degree_edge[k - 1].u_count; - k--; - } - lvertex[i1].out_degree_edge[k].index = j1; - lvertex[i1].out_degree_edge[k].u_count = curr_u_count; - (lvertex[i1].nb_out_degree_edges)++; - /*Tested*/ - if(lvertex[j1].max_nb_out_degree_edgesin_degree_edge[i].index; -#ifdef DEBUGR - cout << " supress the out_degree edge from " << lvertex[i1].index << " to " << ilvertex->index << "\n"; -#endif /**DEBUGR**/ - to_add = -1; - for(k = 0;((k < lvertex[i1].nb_out_degree_edges) && (to_add < 0));k++) - if(lvertex[i1].out_degree_edge[k].index == vertex_to_eliminate) - to_add = k; - if(to_add < 0) - { - cout << "error: Model_Graph not correctly filled in out_degree (edge from " << lvertex[i1].index << " to " << lvertex[vertex_to_eliminate].index << ")\n"; - print_Graph(model_graph); - system("pause"); - exit(EXIT_FAILURE); - } -#ifdef SIMPLIFYS - nb_free++; - set_free_u_list(lvertex[i1].out_degree_edge[to_add].u_count); -#endif /**SIMPLIFY**/ -#ifdef SIMPLIFYT - if(simplification_allowed) - { - free_u_list1[nb_free_u_list1] = lvertex[i1].out_degree_edge[to_add].u_count; - if(nb_free_u_list1 > max_nb_free_u_list1) - max_nb_free_u_list1 = nb_free_u_list1; - nb_free_u_list1++; - } -#endif /**SIMPLIFY**/ - for(k = to_add + 1;k < lvertex[i1].nb_out_degree_edges;k++) - { - lvertex[i1].out_degree_edge[k - 1].index = lvertex[i1].out_degree_edge[k].index; - lvertex[i1].out_degree_edge[k - 1].u_count = lvertex[i1].out_degree_edge[k].u_count; - } - (lvertex[i1].nb_out_degree_edges)--; - } - for(j = 0;j < j_max;j++) - { - /*now supress the in_degree edge from vertex_to_elminate to j1 */ - j1 = ilvertex->out_degree_edge[j].index; -#ifdef DEBUGR - cout << " supress the in_degree edge from " << ilvertex->index << " to " << lvertex[j1].index << "\n"; -#endif /**DEBUGR**/ - to_add = -1; - for(k = 0;((k < lvertex[j1].nb_in_degree_edges) && (to_add < 0));k++) - if(lvertex[j1].in_degree_edge[k].index == vertex_to_eliminate) - to_add = k; - if(to_add < 0) - { - cout << "error: Model_Graph not correctly filled in in_degree (edge from " << lvertex[vertex_to_eliminate].index << " to " << lvertex[j1].index << ")\n"; - print_Graph(model_graph); - system("pause"); - exit(EXIT_FAILURE); - } -#ifdef SIMPLIFYS - set_free_u_list(lvertex[j1].in_degree_edge[to_add].u_count); - nb_free++; -#endif /**SIMPLIFY**/ -#ifdef SIMPLIFYT - if(simplification_allowed) - { - free_u_list1[nb_free_u_list1] = lvertex[j1].in_degree_edge[to_add].u_count; - if(nb_free_u_list1 > max_nb_free_u_list1) - max_nb_free_u_list1 = nb_free_u_list1; - nb_free_u_list1++; - } -#endif /**SIMPLIFY**/ - for(k = to_add + 1;k < lvertex[j1].nb_in_degree_edges;k++) - { - lvertex[j1].in_degree_edge[k - 1].index = lvertex[j1].in_degree_edge[k].index; - lvertex[j1].in_degree_edge[k - 1].u_count = lvertex[j1].in_degree_edge[k].u_count; - } - (lvertex[j1].nb_in_degree_edges)--; - } - if(a_loop) - { - for(i = 0;i < a_loop;i++) - { - i1 = s_i1[i]; - if(nstacked) - { - /*re-use the u_count index eliminated during the loop elmination*/ - k = 0; - while((k <= nb_loop_table) && (loop_table_vertex_index[k] != model_graph->vertex[i1].index)) - k++; - if(loop_table_vertex_index[k] != model_graph->vertex[i1].index) - { - cout << "Error can't find the loop associated to vertex " << model_graph->vertex[i1].index << "\n"; - k = 0; - while((k <= nb_loop_table) && (loop_table_vertex_index[k] != model_graph->vertex[i1].index)) - { - cout << "loop_table_vertex_index[" << k << "]=" << loop_table_vertex_index[k] << " =? model_graph->vertex[" << i1 << "].index=" << model_graph->vertex[i1].index << "\n"; - k++; - } - exit(EXIT_FAILURE); - } - curr_u_count = loop_table_u_count[k] ; - } - else - { -#ifdef SIMPLIFY - curr_u_count = get_free_u_list(dynamic); -#else /**SIMPLIFY**/ - curr_u_count = u_count; - u_count++; -#endif /**SIMPLIFY**/ - } - nb_table_u++; - if(save_direct) - { - usi=6; - SaveCode.write(reinterpret_cast(&usi), sizeof(usi)); - SaveCode.write(reinterpret_cast(&curr_u_count), sizeof(curr_u_count)); - SaveCode.write(reinterpret_cast(&s_i2[i]), sizeof(s_i2[i])); - SaveCode.write(reinterpret_cast(&s_j2[i]), sizeof(s_j2[i])); -#ifdef PRINT_OUT_1 - cout << "u[" << curr_u_count << "]=1/(1-u[" << s_i2[i] << "]*u[" << s_j2[i] << "]) ;\n"; -#endif - } - else - { - table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) ); - table_u = table_u->pNext; - table_u->pNext = NULL; - table_u->type = 6; - table_u->index = curr_u_count; - table_u->op1 = s_i2[i]; - table_u->op2 = s_j2[i]; - } - nb_new++; -#ifdef PRINT_OUT_1 - cout << "u[" << curr_u_count << "]=1/(1-u[" << s_i2[i] << "]*u[" << s_j2[i] << "]) (" << table_u << ");\n"; -#endif /**PRINT_OUT**/ - i5 = curr_u_count; - OK = 0; - for(k = 0;k < lvertex[i1].nb_in_degree_edges;k++) - { - i3 = lvertex[i1].in_degree_edge[k].index; - if(i3 != i1) - { - i4 = lvertex[i1].in_degree_edge[k].u_count; -#ifdef DIRECT_COMPUTE - nb_table_u++; - if(save_direct) - { - usi=7; - SaveCode.write(reinterpret_cast(&usi), sizeof(usi)); - SaveCode.write(reinterpret_cast(&i4), sizeof(i4)); - SaveCode.write(reinterpret_cast(&i5), sizeof(i5)); -#ifdef PRINT_OUT_1 - cout << "u[" << i4 << "]*=u[" << i5 << "]; \n"; -#endif - } - else - { - table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - sizeof(int)); - table_u = table_u->pNext; - table_u->pNext = NULL; - table_u->type = 7; - table_u->index = i4; - table_u->op1 = i5; - } -#ifdef PRINT_OUT_1 - cout << "u[" << i4 << "]*=u[" << i5 << "]; (" << table_u << ")\n"; -#endif /**PRINT_OUT**/ -#endif /**DIRECT_COMPUTE**/ - } - } -#ifdef SIMPLIFY - if(simplification_allowed) - { - set_free_u_list(i5); - nb_free++; - } -#endif /**SIMPLIFY**/ - } - } -#ifdef SIMPLIFYS - if(simplification_allowed) - for(i = 0;i < nb_free_u_list1;i++) - { - set_free_u_list(free_u_list1[i]); - nb_free++; - } -#endif /**SIMPLIFY**/ - /*Change all indexes*/ - free(ilvertex->in_degree_edge); - free(ilvertex->out_degree_edge); - for(i = vertex_to_eliminate + 1;i < size;i++) - lvertex[i - 1] = lvertex[i]; - (model_graph->nb_vertices)--; - OK = 0; - for(i = 0;i < model_graph->nb_vertices;i++) - { - for(j = 0;j < lvertex[i].nb_in_degree_edges;j++) - { - OK = 1; - if(lvertex[i].in_degree_edge[j].index > vertex_to_eliminate) - (lvertex[i].in_degree_edge[j].index)--; - } - for(j = 0;j < lvertex[i].nb_out_degree_edges;j++) - { - if(lvertex[i].out_degree_edge[j].index > vertex_to_eliminate) - (lvertex[i].out_degree_edge[j].index)--; - } - } - if(in_y>max_nb_table_y) - max_nb_table_y=in_y; -#ifdef DIRECT_COMPUTE - table_y[vertex_count].nb = in_y; -#ifdef PRINT_OUT_1 - cout << "y[" << table_y[vertex_count].index << "]="; - for(i = 0;i < table_y[vertex_count].nb;i++) - { - cout << "u[" << table_y[vertex_count].u_index[i] << "]*y[" << table_y[vertex_count].y_index[i] << "]"; - if(i + 1 < table_y[vertex_count].nb) - cout << "+"; - else - cout << "\n"; - } -#endif /**PRINT_OUT**/ - vertex_count++; -#endif /**DIRECT_COMPUTE**/ - } -#ifdef PRINT_OUT - cout << "nb_new=" << nb_new << " nb_free=" << nb_free << "\n"; - cout <<"end of vertex elimination\n"; -#endif - return (OK); -} - - -//================================================================================== -void -SymbolicGaussElimination::Gaussian_Elimination(t_model_graph* model_graph -#ifdef SAVE - , string file_name -#endif - , bool dynamic) -{ - int i, size, j, k, per; - bool OK, interpolate = false; - t_table_u *First_prologue_table_u; - int length_markowitz=0, per_next = 3; - bool try_to_interpolate=false; - int cur_pos=0; - int prologue_nb_table_u=0, first_nb_prologue_save_table_y=0; - int nb_first_u_blck, nb_second_u_blck=0, nb_third_u_blck=0; - int nb_first_y_blck, nb_second_y_blck=0, nb_third_y_blck=0; - - size = model_graph->nb_vertices; -#ifdef PRINT_OUT - cout << "going to open file file_open=" << file_open << " file_name={" << file_name << "}\n"; -#endif - if(file_open) - SaveCode.open((file_name + ".bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate ); - else - SaveCode.open((file_name + ".bin").c_str(), ios::out | ios::binary); - file_open = true; - if(!SaveCode.is_open()) - { - cout << "Error : Can't open file \"" << file_name << ".bin\" for writing\n"; - exit(EXIT_FAILURE); - } -#ifdef PRINT_OUT - print_Graph(model_graph); -#endif /**PRINT_OUT**/ - s_i1 = (int*)malloc(size * sizeof(int)); - s_i2 = (int*)malloc(size * sizeof(int)); - s_j2 = (int*)malloc(size * sizeof(int)); -#ifdef SIMPLIFY - simplification_allowed = 1; - MAX_FREE_U_LIST = size * /*size *//*12*/50; - free_u_list = (int*)malloc(MAX_FREE_U_LIST* sizeof(int)); -#ifdef PRINT_OUT - cout << "free_u_list=" << free_u_list << " length=" << ceil((double)4 / (double)2*(double)u_count)*sizeof(int) << "\n"; - cout << "free_u_list length0=" << (*((short int*)&(*(((char*)free_u_list) - 4))) & ~(3)) << "\n"; - cout << "free_u_list1=(int*)malloc(" << u_count << ");\n"; -#endif - free_u_list1 = (int*)malloc(u_count * sizeof(int)); -#ifdef PRINT_OUT - cout << "after free alloc u_count=" << u_count << "\n"; -#endif -#endif /**SIMPLIFY**/ -#ifdef PRINT_OUT - cout << "periods=" << periods << "\n"; -#endif - if(dynamic) - u_count = (periods + y_kmin + y_kmax) * u_count / (y_kmin + y_kmax + 2); - if(nstacked) - { - nb_loop_table = model_graph->nb_vertices; - loop_table_u_count = (int*)malloc(nb_loop_table * sizeof(int)); - loop_table_vertex_index = (int*)malloc(nb_loop_table * sizeof(int)); - } -#ifdef PRINT_OUT - cout << "after nb_loop_table=" << nb_loop_table << "\n"; - system("pause"); -#endif - - SaveCode.write(reinterpret_cast(&nb_endo), sizeof(nb_endo)); - SaveCode.write(reinterpret_cast(&u_count), sizeof(u_count)); - SaveCode.write(reinterpret_cast(&u_count_init), sizeof(u_count_init)); - - -#ifdef PRINT_OUT - cout << "u_count=" << u_count << "\n"; - //We start with the elimination of all loops in the graph - cout << "going to loop elminate\n"; -#endif - save_direct=false; - if(nstacked) - save_direct=false; - else - { - i=0; - SaveCode.write(reinterpret_cast(&i), sizeof(i)); - SaveCode.write(reinterpret_cast(&i), sizeof(i)); - SaveCode.write(reinterpret_cast(&i), sizeof(i)); - SaveCode.write(reinterpret_cast(&i), sizeof(i)); - } - if(Loop_Elimination(model_graph)) - { -#ifdef PRINT_OUT - cout << "->table_u=" << table_u << "\n"; -#endif - if(nstacked) - { - //if there are some loops eliminated, we have to update it for the next blocks - // So we call the interpolation procedure - //cout << "nb_first_table_u="; - interpolation(model_graph, table_y, 1, false, 0); - } -#ifdef PRINT_OUT - else - { - cout << "nb_first_save_table_u=" << nb_first_save_table_u << " nb_table_u=" << nb_table_u << "\n"; - } -#endif -#ifdef PRINT_OUT - cout << "->table_u=" << table_u << "\n"; -#endif - } - if(nstacked) - { - first_nb_prologue_save_table_y = vertex_count; - nb_prologue_save_table_y=0; - First_prologue_table_u = table_u; - cur_pos=SaveCode.tellp(); - prologue_nb_table_u=nb_table_u; - i=0; - SaveCode.write(reinterpret_cast(&i), sizeof(i)); -#ifdef PRINT_OUT - cout << "--> fixe prologue\n"; -#endif - } -#ifdef PRINT_OUT -#ifdef SIMPLIFY - cout << "1 free_u_list length0=" << (*((short int*)&(*(((char*)free_u_list) - 4))) & ~(3)) << "\n"; -#endif /**SIMPLIFY**/ -#endif -#ifdef DEBUGR - //Check the graph -#ifdef PRINT_OUT - cout << "going to Check Graph\n"; -#endif - Check_Graph(model_graph); -#endif - per = 0; - OK = true; -#ifdef PRINT_OUT - cout << "nb_endo=" << nb_endo << " y_kmin=" << y_kmin << " y_kmax=" << y_kmax << " size=" << size << "\n"; -#endif - int pos = 0; - j = 0; - if(nstacked) - try_to_interpolate = true; -#ifdef PRINT_OUT - cout << "try_to_interpolate=" << try_to_interpolate << "\n"; - cout << "y_kmin=" << y_kmin << " y_kmax=" << y_kmax << "\n"; - cout << "(0) nb_table_u=" << nb_table_u << "\n"; -#endif - if(nstacked) - save_direct=true; - while(OK) - { -#ifdef DEBUGR - Check_Graph(model_graph); -#endif /**DEBUGR**/ - if(!(j % nb_endo)) - length_markowitz = nb_endo; - else - length_markowitz--; - if(nstacked) - { -#ifdef PRINT_OUT - cout << "try_to_interpolate=" << try_to_interpolate << "\n"; -#endif - if(try_to_interpolate) - { -#ifdef PRINT_OUT - cout << j << " % " << nb_endo << " = " << (j % nb_endo) << "\n"; -#endif - if(!(j % nb_endo)) - { - per++; - if(per == y_kmin + 1) - { - nb_prologue_save_table_y=vertex_count-first_nb_prologue_save_table_y; - prologue_save_table_y=(t_table_y*)malloc(nb_prologue_save_table_y*sizeof(*prologue_save_table_y)); - for(i=first_nb_prologue_save_table_y;i(&i), sizeof(i)); - SaveCode.seekp(k); - if(nstacked) - save_direct=false; - first_u_blck = table_u; - nb_first_u_blck= nb_table_u; - nb_first_y_blck= vertex_count; - first_y_blck = vertex_count; -#ifdef PRINT_OUT - cout << "(1) nb_table_u=" << nb_table_u << "\n"; - cout << "first_u_blck=" << first_u_blck << "\n"; - system("pause"); -#endif - } - else if(per == y_kmin + 2) - { - second_u_blck = table_u; - nb_second_u_blck= nb_table_u; - second_y_blck = vertex_count; - nb_second_y_blck= vertex_count; -#ifdef PRINT_OUT - cout << "(2) nb_table_u=" << nb_table_u << "\n"; - cout << "second_u_blck=" << second_u_blck << "\n"; - system("pause"); -#endif - } - else if(per >= y_kmin + 3) - { - third_u_blck = table_u; - nb_third_u_blck= nb_table_u; - third_y_blck = vertex_count; - nb_third_y_blck= vertex_count; -#ifdef PRINT_OUT - cout << "(3) nb_table_u=" << nb_table_u << "\n"; - cout << "third_u_blck=" << third_u_blck << "\n"; - system("pause"); -#endif - } - } - } - OK = Vertex_Elimination(model_graph, pos, &interpolate, length_markowitz, dynamic); - } - else - { -#ifdef PRINT_OUT - cout << "j=" << j << "\n"; -#endif - OK = Vertex_Elimination(model_graph, pos, &interpolate, length_markowitz, dynamic); - } -#ifdef SIMPLIFY -#endif /**SIMPLIFY**/ - j++; -#ifdef SIMPLIFY - if(!(j % nb_endo)) - nb_free_u_list = 0; -#endif /**SIMPLIFY**/ - if(nstacked) - { -#ifdef PRINT_OUT - cout << "try_to_interpolate=" << try_to_interpolate << "\n"; -#endif - if(try_to_interpolate) - { -#ifdef PRINT_OUT - cout << j << " % " << nb_endo << " = " << (j % nb_endo) << "\n"; -#endif - if(!((j) % nb_endo)) - { -#ifdef PRINT_OUT - cout << "per (" << per << ") >= y_kmin (" << y_kmin << ")+3\n"; -#endif - if( per >= y_kmin + 3) - { - // pctimer_t t1, t2; -#ifdef PRINT_OUT - cout << "bef check regularity\n"; - system("pause"); -#endif - if(Check_Regularity(first_u_blck, second_u_blck, third_u_blck)) - { -#ifdef PRINT_OUT - cout << "af check regularity OK \n"; - system("pause"); - cout << "nb_first_save_table_u=" << nb_first_save_table_u << "\n"; -#endif -#ifdef PRINT_OUT - cout << "table_u=" << table_u << " table_u->pNext=" << table_u->pNext << "\n"; - cout << "(3) nb_table_u=" << nb_table_u << "\n"; -#endif - // t1 = pctimer(); - - k=SaveCode.tellp(); -#ifdef PRINT_OUT - cout << "middle_count_loop= " << middle_count_loop << "\n"; -#endif - i=0; - SaveCode.write(reinterpret_cast(&i), sizeof(i)); - table_u = interpolation(model_graph, table_y, per_next /*+1*/ + 2, true, per - y_kmin - 3); -#ifdef PRINT_OUT - cout << "middle_count_loop= " << middle_count_loop << "\n"; -#endif - i=SaveCode.tellp(); - SaveCode.seekp(k); - SaveCode.write(reinterpret_cast(&middle_count_loop), sizeof(middle_count_loop)); - SaveCode.seekp(i); - i=0; - SaveCode.write(reinterpret_cast(&i), sizeof(i)); //last_u - OK = false; - vertex_count -= nb_endo; - // t2 = pctimer(); -#ifdef PRINT_OUT - // cout << "interpolate=" << 1000*(t2 - t1) << "\n"; - cout << "table_u=" << table_u << " table_u->pNext=" << table_u->pNext << "\n"; -#endif -#ifdef SAVE - last_save_table_u = table_u; - nb_last_save_table_y = vertex_count; -#endif /**SAVE**/ - } - else - { - nb_prologue_save_table_y=nb_second_y_blck-first_nb_prologue_save_table_y; - prologue_save_table_y=(t_table_y*)malloc(nb_prologue_save_table_y*sizeof(*prologue_save_table_y)); - for(i=first_nb_prologue_save_table_y;i(&i), sizeof(i)); - SaveCode.seekp(k); - write_to_file_table_u_b(first_u_blck,second_u_blck, &nb_prologue_save_table_u,true); - //cout << "nb_prologue_save_table_u(1)=" << nb_prologue_save_table_u << "\n"; - nb_first_u_blck=nb_second_u_blck; - nb_second_u_blck=nb_third_u_blck+1; - nb_first_y_blck=nb_second_y_blck; - nb_second_y_blck=nb_third_y_blck; - first_u_blck = second_u_blck; - second_u_blck = third_u_blck; - first_y_blck = second_y_blck; - second_y_blck = third_y_blck; - } - } - } -#ifdef DEBUGR - cout << 100*j / size << "%\n"; -#endif /**DEBUGR**/ - } - } - } -#ifdef PRINT_OUT - cout << "nstacked=" << nstacked << "\n"; - cout << "nb_first_save_table_y=" << nb_first_save_table_y << "\n"; - cout << "nb_prologue_save_table_y=" << nb_prologue_save_table_y << "\n"; - cout << "nb_middle_save_table_y=" << nb_middle_save_table_y << "\n"; - cout << "nb_last_save_table_y=" << nb_last_save_table_y << "\n"; -#endif - if((nstacked)&&(!nb_last_save_table_y)) - { - cout << "not synchronized per=" << per << "\n"; - exit(EXIT_FAILURE); - } - -#ifdef PRINT_OUT - cout << "end vertex supress\n"; -#endif -#ifdef SAVE - if(!nstacked) - table_u->pNext = NULL; - if(nstacked) - { - table_u = NULL; - last_save_table_u = NULL; - nb_last_save_table_y = 0; - OK = true; - } - else - { - SaveCode.write(reinterpret_cast(&nb_table_u), sizeof(nb_table_u)); - write_to_file_table_u_b(First_table_u->pNext, table_u->pNext, &nb_last_save_table_u, false ); - nb_last_save_table_y = vertex_count; - last_save_table_y=(t_table_y*)malloc(nb_last_save_table_y*sizeof(*last_save_table_y)); - for(i=0;iwrite prologue\n"; -#endif -#ifdef PRINT_OUT - cout << "-->write first\n"; -#endif -#ifdef PRINT_OUT - cout << "middle_count_loop=" << middle_count_loop << "\n"; - cout << "-->write middle\n"; -#endif -#ifdef PRINT_OUT - cout << "-->write last\n"; -#endif - nb_last_save_table_u = i; -#ifdef PRINT_OUT - cout << "nb_prologue_save_table_y=" << nb_prologue_save_table_y << "\n"; -#endif - write_to_file_table_y( prologue_save_table_y, NULL, nb_prologue_save_table_y, 0); -#ifdef PRINT_OUT - cout << "nb_first_save_table_y=" << nb_first_save_table_y << "\n"; -#endif - write_to_file_table_y( first_save_table_y, NULL, nb_first_save_table_y, 0); -#ifdef PRINT_OUT - cout << "nb_middle_save_table_y=" << nb_first_save_table_y << "\n"; -#endif - write_to_file_table_y( middle_save_table_y, middle_save_i_table_y, nb_middle_save_table_y, nb_middle_save_table_y); -#ifdef PRINT_OUT - cout << "//// last_save_table_y ===\n"; - cout << "nb_last_save_table_y=" << nb_last_save_table_y << "\n"; -#endif - write_to_file_table_y( last_save_table_y, NULL, nb_last_save_table_y, 0); - SaveCode.close(); -#endif /**SAVE**/ -#ifdef SIMPLIFY -#ifdef PRINT_OUT - cout << "try_to_interpolate=" << try_to_interpolate << "\n"; -#endif - free(free_u_list); - free(free_u_list1); -#endif /**SIMPLIFY**/ - free(s_i1); - free(s_i2); - free(s_j2); -} - -void -SymbolicGaussElimination::file_is_open() -{ - file_open=true; - //file_is_open1(); -} - -void -SymbolicGaussElimination::SGE_compute(Model_Block *ModelBlock, int blck, bool dynamic, string file_name, int endo_nbr) -{ - t_model_graph *model_graph; - int block_u_count, nb_table_y; - // pctimer_t t1, t2; - int i; - int mean_var_in_equation; - - init_glb(); - model_graph = (t_model_graph*)malloc(sizeof(*model_graph)); - nstacked = dynamic; -#ifdef PRINT_OUT - periods = ModelBlock->Periods; - cout << "nstacked=" << nstacked << "\n"; - // t1 = pctimer(); - cout << "periods=" << periods << "\n"; -#endif - - u_count = ModelBlock_Graph(ModelBlock, blck, dynamic, model_graph, endo_nbr, &block_u_count, &starting_vertex, &periods, &nb_table_y, &mean_var_in_equation); - if(dynamic) - { - cout << "Mean endogenous variable per equation: " << mean_var_in_equation << ", density indicator=" << double(mean_var_in_equation)/endo_nbr*100 << "%\n"; - cout << "Coding the model..."; - } - Time = periods; -#ifdef PRINT_OUT - // t2 = pctimer(); - // cout << "/* ModelBlock_Graph : " << 1000*(t2 - t1) << " milliseconds u_count : " << u_count << "*/\n"; -#endif - int size = model_graph->nb_vertices; - nb_endo = ModelBlock->Block_List[blck].Size; - y_kmin = ModelBlock->Block_List[blck].Max_Lag; - y_kmax = ModelBlock->Block_List[blck].Max_Lead; - periods = ModelBlock->Periods; - if(periods<=3) - nstacked=false; - //cout << "periods=" << periods << " y_kmin=" << y_kmin << " y_kmax=" << y_kmax << "\n"; - u_count_init = u_count; -#ifdef PRINT_OUT - cout << "size : " << size << "\n"; -#endif - table_y = (t_table_y*)malloc((size + 1) * sizeof(t_table_y)); - for(i = 0;i <= size;i++) - { - table_y[i].u_index = (int*)malloc(nb_table_y * sizeof(int)); - table_y[i].y_index = (int*)malloc(nb_table_y * sizeof(int)); - table_y[i].index = i; - table_y[i].nb = 0; - } - table_u = (t_table_u*)malloc(sizeof(*table_u)-3*sizeof(int)); - table_u->pNext = NULL; - First_table_u = table_u; -#ifdef PRINT_OUT - cout << "oł en est-on ? nb_table_y=" << nb_table_y << "\n"; - system("pause"); - cout << "u_count_init=" << u_count_init << "\n"; - cout << "table_u=" << table_u << "\n"; - t1 = pctimer(); -#endif - Gaussian_Elimination(model_graph, file_name, dynamic); - free_model_graph(model_graph); -#ifdef PRINT_OUT - cout << "max_nb_table_y=" << max_nb_table_y << "\n"; - cout << "max_nb_in_degree_edges=" << max_nb_in_degree_edges << "\n"; - cout << "max_nb_out_degree_edges=" << max_nb_out_degree_edges << "\n"; - t2 = pctimer(); - cout << "/* Gaussian Elimination : " << 1000*(t2 - t1) << " milliseconds u_count : " << u_count << "*/\n"; -#endif -} - diff --git a/preprocessor/include/ExprNode.hh b/preprocessor/include/ExprNode.hh index 9fb52b9df..29cc50aaa 100644 --- a/preprocessor/include/ExprNode.hh +++ b/preprocessor/include/ExprNode.hh @@ -159,7 +159,7 @@ public: }; virtual double eval(const eval_context_type &eval_context) const throw (EvalException) = 0; - virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const = 0; + virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const = 0; }; //! Object used to compare two nodes (using their indexes) @@ -186,7 +186,7 @@ public: virtual void collectExogenous(set > &result) const; virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); - virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; + virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; }; //! Symbol or variable node @@ -214,7 +214,7 @@ public: map_idx_type &map_idx) const; virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); - virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; + virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; }; //! Unary operator node @@ -243,7 +243,7 @@ public: virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); - virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; + virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; }; //! Binary operator node @@ -273,7 +273,7 @@ public: virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); - virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; + virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; virtual NodeID get_arg1() { return(arg1);}; virtual NodeID get_arg2() { return(arg2);}; }; @@ -310,7 +310,7 @@ public: virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); - virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; + virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; }; //! Unknown function node @@ -337,7 +337,7 @@ public: virtual void collectExogenous(set > &result) const; virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); - virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; + virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const; }; //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model diff --git a/preprocessor/include/ModFile.hh b/preprocessor/include/ModFile.hh index 3e0ed90c6..f673a15c9 100644 --- a/preprocessor/include/ModFile.hh +++ b/preprocessor/include/ModFile.hh @@ -23,6 +23,8 @@ using namespace std; #include +#include + #include "SymbolTable.hh" #include "NumericalConstants.hh" diff --git a/preprocessor/include/ModelTree.hh b/preprocessor/include/ModelTree.hh index 3f4ff0f83..9a3b14c34 100644 --- a/preprocessor/include/ModelTree.hh +++ b/preprocessor/include/ModelTree.hh @@ -32,7 +32,6 @@ using namespace std; #include "NumericalConstants.hh" #include "DataTree.hh" #include "BlockTriangular.hh" -#include "SymbolGaussElim.hh" //! The three in which ModelTree can work enum ModelTreeMode @@ -85,7 +84,7 @@ private: //! Write derivative of an equation w.r. to a variable void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, SymbolType type) const; //! Write derivative code of an equation w.r. to a variable - void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type map_idx) const; + void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const; //! Computes temporary terms void computeTemporaryTerms(int order); void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock); @@ -109,7 +108,7 @@ private: //! Writes the Block reordred structure of the static model in M output void writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const string &static_basename) const; //! Writes the code of the Block reordred structure of the model in virtual machine bytecode - void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type) const; + void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type, map_idx_type map_idx) const; //! Writes static model file (Matlab version) void writeStaticMFile(const string &static_basename) const; //! Writes static model file (C version) @@ -173,7 +172,7 @@ public: BlockTriangular block_triangular; //! Adds informations for simulation in a binary file void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, - const int &num, int &u_count_int, bool &file_open) const; + const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const; //! Returns the number of equations in the model int equation_number() const; }; diff --git a/tests/ferhat/ls2003.mod b/tests/ferhat/ls2003.mod index ef4366157..bfd229252 100644 --- a/tests/ferhat/ls2003.mod +++ b/tests/ferhat/ls2003.mod @@ -17,8 +17,8 @@ rho_ys = 0.9; rho_pies = 0.7; -//model(sparse_dll,cutoff=1e-17); -model(sparse); +model(sparse_dll); +//model(sparse); //model; y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1); pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s; @@ -84,7 +84,7 @@ options_.maxit_=100; steady; model_info; -check; +//check; shocks; var e_q; @@ -99,5 +99,5 @@ rplot ww; rplot A; rplot pie; -stoch_simul(periods=200,order=1); +//stoch_simul(periods=200,order=1); diff --git a/tests/ferhat/ramst.mod b/tests/ferhat/ramst.mod index 32c088065..fdf612aa6 100644 --- a/tests/ferhat/ramst.mod +++ b/tests/ferhat/ramst.mod @@ -10,8 +10,8 @@ aa=0.5; //model(sparse); -//model(sparse_dll); -model; +model(sparse_dll); +//model; c + k - aa*x*k(-1)^alph - (1-delt)*k(-1); c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam); end;