From d2dfda62d6d44631a6c8efdad6478fbc84dee8ef Mon Sep 17 00:00:00 2001 From: ferhat Date: Fri, 10 Jul 2009 15:10:11 +0000 Subject: [PATCH] - Correction of several bugs - Reduction of the SparseMatrix.cc code git-svn-id: https://www.dynare.org/svn/dynare/trunk@2835 ac1d8469-bf42-47a9-8791-bf33cf982152 --- mex/sources/simulate/Interpreter.cc | 295 +++---- mex/sources/simulate/Mem_Mngr.cc | 9 +- mex/sources/simulate/SparseMatrix.cc | 1069 ++------------------------ mex/sources/simulate/simulate.cc | 207 +++-- 4 files changed, 357 insertions(+), 1223 deletions(-) diff --git a/mex/sources/simulate/Interpreter.cc b/mex/sources/simulate/Interpreter.cc index 3db1879a4..bb7205ec4 100644 --- a/mex/sources/simulate/Interpreter.cc +++ b/mex/sources/simulate/Interpreter.cc @@ -40,9 +40,11 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub solve_tolf=solve_tolf_arg; size_of_direction=size_of_direction_arg; slowc=slowc_arg; + slowc_save = slowc; y_decal=y_decal_arg; markowitz_c=markowitz_c_arg; filename=filename_arg; + T=NULL; //GaussSeidel=true; } @@ -52,6 +54,7 @@ Interpreter::pow1(double a, double b) double r=pow_(a,b); if (isnan(r) || isinf(r)) { + //mexPrintf("pow(%f, %f)=%f\n",a, b, r); max_res=res1=res2=r; return(r); } @@ -65,7 +68,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ { int var, lag, op; double v1, v2; - char cc; + char/*uint8_t*/ cc; bool go_on=true; double *ll; while (go_on) @@ -85,7 +88,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ switch (get_code_char) { case eParameter : - var=get_code_int + var=get_code_int; #ifdef DEBUGC if(Block_Count==2) { @@ -96,8 +99,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ Stack.push(params[var]); break; case eEndogenous : - var=get_code_int - lag=get_code_int + var=get_code_int; + lag=get_code_int; #ifdef DEBUGC if(Block_Count==2) { @@ -110,8 +113,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ #ifdef DEBUGC mexPrintf("Exogenous\n"); #endif - var=get_code_int - lag=get_code_int + var=get_code_int; + lag=get_code_int; #ifdef DEBUGC if(Block_Count==2) { @@ -125,8 +128,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ mexPrintf("ExogenousDet\n"); #endif - var=get_code_int - lag=get_code_int + var=get_code_int; + lag=get_code_int; #ifdef DEBUGC mexPrintf(" x(det)[%d]=%f\n",it_+lag+var*nb_row_xd,x[it_+lag+var*nb_row_xd]); mexEvalString("drawnow;"); @@ -134,12 +137,12 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ Stack.push(x[it_+lag+var*nb_row_xd]); break; default: - mexPrintf("Unknown vraibale type\n"); + mexPrintf("Unknown variable type\n"); } break; case FLDT : //load a temporary variable in the processor - var=get_code_int + var=get_code_int; #ifdef DEBUGC mexPrintf("FLDT %d\n",var); mexEvalString("drawnow;"); @@ -152,7 +155,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ mexPrintf("FLDU\n"); mexEvalString("drawnow;"); #endif - var=get_code_int + var=get_code_int; var+=Per_u_; Stack.push(u[var]); break; @@ -162,7 +165,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ mexPrintf("FLDR\n"); mexEvalString("drawnow;"); #endif - var=get_code_int + var=get_code_int; Stack.push(r[var]); break; case FLDZ : @@ -177,7 +180,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ //load a numerical constant in the processor /*asm("fld\n\t" "fstp %%st" : "=t" (ll) : "0" ((double)(*Code)));*/ - ll=get_code_pdouble + ll=get_code_pdouble; #ifdef DEBUGC mexPrintf("FLDC %f\n",*ll); mexEvalString("drawnow;"); @@ -193,7 +196,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ switch (get_code_char) { case eParameter : - var=get_code_int + var=get_code_int; params[var] = Stack.top(); #ifdef DEBUGC mexPrintf("FSTP params[%d]=%f\n", var, params[var]); @@ -205,8 +208,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ #ifdef DEBUGC mexPrintf("FSTP Endogenous\n"); #endif - var=get_code_int - lag=get_code_int + var=get_code_int; + lag=get_code_int; #ifdef DEBUGC //mexPrintf("y[%d(it_=%d, lag=%d, y_size=%d, var=%d)](%d)=",(it_+lag)*y_size+var,it_, lag, y_size, var, Stack.size()); mexPrintf("y[it_=%d, lag=%d, y_size=%d, var=%d, block=%d)]=",it_, lag, y_size, var+1, Block_Count+1); @@ -223,9 +226,9 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ #ifdef DEBUGC mexPrintf("Exogenous\n"); #endif - var=get_code_int - var=get_code_int - lag=get_code_int + var=get_code_int; + var=get_code_int; + lag=get_code_int; x[it_+lag+var*nb_row_x] = Stack.top(); Stack.pop(); break; @@ -233,10 +236,9 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ #ifdef DEBUGC mexPrintf("ExogenousDet\n"); #endif - var=get_code_int - - var=get_code_int - lag=get_code_int + var=get_code_int; + var=get_code_int; + lag=get_code_int; x[it_+lag+var*nb_row_xd] = Stack.top(); Stack.pop(); break; @@ -246,7 +248,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ break; case FSTPT : //load a temporary variable in the processor - var=get_code_int + var=get_code_int; #ifdef DEBUGC mexPrintf("FSTPT T[(var=%d, it_=%d, periods=%d, y_kmin=%d, y_kmax=%d)%d]=", var, it_, periods, y_kmin, y_kmax, var*(periods+y_kmin+y_kmax)+it_); mexEvalString("drawnow;"); @@ -260,7 +262,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ break; case FSTPU : //load u variable in the processor - var=get_code_int + var=get_code_int; var+=Per_u_; #ifdef DEBUGC mexPrintf("FSTPU u[%d]",var); @@ -275,7 +277,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ break; case FSTPR : //load u variable in the processor - var=get_code_int + var=get_code_int; r[var] = Stack.top(); #ifdef DEBUGC mexPrintf("FSTPR residual[%d]=%f\n",var,r[var]); @@ -285,7 +287,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ break; case FSTPG : //load u variable in the processor - var=get_code_int + var=get_code_int; g1[var] = Stack.top(); #ifdef DEBUGC mexPrintf("FSTPG g1[%d)=%f\n",var,g1[var]); @@ -298,8 +300,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ mexPrintf("FBINARY\n"); mexEvalString("drawnow;"); #endif - op=get_code_int - v2=Stack.top(); + op=get_code_int; + v2=Stack.top(); Stack.pop(); v1=Stack.top(); Stack.pop(); @@ -376,11 +378,11 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ #endif break; case oPower: - Stack.push(pow1(v1, v2)); #ifdef DEBUGC mexPrintf("pow(%f, %f)\n",v1,v2); mexEvalString("drawnow;"); #endif + Stack.push(pow1(v1, v2)); break; case oMax: Stack.push(max(v1, v2)); @@ -407,7 +409,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ mexPrintf("FUNARY\n"); mexEvalString("drawnow;"); #endif - op=get_code_int + op=get_code_int; v1=Stack.top(); Stack.pop(); switch (op) @@ -503,7 +505,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ mexEvalString("drawnow;"); #endif break; - case oAcosh: +/* case oAcosh: Stack.push(acosh(v1)); #ifdef DEBUGC mexPrintf("acosh\n"); @@ -523,7 +525,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ mexPrintf("atanh\n"); mexEvalString("drawnow;"); #endif - break; + break;*/ case oSqrt: Stack.push(sqrt(v1)); #ifdef DEBUGC @@ -550,8 +552,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ case FENDBLOCK : //it's the block end #ifdef DEBUGC - mexPrintf("FENDBLOCK\n"); - mexEvalString("drawnow;"); + //mexPrintf("FENDBLOCK\n"); + //mexEvalString("drawnow;"); #endif //Block[Block_Count].end=get_code_pos; go_on=false; @@ -569,7 +571,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ mexPrintf("FOK\n"); mexEvalString("drawnow;"); #endif - op=get_code_int + op=get_code_int; #ifdef DEBUGC mexPrintf("var=%d\n",op); #endif @@ -626,10 +628,12 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba { case EVALUATE_FORWARD : //case EVALUATE_FORWARD_R : -#ifdef DEBUGC +//#ifdef DEBUGC mexPrintf("EVALUATE_FORWARD\n"); -#endif + mexEvalString("drawnow;"); +//#endif begining=get_code_pointer; + //mexPrintf("y_kmin=%d periods=%d",y_kmin, periods); for (it_=y_kmin;it_=y_kmin;it_--) { @@ -670,9 +662,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba } break; case SOLVE_FORWARD_SIMPLE : -#ifdef DEBUGC +//#ifdef DEBUGC mexPrintf("SOLVE_FORWARD_SIMPLE\n"); -#endif + mexEvalString("drawnow;"); +//#endif g1=(double*)mxMalloc(size*size*sizeof(double)); r=(double*)mxMalloc(size*sizeof(double)); begining=get_code_pointer; @@ -707,9 +700,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba mxFree(r); break; case SOLVE_BACKWARD_SIMPLE : -#ifdef DEBUGC +//#ifdef DEBUGC mexPrintf("SOLVE_BACKWARD_SIMPLE\n"); -#endif + mexEvalString("drawnow;"); +//#endif g1=(double*)mxMalloc(size*size*sizeof(double)); r=(double*)mxMalloc(size*sizeof(double)); begining=get_code_pointer; @@ -744,85 +738,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba mxFree(g1); mxFree(r); break; - /*case SOLVE_TWO_BOUNDARIES_SIMPLE : -#ifdef DEBUGC - mexPrintf("SOLVE_TWO_BOUNDARIES_SIMPLE\n"); -#endif - is_linear=get_code_bool; - max_lag_plus_max_lead_plus_1=get_code_int; - symbol_table_endo_nbr=get_code_int; - Block_List_Max_Lag=get_code_int; - Block_List_Max_Lead=get_code_int; - Read_file(file_name, periods, max_lag_plus_max_lead_plus_1, 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); - g1=(double*)mxMalloc(size*size*sizeof(double)); - r=(double*)mxMalloc(size*sizeof(double)); - begining=get_code_pointer; - if (!is_linear) - { - cvg=false; - iter=0; - while (!(cvg||(iter>maxit_))) - { - res2=0; - res1=0; - max_res=0; - for (it_=y_kmin;it_ in SOLVE_FORWARD_COMPLETE : OK\n"); + /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK it_%d\n",it_); 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); @@ -911,7 +831,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba //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); cvg=false; + /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : OK it_%d\n",it_); + mexEvalString("drawnow;");*/ simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter); + /*mexPrintf("End of simulate_NG\n"); + mexEvalString("drawnow;");*/ //simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter); } /*mexPrintf("solve forward complete simulation\n"); @@ -923,15 +847,19 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba mexPrintf("\n"); }*/ } + /*mexPrintf("end of forward solve\n"); + mexEvalString("drawnow;");*/ + mxFree(index_equa); memset(direction,0,size_of_direction); mxFree(g1); mxFree(r); mxFree(u); break; case SOLVE_BACKWARD_COMPLETE : -#ifdef DEBUGC +//#ifdef DEBUGC mexPrintf("SOLVE_BACKWARD_COMPLETE\n"); -#endif + mexEvalString("drawnow;"); +//#endif is_linear=get_code_bool; max_lag_plus_max_lead_plus_1=get_code_int; symbol_table_endo_nbr=get_code_int; @@ -939,7 +867,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba Block_List_Max_Lead=get_code_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); - u_count_int=get_code_int + 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); @@ -1011,36 +939,36 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba #if GNUVER >= 432 //mexPrintf("omp_get_max_threads=%d\n",omp_get_max_threads()); #endif -#ifdef DEBUGC +//#ifdef DEBUGC mexPrintf("SOLVE_TWO_BOUNDARIES_COMPLETE\n"); mexEvalString("drawnow;"); -#endif +//#endif is_linear=get_code_bool; #ifdef DEBUGC mexPrintf("is_linear=%d\n",is_linear); mexEvalString("drawnow;"); #endif - max_lag_plus_max_lead_plus_1=get_code_int + max_lag_plus_max_lead_plus_1=get_code_int; #ifdef DEBUGC mexPrintf("max_lag_plus_max_lead_plus_1=%d\n",max_lag_plus_max_lead_plus_1); mexEvalString("drawnow;"); #endif - symbol_table_endo_nbr=get_code_int + symbol_table_endo_nbr=get_code_int; #ifdef DEBUGC mexPrintf("symbol_table_endo_nbr=%d\n",symbol_table_endo_nbr); mexEvalString("drawnow;"); #endif - Block_List_Max_Lag=get_code_int + Block_List_Max_Lag=get_code_int; #ifdef DEBUGC mexPrintf("Block_List_Max_Lag=%d\n",Block_List_Max_Lag); mexEvalString("drawnow;"); #endif - Block_List_Max_Lead=get_code_int + Block_List_Max_Lead=get_code_int; #ifdef DEBUGC mexPrintf("Block_List_Max_Lead=%d\n",Block_List_Max_Lead); mexEvalString("drawnow;"); #endif - u_count_int=get_code_int + u_count_int=get_code_int; #ifdef DEBUGC mexPrintf("u_count_int=%d\n",u_count_int); mexPrintf("periods=%d\n",periods); @@ -1113,7 +1041,6 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba if (isnan(res1)||isinf(res1)) { memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin)); - //GaussSeidel=false; break; } for (i=0; i< size; i++) @@ -1153,6 +1080,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba /*mexPrintf("u_count=%d &u_count=%x\n",u_count,&u_count); mexEvalString("drawnow;");*/ simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter); + /*mexPrintf("after simulate_NG1\n"); + mexEvalString("drawnow;");*/ /*mexPrintf("aft simulate_NG1\n"); mexEvalString("drawnow;");*/ } @@ -1195,7 +1124,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba res1=res2=max_res=0;max_res_idx=0; cvg = false; if(Gaussian_Elimination) - simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter); + { + simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter); + /*mexPrintf("after simulate_NG1\n"); + mexEvalString("drawnow;");*/ + } else { #ifdef LINBCG @@ -1217,11 +1150,16 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba #endif //mxFree(g1); + /*mexPrintf("end of simulate_a_block\n"); + mexEvalString("drawnow;");*/ mxFree(r); mxFree(y_save); mxFree(u); mxFree(index_vara); + mxFree(index_equa); memset(direction,0,size_of_direction); + /*mexPrintf("after free\n"); + mexEvalString("drawnow;");*/ //GaussSeidel=false; break; default: @@ -1230,6 +1168,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba mexEvalString("drawnow;"); mexErrMsgTxt("End of simulate"); } + /*mexPrintf("finish simulate_a_block\n"); + mexEvalString("drawnow;");*/ } void @@ -1237,6 +1177,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename) { ifstream CompiledCode; int Code_Size, var; + //printf("open(%s)\n",(file_name + ".cod").c_str()); //First read and store inn memory the code CompiledCode.open((file_name + ".cod").c_str(),std::ios::in | std::ios::binary| std::ios::ate); if (!CompiledCode.is_open()) @@ -1279,77 +1220,89 @@ Interpreter::compute_blocks(string file_name, string bin_basename) { case FBEGINBLOCK : //it's a new block + Block_Count++; Block_type lBlock; Block.clear(); Block_Contain.clear(); Block_contain_type lBlock_Contain; -#ifdef DEBUGC +//#ifdef DEBUGC mexPrintf("FBEGINBLOCK\n"); mexEvalString("drawnow;"); -#endif - lBlock.begin=get_code_pos-(long int)(Init_Code); -#ifdef DEBUGC +//#endif + /*uint64_t *Init_Code_64; + memcpy(Init_Code_64, Init_Code, sizeof(*Init_Code_64));*/ + lBlock.begin=get_code_pos-(uint64_t*)Init_Code; +//#ifdef DEBUGC mexPrintf("Block[%d].begin=%d\n",Block_Count, lBlock.begin); mexEvalString("drawnow;"); -#endif - lBlock.size=get_code_int -#ifdef DEBUGC +//#endif + lBlock.size=get_code_int; +//#ifdef DEBUGC mexPrintf("Block[Block_Count].size=%d\n",lBlock.size); mexEvalString("drawnow;"); -#endif - lBlock.type=get_code_int -#ifdef DEBUGC +//#endif + lBlock.type=get_code_int; +//#ifdef DEBUGC mexPrintf("Block[Block_Count].type=%d\n",lBlock.type); mexEvalString("drawnow;"); -#endif +//#endif Block.push_back(lBlock); for (int i=0;i. */ - +#include #include "Mem_Mngr.hh" Mem_Mngr::Mem_Mngr() @@ -82,6 +82,7 @@ Mem_Mngr::mxMalloc_NZE() mexPrintf("CHUNK_BLCK_SIZE=%d\n",CHUNK_BLCK_SIZE); #endif NZE_Mem=(NonZeroElem*)mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); + //mexPrintf("in mxMalloc NZE_Mem=%x CHUNK_heap_pos=%d CHUNK_BLCK_SIZE=%d Nb_CHUNK=%d\n",NZE_Mem, CHUNK_heap_pos, CHUNK_BLCK_SIZE, Nb_CHUNK); if(!NZE_Mem) { mexPrintf("Not enough memory available\n"); @@ -123,7 +124,7 @@ Mem_Mngr::mxFree_NZE(void* pos) mexPrintf("NZE_Mem_add[i*CHUNK_BLCK_SIZE]=%d\n",NZE_Mem_add[i*CHUNK_BLCK_SIZE]); mexEvalString("drawnow;"); }*/ - gap=((long int)(pos)-(long int)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem); + gap=((uint64_t)(pos)-(uint64_t)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem); if ((gap=0)) break; } @@ -277,8 +278,12 @@ void Mem_Mngr::Free_All() { int i; + /*mexPrintf("Nb_CHUNK=%d\n",Nb_CHUNK); + mexEvalString("drawnow;");*/ for (i=0;ic_index!=c /*&& first->NZE_R_N*/) + while (first->c_index!=c) { -#ifdef PRINT_OUT - mexPrintf("looking not CRS [%d, %d]\n",first->r_index,first->c_index); - mexEvalString("drawnow;"); -#endif firsta=first; first=first->NZE_R_N; } @@ -186,12 +182,8 @@ void SparseMatrix::Delete(const int r,const int c, const int Size) #endif first=FNZE_C[c]; firsta=NULL; - while (first->r_index!=r /*&& first->NZE_C_N*/) + while (first->r_index!=r) { -#ifdef PRINT_OUT - mexPrintf("looking not CSS [%d, %d]\n",first->r_index,first->c_index); - mexEvalString("drawnow;"); -#endif firsta=first; first=first->NZE_C_N; } @@ -945,6 +937,7 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i mexPrintf("at %d same construction\n",beg_t); else mexPrintf("at %d different construction\n",beg_t); + mexEvalString("drawnow;"); #endif // the same pivot for all remaining periods if (OK) @@ -959,17 +952,17 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i } if (OK) { -#ifdef WRITE_u +/*#ifdef WRITE_u long int i_toto=0; fstream toto; toto.open("compare_s.txt", std::ios::out); -#endif +#endif*/ if (max_save_ops_first>=u_count_alloc) { u_count_alloc+=5*u_count_alloc_save; -#ifdef MEM_ALLOC_CHK +/*#ifdef MEM_ALLOC_CHK mexPrintf("u=(double*)mxRealloc(u,u_count_alloc*sizeof(double))=%d, t=%d, omp_get_thread_num()=%d\n",u_count_alloc,t,omp_get_thread_num()); -#endif +#endif*/ u=(double*)mxRealloc(u,u_count_alloc*sizeof(double)); if (!u) { @@ -978,48 +971,29 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i mexErrMsgTxt("Exit from Dynare"); } } - /*#pragma omp parallel num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) - mexPrintf("omp_in_parallel=%hd, omp_get_thread_num=%d, omp_get_num_threads=%d, t=%d\n",omp_in_parallel(), omp_get_thread_num(), omp_get_num_threads(), t);*/ - //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(i,j, save_op_s, index_d, r) schedule(dynamic) for (t=1;tfirst+t*diff1[j]; - /*if(i==0) - mexPrintf("t=%d\n",t);*/ switch (save_op_s->operat) { case IFLD : r=u[index_d]; -#ifdef PRINT_u - mexPrintf("FLD u[%d] (%f)\n",index_d,u[index_d]); -#endif i+=2; break; case IFDIV : u[index_d]/=r; -#ifdef PRINT_u - mexPrintf("FDIV u[%d](%f)/=r(%f)=(%f)\n",index_d,u[index_d],r,u[index_d]); -#endif i+=2; break; case IFSUB : u[index_d]-=u[save_op_s->second+t*diff2[j]]*r; -#ifdef PRINT_u - mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,save_op_s->second+t*diff2[j],u[save_op_s->second+t*diff2[j]],r,u[index_d] ); -#endif i+=3; break; case IFLESS: u[index_d]=-u[save_op_s->second+t*diff2[j]]*r; -#ifdef PRINT_u - mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,save_op_s->second+t*diff2[j],u[save_op_s->second+t*diff2[j]],r,u[index_d] ); -#endif i+=3; break; } @@ -1027,13 +1001,7 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i if (index_d+3>=u_count_alloc) { u_count_alloc+=2*u_count_alloc_save; -#ifdef MEM_ALLOC_CHK - mexPrintf("u=(double*)mxRealloc(u,u_count_alloc*sizeof(double))=%d, t=%d, omp_get_thread_num()=%d\n",u_count_alloc,t,omp_get_thread_num()); -#endif u=(double*)mxRealloc(u,u_count_alloc*sizeof(double)); -#ifdef MEM_ALLOC_CHK - mexPrintf("ok\n"); -#endif if (!u) { mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n",u_count_alloc*sizeof(double)); @@ -1044,6 +1012,10 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i } } int t1=periods-beg_t-max(y_kmax,y_kmin); +#ifdef PROFILER + mexPrintf("first step done\n"); + mexEvalString("drawnow;"); +#endif //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(t, i,j, save_op_s, index_d, r) schedule(dynamic) for (t=t1;t0) - period = y_decal * y_size; - //pctimer_t t1 = pctimer(); - clock_t t1 = clock(); - double uu, yy; - //char tmp_s[150]; - //mexPrintf("period=%d\n",period); - //mexPrintf("Direct_Simulate\n"); -#ifdef PRINT_OUT - for (j = 0;j < it_ -y_kmin;j++) - { - for (i = 0;i < u_size;i++) - { - mexPrintf("u[%d]=%f ", j*u_size + i, double(u[j*u_size + i])); - } - mexPrintf("\n"); - } -#endif - if (/*nb_first_table_u*/print_it > 0) - { - first_count_loop = it_; - //mexPrintf("nb_prologue_table_y=%d Size%d\n",nb_prologue_table_y,Size); - /*s_middle_count_loop = it_ -y_kmin - middle_count_loop + 1;*/ - s_middle_count_loop = it_ - y_kmin - (nb_prologue_table_y / Size); -#ifdef PRINT_OUT - mexPrintf("nb_first_table_u=%d first_count_loop=%d s_middle_count_loop=%d\n",nb_first_table_u,first_count_loop, s_middle_count_loop); - mexPrintf("y_kmin=%d y_kmax=%d \n",y_kmin,y_kmax); - mexPrintf("middle_count_loop=%d\n",middle_count_loop); - mexPrintf("it_=%d\n",it_); -#endif -//#ifdef PRINT_OUT - mexPrintf("----------------------------------------------------------------------\n"); - mexPrintf(" Simulate iteration° %d \n",iter+1); - mexPrintf(" max. error=%.10e \n",double(max_res)); - mexPrintf(" sqr. error=%.10e \n",double(res2)); - mexPrintf(" abs. error=%.10e \n",double(res1)); - mexPrintf("----------------------------------------------------------------------\n"); -//endif -#ifdef PRINT_OUT - mexPrintf("first_count\n"); - mexPrintf("(first_count_loop=%d - y_kmin=%d)=%d\n",first_count_loop,y_kmin, first_count_loop-y_kmin); -#endif - } - nop = 0; - - for (j = 0 ;j < first_count_loop - y_kmin;j++) - { -#ifdef PRINT_OUT_b - mexPrintf("------------------------------------------------------\n"); - mexPrintf("j=%d\n",j); -#endif - first_table_u = F_first_table_u->pNext; - first_i_table_u = F_first_i_table_u->pNext; - for (i = 0;i < nb_first_table_u;i++) - { - switch (first_table_u->type) - { - case 1: - u[first_table_u->index + j*first_i_table_u->index] = u[first_table_u->op1 + j * first_i_table_u->op1] * u[first_table_u->op2 + j * first_i_table_u->op2]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=u[%d]*u[%d]=%f\n", first_table_u->index + j*first_i_table_u->index , first_table_u->op1 + j*first_i_table_u->op1, first_table_u->op2 + j*first_i_table_u->op2, double(u[first_table_u->index + j*first_i_table_u->index])); -#endif - break; - case 2: - u[first_table_u->index + j*first_i_table_u->index] += u[first_table_u->op1 + j * first_i_table_u->op1] * u[first_table_u->op2 + j * first_i_table_u->op2]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]+=u[%d]*u[%d]=%f\n" , first_table_u->index + j*first_i_table_u->index, first_table_u->op1 + j*first_i_table_u->op1, first_table_u->op2 + j*first_i_table_u->op2, double(u[first_table_u->index + j*first_i_table_u->index])); -#endif - break; - case 3: - u[first_table_u->index + j*first_i_table_u->index] = 1 / ( -u[first_table_u->op1 + j * first_i_table_u->op1]); -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=1/(-u[%d])=%f\n", first_table_u->index + j*first_i_table_u->index, first_table_u->op1 + j*first_i_table_u->op1, double(u[first_table_u->index + j*first_i_table_u->index])); -#endif - break; - case 5: - Stack.push(u[first_table_u->index + j*first_i_table_u->index]); -#ifdef PRINT_OUT_b - mexPrintf("push(u[%d])\n", first_table_u->index + j*first_i_table_u->index); -#endif - break; - case 6: - u[first_table_u->index + j*first_i_table_u->index] = 1 / (1 - u[first_table_u->op1 + j * first_i_table_u->op1] * u[first_table_u->op2 + j * first_i_table_u->op2]); -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=1/(1-u[%d]*u[%d])=%f\n", first_table_u->index + j*first_i_table_u->index, first_table_u->op1 + j*first_i_table_u->op1, first_table_u->op2 + j*first_i_table_u->op2, double(u[first_table_u->index + j*first_i_table_u->index])); -#endif - break; - case 7: - u[first_table_u->index + j*first_i_table_u->index] *= u[first_table_u->op1 + j * first_i_table_u->op1]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]*=u[%d]=%f\n", first_table_u->index + j*first_i_table_u->index, first_table_u->op1 + j*first_i_table_u->op1, double(u[first_table_u->index + j*first_i_table_u->index])); -#endif - break; - } - if (isnan(u[first_table_u->index+ j*first_i_table_u->index]) || isinf(u[first_table_u->index+ j*first_i_table_u->index])) - { - mexPrintf("Error during the computation of u[%d] at time %d (in first_table_u) (operation type %d)",first_table_u->index,j,int(first_table_u->type)); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - else if (fabs(u[first_table_u->index+ j*first_i_table_u->index])>very_big) - { - mexPrintf("(first) big u[%d]=%f>%f in type=%d",first_table_u->index+ j*first_i_table_u->index, double(u[first_table_u->index+ j*first_i_table_u->index]),very_big,first_table_u->type); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - first_table_u = first_table_u->pNext; - first_i_table_u = first_i_table_u->pNext; - nop++; - } - } -#ifdef PRINT_OUT_b - mexPrintf("//////////////////////////////////////////////////////////////\n"); - mexPrintf("prologue\n"); -#endif - //int nb_prologue_push=0; - prologue_table_u = F_prologue_table_u->pNext; - for (i = 0;i < nb_prologue_table_u ;i++) - { - switch (prologue_table_u->type) - { - case 1: - u[prologue_table_u->index ] = u[prologue_table_u->op1 ] * u[prologue_table_u->op2 ]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=u[%d]*u[%d]=%f\n", prologue_table_u->index , prologue_table_u->op1 , prologue_table_u->op2 , double(u[prologue_table_u->index ])); -#endif - break; - case 2: - u[prologue_table_u->index ] += u[prologue_table_u->op1 ] * u[prologue_table_u->op2 ]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]+=u[%d]*u[%d]=%f\n" , prologue_table_u->index , prologue_table_u->op1 , prologue_table_u->op2 , double(u[prologue_table_u->index ])); -#endif - break; - case 3: - u[prologue_table_u->index ] = 1 / ( -u[prologue_table_u->op1 ]); -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=1/(-u[%d])=%f\n", prologue_table_u->index, prologue_table_u->op1, double(u[prologue_table_u->index])); -#endif - break; - case 5: - //nb_prologue_push++; - Stack.push(u[prologue_table_u->index]); -#ifdef PRINT_OUT_b - mexPrintf("push(u[%d])\n", prologue_table_u->index ); -#endif - break; - case 6: - u[prologue_table_u->index ] = 1 / (1 - u[prologue_table_u->op1] * u[prologue_table_u->op2]); -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=1/(1-u[%d]*u[%d])=%f\n", prologue_table_u->index, prologue_table_u->op1, prologue_table_u->op2, double(u[prologue_table_u->index])); -#endif - break; - case 7: - u[prologue_table_u->index] *= u[prologue_table_u->op1]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]*=u[%d]=%f\n", prologue_table_u->index, prologue_table_u->op1, double(u[prologue_table_u->index])); -#endif - break; - } - if (isnan(u[prologue_table_u->index]) || isinf(u[prologue_table_u->index])) - { - mexPrintf("Error during the computation of u[%d] type=%d (in prologue_table_u)",prologue_table_u->index, prologue_table_u->type); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - else if (fabs(u[prologue_table_u->index])>very_big) - { - mexPrintf("(prologue) big u[%d]=%f>%f type=%d",prologue_table_u->index, double(u[prologue_table_u->index]), very_big, prologue_table_u->type); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - prologue_table_u = prologue_table_u->pNext; - nop++; - } -#ifdef PRINT_OUT_b - mexPrintf("middle_u (s_middle_count_loop=%d\n", s_middle_count_loop); -#endif - //int nb_middle_push=0; - for (j = 0;j < s_middle_count_loop /*- y_kmin*/;j++) - { - //cout << "j=" << j << "\n"; -#ifdef PRINT_OUT_b - mexPrintf("-----------------------------------------------------------------\n"); -#endif - middle_table_u = F_middle_table_u->pNext; - middle_i_table_u = F_middle_i_table_u->pNext; - for (i = 0;i < nb_middle_table_u;i++) - { - switch (middle_table_u->type) - { - case 1: - u[middle_table_u->index + j*middle_i_table_u->index] = u[middle_table_u->op1 + j * middle_i_table_u->op1] * u[middle_table_u->op2 + j * middle_i_table_u->op2]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d+%d*%d=%d]=u[%d]*u[%d]=%f\n", middle_table_u->index, j, middle_i_table_u->index, middle_table_u->index + j*middle_i_table_u->index, middle_table_u->op1 + j*middle_i_table_u->op1, middle_table_u->op2 + j*middle_i_table_u->op2, double(u[middle_table_u->index + j*middle_i_table_u->index])); -#endif - break; - case 2: - u[middle_table_u->index + j*middle_i_table_u->index] += u[middle_table_u->op1 + j * middle_i_table_u->op1] * u[middle_table_u->op2 + j * middle_i_table_u->op2]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d+%d*%d=%d]+=u[%d]*u[%d]=%f\n" , middle_table_u->index, j, middle_i_table_u->index , middle_table_u->index + j*middle_i_table_u->index , middle_table_u->op1 + j*middle_i_table_u->op1, middle_table_u->op2 + j*middle_i_table_u->op2, double(u[middle_table_u->index + j*middle_i_table_u->index])); -#endif - break; - case 3: - u[middle_table_u->index + middle_i_table_u->index] = 1 / ( -u[middle_table_u->op1 + j * middle_i_table_u->op1]); -#ifdef PRINT_OUT_b - mexPrintf("u[%d+%d*%d=%d]=1/(-u[%d])=%f\n", middle_table_u->index, j, middle_i_table_u->index, middle_table_u->index + j*middle_i_table_u->index, middle_table_u->op1 + j*middle_i_table_u->op1, double(u[middle_table_u->index + j*middle_i_table_u->index])); -#endif - break; - case 5: -#ifdef PRINT_OUT_b - mexPrintf("push(u[%d+%d*%d=%d])\n", middle_table_u->index, j, middle_i_table_u->index, middle_table_u->index + j*middle_i_table_u->index); -#endif - //nb_middle_push++; - Stack.push(u[middle_table_u->index + j*middle_i_table_u->index]); - break; - case 6: - u[middle_table_u->index + j*middle_i_table_u->index] = 1 / (1 - u[middle_table_u->op1 + j * middle_i_table_u->op1] * u[middle_table_u->op2 + j * middle_i_table_u->op2]); -#ifdef PRINT_OUT_b - mexPrintf("u[%d+%d*%d=%d]=1/(1-u[%d]*u[%d])=%f\n", middle_table_u->index, j, middle_i_table_u->index, middle_table_u->index + j*middle_i_table_u->index, middle_table_u->op1 + j*middle_i_table_u->op1, middle_table_u->op2 + j*middle_i_table_u->op2, double(u[middle_table_u->index + j*middle_i_table_u->index])); -#endif - break; - case 7: - u[middle_table_u->index + j*middle_i_table_u->index] *= u[middle_table_u->op1 + j * middle_i_table_u->op1]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d+%d*%d=%d]*=u[%d]=%f\n", middle_table_u->index, j, middle_i_table_u->index, middle_table_u->index + j*middle_i_table_u->index, middle_table_u->op1 + j*middle_i_table_u->op1, double(u[middle_table_u->index + j*middle_i_table_u->index])); -#endif - break; - } - if (isnan(u[middle_table_u->index+ j*middle_i_table_u->index]) || isinf(u[middle_table_u->index+ j*middle_i_table_u->index])) - { - mexPrintf("Error during the computation of u[%d] at time %d type=%d (in middle_table_u)",middle_table_u->index,j,middle_table_u->type); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - else if (fabs(u[middle_table_u->index+ j*middle_i_table_u->index])>very_big) - { - mexPrintf("(middle) big u[%d]=%f>%f type=%d",middle_table_u->index+ j*middle_i_table_u->index, double(u[middle_table_u->index+ j*middle_i_table_u->index]), very_big,middle_table_u->type); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - middle_table_u = middle_table_u->pNext; - middle_i_table_u = middle_i_table_u->pNext; - nop++; - } - } -#ifdef PRINT_OUT - mexPrintf("last_u\n"); -#endif - last_table_u = F_last_table_u->pNext; - for (i = 0;i < nb_last_table_u ;i++) - { - switch (last_table_u->type) - { - case 1: - u[last_table_u->index] = u[last_table_u->op1] * u[last_table_u->op2]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=u[%d]*u[%d]=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, double(u[last_table_u->index])); -#endif - break; - case 2: - u[last_table_u->index] += u[last_table_u->op1] * u[last_table_u->op2]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]+=u[%d]*u[%d]=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, double(u[last_table_u->index])); -#endif - break; - case 3: - u[last_table_u->index] = 1 / ( -u[last_table_u->op1]); -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=1/(-u[%d])=%f\n", last_table_u->index, last_table_u->op1, double(u[last_table_u->index])); -#endif - break; - case 5: - Stack.push(u[last_table_u->index]); -#ifdef PRINT_OUT_b - mexPrintf("push(u[%d])\n", last_table_u->index); -#endif - break; - case 6: - u[last_table_u->index] = 1 / (1 - u[last_table_u->op1] * u[last_table_u->op2]); -#ifdef PRINT_OUT_b - mexPrintf("u[%d]=1/(1-u[%d]*u[%d])=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, double(u[last_table_u->index])); -#endif - break; - case 7: - u[last_table_u->index] *= u[last_table_u->op1]; -#ifdef PRINT_OUT_b - mexPrintf("u[%d]*=u[%d]=%f\n", last_table_u->index, last_table_u->op1, double(u[last_table_u->index])); -#endif - break; - } - if (isnan(u[last_table_u->index]) || isinf(u[last_table_u->index])) - { - mexPrintf("Error during the computation of u[%d] type=%d (in last_table_u)",last_table_u->index,last_table_u->type); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - else if (fabs(u[last_table_u->index])>very_big) - { - mexPrintf("(last) big u[%d]=%f>%f type=%d",last_table_u->index, double(u[last_table_u->index]),very_big,last_table_u->type); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - last_table_u = last_table_u->pNext; - nop++; - } - res1 = res2 =max_res = 0; - for (i = nb_last_table_y - 1;i >= 0;i--) - { - k = last_table_y[i].index; - yy = 0; - /*y[period + k] = 0;*/ -#ifdef PRINT_OUT_y - //mexPrintf("it_=%d\n", it_); - mexPrintf("->y[it_*y_size+%d]=y[%d]=", k, period + k); -#endif - for (j = last_table_y[i].nb - 1;j >= 0;j--) - { - uu = Stack.top(); - Stack.pop(); - m = last_table_y[i].y_index[j]; -#ifdef PRINT_OUT_y - if (j > 0) - { - if (m >= 0) - mexPrintf("u[%d](%f)*y[%d](%f)+", last_table_y[i].u_index[j], double(uu), m + period, double(y[period + m])); - else - mexPrintf("u[%d](%f)+", last_table_y[i].u_index[j], double(uu)); - } - else - { - if (m >= 0) - mexPrintf("u[%d](%f)*y[%d](%f)", last_table_y[i].u_index[j], double(uu), m + period, double(y[period + m])); - else - mexPrintf("u[%d](%f)", last_table_y[i].u_index[j], double(uu)); - } -#endif - if (m >= 0) - yy/*y[period + k]*/ += uu * y[period + m]; - else - yy/*y[period + k]*/ += uu; - } - if (isnan(yy) || isinf(yy)) - { - mexPrintf("Error during the computation of y[%d] at time %d (in last_table_u)",k,period); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - /*if(((k-73) % y_size)==0) - mexPrintf("y[it_*y_size +73]=%f \n",yy);*/ - err = fabs(yy - y[period + k]); - res1 += err; - if (max_resy_decal;j--) - { -#ifdef PRINT_OUT_y - mexPrintf("(per)j=%d y_decal=%d deca=%d in y compute\n",j,y_decal,deca); -#endif - for (i = nb_middle_table_y - 1;i >= 0;i--) - { - k = middle_table_y[i].index + (j-1) * middle_i_table_y[i].index; - yy = 0; -#ifdef PRINT_OUT_y - mexPrintf("(0)y[%d]=", k ); -#endif - for (l = middle_table_y[i].nb - 1;l >= 0;l--) - { - uu = Stack.top(); - Stack.pop(); - m = middle_table_y[i].y_index[l] + (j/*-deca*/-1) * middle_i_table_y[i].y_index[l]; -#ifdef PRINT_OUT_y - if (/*m*/middle_table_y[i].y_index[l] >= 0) - { - m1 = middle_table_y[i].u_index[l] + (j-deca-1) * middle_i_table_y[i].u_index[l]; - if (l > 0) - mexPrintf("u[%d](%f)*y[%d](%f)+", m1, double(uu), m, double(y[m])); - else - mexPrintf("u[%d](%f)*y[%d](%f)", m1, double(uu), m, double(y[m])); - } - else - { - m1 = middle_table_y[i].u_index[l] + (j-deca-1) * middle_i_table_y[i].u_index[l]; - if (l > 0) - mexPrintf("u[%d](%f)*y[%d](%f)+", m1, double(uu), m, 1.0); - else - mexPrintf("u[%d](%f)*y[%d](%f)", m1, double(uu), m, 1.0); - } -#endif - if (/*m*/middle_table_y[i].y_index[l] >= 0) - yy += uu * y[m]; - else - yy += uu; - } - //mexPrintf("y[%d]=%f\n",k,yy); - if (isnan(yy) || isinf(yy)) - { - mexPrintf("Error during the computation of y[%d] at time %d (in middle_table_u)",middle_table_y[i].index % y_size,j+middle_table_y[i].index / y_size); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - /*if(((k-73) % y_size)==0) - mexPrintf("y[it_*y_size +73]=%f \n",yy);*/ - err = fabs(yy - y[k]); - res1 += err; - if (max_res= 0;i--) - { - k = prologue_table_y[i].index; - yy = 0; -#ifdef PRINT_OUT_y - mexPrintf("(1)y[%d]=", k+y_decal*y_size); -#endif - for (j = prologue_table_y[i].nb - 1;j >= 0;j--) - { - //nb_prologue_pop++; - uu = Stack.top(); - Stack.pop(); -#ifdef PRINT_OUT_y - if (prologue_table_y[i].y_index[j] >= 0) - { - if (j > 0) - mexPrintf("u[%d](%f)*y[%d](%f)+", prologue_table_y[i].u_index[j], double(uu), prologue_table_y[i].y_index[j], double(y[prologue_table_y[i].y_index[j]])); - else - mexPrintf("u[%d](%f)*y[%d](%f)", prologue_table_y[i].u_index[j], double(uu), prologue_table_y[i].y_index[j], double(y[prologue_table_y[i].y_index[j]])); - } - else - { - if (j > 0) - mexPrintf("u[%d](%f)*y[%d](%f)+", prologue_table_y[i].u_index[j], double(uu), prologue_table_y[i].y_index[j], 1.0); - else - mexPrintf("u[%d](%f)*y[%d](%f)", prologue_table_y[i].u_index[j], double(uu), prologue_table_y[i].y_index[j], 1.0); - } -#endif - if (prologue_table_y[i].y_index[j] >= 0) - yy += uu * y[prologue_table_y[i].y_index[j]+y_decal*y_size]; - else - yy += uu; - } - if (isnan(yy) || isinf(yy)) - { - mexPrintf("Error during the computation of y[%d] at time %d (in prologue_table_u)",k % y_size, k / y_size); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - err = fabs(yy - y[k+y_decal*y_size]); - res1 += err; - if (max_res= 0;i--) - { - k = first_table_y[i].index; - yy = 0; -#ifdef PRINT_OUT_y - mexPrintf("(2)y[%d]=", k); -#endif - for (j = first_table_y[i].nb - 1;j >= 0;j--) - { - uu = Stack.top(); - Stack.pop(); -#ifdef PRINT_OUT_y - if (j > 0) - mexPrintf("u[%d](%f)*y[%d](%f)+", first_table_y[i].u_index[j], double(uu), first_table_y[i].y_index[j], double(y[first_table_y[i].y_index[j]])); - else - mexPrintf("u[%d](%f)*y[%d](%f)", first_table_y[i].u_index[j], double(uu), first_table_y[i].y_index[j], double(y[first_table_y[i].y_index[j]])); -#endif - if (m >= 0) - yy += uu * y[first_table_y[i].y_index[j]]; - else - yy += uu; - } - if (isnan(yy) || isinf(yy)) - { - mexPrintf("Error during the computation of y[%d] (in first_table_u)",k); - mexEvalString("st=fclose('all');clear all;"); - filename+=" stopped"; - mexErrMsgTxt(filename.c_str()); - } - err = fabs(yy - y[k]); - res1 += err; - if (max_res 0*/print_it) - { - mexPrintf("res1 = %.10e\n", double(res1)); - mexPrintf("res2 = %.10e\n", double(res2)); - mexPrintf("max_res = %.10e\n", double(max_res)); - mexPrintf("(**%f milliseconds u_count : %d nop : %d **)\n", 1000.0*(double(t2) - double(t1))/double(CLOCKS_PER_SEC), u_count, nop); - mexEvalString("drawnow;"); - } -} - @@ -2128,15 +1550,21 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, bool one; #endif #ifdef PROFILER - long int ndiv=0, nsub=0, ncomp=0, nmul=0; - double tinsert=0, tdelete=0, tpivot=0, tbigloop=0; - clock_t td1; - int nbpivot=0, nbdiv=0, nbless=0, nbpivot_it=0, nbRealloc=0, insert=0; + //long int ndiv=0, nsub=0, ncomp=0, nmul=0; + //double tinsert=0, tdelete=0, tpivot=0, tbigloop=0; + //clock_t td1; + int nbpivot=0, nbpivot_it=0; + //int nbdiv=0, nbless=0, nbRealloc=0, insert=0; #endif - + /*mexPrintf("begining\n"); + mexEvalString("drawnow;");*/ if (cvg) return(0); + /*mexPrintf("begining after cvg Size=%d\n", Size); + mexEvalString("drawnow;");*/ Simple_Init(it_, y_kmin, y_kmax, Size, IM_i); + /*mexPrintf("begining after Simple_Init\n"); + mexEvalString("drawnow;");*/ if (isnan(res1) || isinf(res1)) { if (slowc_save<1e-8) @@ -2155,23 +1583,34 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, iter--; return(0); } + //mexPrintf("begining after nan\n"); + //mexEvalString("drawnow;"); //mexPrintf("before the main loop y_size=%d\n",y_size); for (i=0;ir_index=%d \n",first->r_index); + //mexPrintf("line_done[%d]=%d \n",first->r_index,line_done[first->r_index]); + //mexEvalString("drawnow;"); if (!line_done[first->r_index]) { + //mexPrintf("first->u_index=%d \n",first->u_index); + //mexEvalString("drawnow;"); k=first->u_index; int jj=first->r_index; int NRow_jj=NRow(jj); @@ -2182,6 +1621,8 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, #ifdef MARKOVITZ piv_v[l]=u[k]; + //mexPrintf("piv_v[%d]=%f\n",l, piv_v[l]); + //mexEvalString("drawnow;"); double piv_fabs=fabs(u[k]); pivj_v[l]=jj; pivk_v[l]=k; @@ -2224,6 +1665,8 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, } first=first->NZE_C_N; } + /*mexPrintf("pivot found piv_v[0]=%f pivk_v[0]=%d l=%d one=%d\n", piv_v[0], pivk_v[0], l, one); + mexEvalString("drawnow;");*/ #ifdef MARKOVITZ double markovitz=0, markovitz_max=-9e70; if (!one) @@ -2245,12 +1688,14 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, 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; + //mexPrintf("stored\n"); } } } @@ -2269,6 +1714,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, filename+=" stopped"; mexErrMsgTxt(filename.c_str()); } + //mexPrintf("piv=%f\n",piv); /*divide all the non zeros elements of the line pivj by the max_pivot*/ int nb_var=At_Row(pivj,&first); for (j=0;jr_index; + //mexPrintf("j=%d row=%d line_done[row]=%d\n",j, row, line_done[row]); 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); + //mexPrintf("nb_var_sub=%d nb_var_piv=%d\n",nb_var_sub,nb_var_piv); 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; @@ -2360,8 +1808,11 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, } else first=first->NZE_C_N; + /*first=first->NZE_R_N;*/ } } + /*mexPrintf("before bcksub\n"); + mexEvalString("drawnow;");*/ double slowc_lbx=slowc, res1bx; //mexPrintf("before bksub it_=%d\n",it_); for (i=0;i0) { mexPrintf("Sim : %f ms\n",(1000.0*(double(clock())-double(time00)))/double(CLOCKS_PER_SEC)); @@ -2428,12 +1880,15 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax { mexPrintf("slowc_save=%g\n", slowc_save); for(j=0;j0)) { u_count=save_u_count; - End(Size); + } #endif + End(Size); if (print_it) { //pctimer_t t2 = pctimer(); @@ -3505,399 +2970,3 @@ SparseMatrix::fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus u_count_init=max_lag_plus_max_lead_plus_1; } - -void -SparseMatrix::read_file_table_u(t_table_u **table_u, t_table_u **F_table_u, t_table_u **i_table_u, t_table_u **F_i_table_u, int *nb_table_u, bool i_to_do, bool shifting, int *nb_add_u_count, int y_kmin, int y_kmax, int u_size) -{ - char type; - int i; - i=SaveCode.tellp(); -#ifdef PRINT_OUT - mexPrintf("SaveCode.tellp()=%d\n",i); -#endif - SaveCode.read(reinterpret_cast(nb_table_u), sizeof(*nb_table_u)); -#ifdef PRINT_OUT - mexPrintf("->*nb_table_u=%d\n", *nb_table_u); -#endif - *table_u = (t_table_u*)mxMalloc(sizeof(t_table_u) - 2 * sizeof(int)); - *F_table_u = *table_u; - if (i_to_do) - { -#ifdef PRINT_OUT - mexPrintf("=>i_table\n"); -#endif - (*i_table_u) = (t_table_u*)mxMalloc(sizeof(t_table_u) - 2 * sizeof(int)); - (*F_i_table_u) = (*i_table_u); - } - for (i = 0;i < *nb_table_u;i++) - { - SaveCode.read(reinterpret_cast(&type), sizeof(type)); - switch (type) - { - case 3: - case 7: - (*table_u)->pNext = (t_table_u*)mxMalloc(sizeof(t_table_u) - sizeof(int)); - (*table_u) = (*table_u)->pNext; - (*table_u)->type = type; - SaveCode.read(reinterpret_cast(&(*table_u)->index), sizeof((*table_u)->index)); - if ((*table_u)->index>max_u) - max_u=(*table_u)->index; - if ((*table_u)->indexindex; - SaveCode.read(reinterpret_cast(&(*table_u)->op1), sizeof((*table_u)->op1)); - if (shifting) - { - (*table_u)->index -= y_kmin * u_size; - if ((*table_u)->index>max_u) - max_u=(*table_u)->index; - if ((*table_u)->indexindex; - (*table_u)->op1 -= y_kmin * u_size; - } -#ifdef PRINT_OUT - - if ((*table_u)->type == 3) - mexPrintf("u[%d]=-1/u[%d]\n", (*table_u)->index, (*table_u)->op1); - else - mexPrintf("u[%d]*=u[%d]\n", (*table_u)->index, (*table_u)->op1); -#endif - if (i_to_do) - { - (*i_table_u)->pNext = (t_table_u*)mxMalloc(sizeof(t_table_u) - sizeof(int)); - (*i_table_u) = (*i_table_u)->pNext; - (*i_table_u)->type = type; - SaveCode.read(reinterpret_cast(&(*i_table_u)->index), sizeof((*i_table_u)->index)); - SaveCode.read(reinterpret_cast(&(*i_table_u)->op1), sizeof((*i_table_u)->op1)); -#ifdef FIXE - (*i_table_u)->index = u_size; - (*i_table_u)->op1 = u_size; -#endif -#ifdef PRINT_OUT - if ((*i_table_u)->type == 3) - mexPrintf("i u[%d]=1/(1-u[%d])\n", (*i_table_u)->index, (*i_table_u)->op1); - else - mexPrintf("i u[%d]*=u[%d]\n", (*i_table_u)->index, (*i_table_u)->op1); -#endif - } - break; - case 1: - case 2: - case 6: - (*table_u)->pNext = (t_table_u*)mxMalloc(sizeof(t_table_u)); - (*table_u) = (*table_u)->pNext; - (*table_u)->type = type; - SaveCode.read(reinterpret_cast(&(*table_u)->index), sizeof((*table_u)->index)); - if ((*table_u)->index>max_u) - max_u=(*table_u)->index; - if ((*table_u)->indexindex; - SaveCode.read(reinterpret_cast(&(*table_u)->op1), sizeof((*table_u)->op1)); - SaveCode.read(reinterpret_cast(&(*table_u)->op2), sizeof((*table_u)->op2)); - if (shifting) - { - (*table_u)->index -= y_kmin * u_size; - if ((*table_u)->index>max_u) - max_u=(*table_u)->index; - if ((*table_u)->indexindex; - (*table_u)->op1 -= y_kmin * u_size; - (*table_u)->op2 -= y_kmin * u_size; - } - if ((*table_u)->type == 1) - { -#ifdef PRINT_OUT - mexPrintf("u[%d]=u[%d]*u[%d]\n", (*table_u)->index, (*table_u)->op1, (*table_u)->op2); -#endif - if (i_to_do) - (*nb_add_u_count)++; - } -#ifdef PRINT_OUT - else if ((*table_u)->type == 2) - mexPrintf("u[%d]+=u[%d]*u[%d]\n", (*table_u)->index, (*table_u)->op1, (*table_u)->op2); - else - mexPrintf("u[%d]=1/(1-u[%d]*u[%d])\n", (*table_u)->index, (*table_u)->op1, (*table_u)->op2); -#endif - if (i_to_do) - { - (*i_table_u)->pNext = (t_table_u*)mxMalloc(sizeof(t_table_u)); - (*i_table_u) = (*i_table_u)->pNext; - (*i_table_u)->type = type; - SaveCode.read(reinterpret_cast(&(*i_table_u)->index), sizeof((*i_table_u)->index)); - SaveCode.read(reinterpret_cast(&(*i_table_u)->op1), sizeof((*i_table_u)->op1)); - SaveCode.read(reinterpret_cast(&(*i_table_u)->op2), sizeof((*i_table_u)->op2)); -#ifdef FIXE - (*i_table_u)->index = u_size; - (*i_table_u)->op1 = u_size; - (*i_table_u)->op2 = u_size; -#endif -#ifdef PRINT_OUT - if ((*i_table_u)->type == 1) - mexPrintf("i u[%d]=u[%d]*u[%d]\n", (*i_table_u)->index, (*i_table_u)->op1, (*i_table_u)->op2); - else if ((*i_table_u)->type == 2) - mexPrintf("i u[%d]+=u[%d]*u[%d]\n", (*i_table_u)->index, (*i_table_u)->op1, (*i_table_u)->op2); - else - mexPrintf("i u[%d]=1/(1-u[%d]*u[%d])\n", (*i_table_u)->index, (*i_table_u)->op1, (*i_table_u)->op2); -#endif - } - break; - case 5: - (*table_u)->pNext = (t_table_u*)mxMalloc(sizeof(t_table_u) - 2 * sizeof(int)); - (*table_u) = (*table_u)->pNext; - (*table_u)->type = type; - SaveCode.read(reinterpret_cast(&(*table_u)->index), sizeof((*table_u)->index)); - if ((*table_u)->index>max_u) - max_u=(*table_u)->index; - if ((*table_u)->indexindex; - if (shifting) - { - (*table_u)->index -= y_kmin * u_size; - if ((*table_u)->index>max_u) - max_u=(*table_u)->index; - if ((*table_u)->indexindex; - } -#ifdef PRINT_OUT - mexPrintf("push(u[%d])\n", (*table_u)->index); -#endif - if (i_to_do) - { - (*i_table_u)->pNext = (t_table_u*)mxMalloc(sizeof(t_table_u) - 2 * sizeof(int)); - (*i_table_u) = (*i_table_u)->pNext; - (*i_table_u)->type = type; - SaveCode.read(reinterpret_cast(&(*i_table_u)->index), sizeof((*i_table_u)->index)); -#ifdef FIXE - (*i_table_u)->index = u_size; -#endif -#ifdef PRINT_OUT - mexPrintf("i push(u[%d])\n", (*i_table_u)->index); -#endif - } - break; - } - } -} - -void -SparseMatrix::read_file_table_y(t_table_y **table_y, t_table_y **i_table_y, int *nb_table_y, bool i_to_do, bool shifting, int y_kmin, int y_kmax, int u_size, int y_size) -{ - int i, k; - SaveCode.read(reinterpret_cast(nb_table_y), sizeof(*nb_table_y)); -#ifdef PRINT_OUT - mexPrintf("nb_table_y=%d\n", *nb_table_y); - //mexPrintf("y_size=%d, u_size=%d, y_kmin=%d, y_kmax=%d\n", y_size, u_size, y_kmin, y_kmax); -#endif - (*table_y) = (t_table_y*)mxMalloc((*nb_table_y) * sizeof(t_table_y)); - for (i = 0;i < *nb_table_y;i++) - { - SaveCode.read(reinterpret_cast(&((*table_y)[i].index)), sizeof((*table_y)[i].index)); - SaveCode.read(reinterpret_cast(&((*table_y)[i].nb)), sizeof((*table_y)[i].nb)); - if (shifting) - (*table_y)[i].index -= y_kmin * y_size; -#ifdef PRINT_OUT - mexPrintf("table_y[i].nb=%d\n", (*table_y)[i].nb); - mexPrintf("y[%d]=", (*table_y)[i].index); -#endif - (*table_y)[i].u_index = (int*)mxMalloc((*table_y)[i].nb * sizeof(int)); - (*table_y)[i].y_index = (int*)mxMalloc((*table_y)[i].nb * sizeof(int)); - for (k = 0;k < (*table_y)[i].nb;k++) - { - SaveCode.read(reinterpret_cast(&((*table_y)[i].u_index[k])), sizeof((*table_y)[i].u_index[k])); - SaveCode.read(reinterpret_cast(&((*table_y)[i].y_index[k])), sizeof((*table_y)[i].y_index[k])); - if (shifting) - { - (*table_y)[i].u_index[k] -= y_kmin * u_size; - if (((*table_y)[i].y_index[k] > y_size*y_kmin) && ((*table_y)[i].y_index[k] < y_size*(2*y_kmin + y_kmax + 2))) - { - (*table_y)[i].y_index[k] -= y_kmin * y_size; - } - } -#ifdef PRINT_OUT - if (k < (*table_y)[i].nb - 1) - mexPrintf("u[%d]*y[%d]+", (*table_y)[i].u_index[k], (*table_y)[i].y_index[k]); - else - mexPrintf("u[%d]*y[%d]\n", (*table_y)[i].u_index[k], (*table_y)[i].y_index[k]); -#endif - } - } -#ifdef PRINT_OUT - mexPrintf("*nb_table_y=%d\n", *nb_table_y); -#endif - if (i_to_do) - { - *i_table_y = (t_table_y*)mxMalloc((*nb_table_y) * sizeof(t_table_y)); - for (i = 0;i < *nb_table_y;i++) - { - SaveCode.read(reinterpret_cast(&((*i_table_y)[i].index)), sizeof((*i_table_y)[i].index)); - SaveCode.read(reinterpret_cast(&((*i_table_y)[i].nb)), sizeof((*i_table_y)[i].nb)); -#ifdef PRINT_OUT - mexPrintf("(*i_table_y)[i].nb=%d\n", (*i_table_y)[i].nb); - mexPrintf("y_i[%d]=", (*i_table_y)[i].index); -#endif - (*i_table_y)[i].u_index = (int*)mxMalloc((*i_table_y)[i].nb * sizeof(int)); - (*i_table_y)[i].y_index = (int*)mxMalloc((*i_table_y)[i].nb * sizeof(int)); - for (k = 0;k < (*i_table_y)[i].nb;k++) - { - SaveCode.read(reinterpret_cast(&((*i_table_y)[i].u_index[k])), sizeof((*i_table_y)[i].u_index[k])); - SaveCode.read(reinterpret_cast(&((*i_table_y)[i].y_index[k])), sizeof((*i_table_y)[i].y_index[k])); -#ifdef PRINT_OUT - if (k < (*i_table_y)[i].nb - 1) - mexPrintf("u[%d]*y[%d]+", (*i_table_y)[i].u_index[k], (*i_table_y)[i].y_index[k]); - else - mexPrintf("u[%d]*y[%d]\n", (*i_table_y)[i].u_index[k], (*i_table_y)[i].y_index[k]); -#endif - } - } - } -} - - - - - -void -SparseMatrix::Read_file(std::string file_name, int periods, int u_size, int y_size, int y_kmin, int y_kmax, int &nb_endo, int &u_count, int &u_count_init, double *u) -{ - int nb_add_u_count = 0; - filename=file_name; - //mexPrintf("-------------------------------------------------\n"); -#ifdef PRINT_OUT - mexPrintf("min_u(initial)=%d\n",min_u); - mexPrintf("%s\n", file_name.c_str()); -#endif - if (!SaveCode.is_open()) - { -#ifdef PRINT_OUT - mexPrintf("file opened\n"); -#endif - SaveCode.open((file_name + ".bin").c_str(), std::ios::in | std::ios::binary); - if (!SaveCode.is_open()) - { - mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + ".bin").c_str()); - mexEvalString("st=fclose('all');clear all;"); - mexErrMsgTxt("Exit from Dynare"); - } -#ifdef PRINT_OUT - mexPrintf("done\n"); -#endif - } - SaveCode.read(reinterpret_cast(&nb_endo), sizeof(nb_endo)); - SaveCode.read(reinterpret_cast(&u_count), sizeof(u_count)); - SaveCode.read(reinterpret_cast(&u_count_init), sizeof(u_count_init)); -#ifdef PRINT_OUT - mexPrintf("nb_endo=%d\n", nb_endo); - mexPrintf("u_count=%d\n", u_count); - mexPrintf("u_count_init=%d\n", u_count_init); - //mexPrintf("first table_u\n"); -#endif - read_file_table_u(&first_table_u, &F_first_table_u, &first_i_table_u, &F_first_i_table_u, &nb_first_table_u, true, false, &nb_add_u_count, y_kmin, y_kmax, u_size); -#ifdef PRINT_OUT - mexPrintf("nb_first_table_u=%d\n", nb_first_table_u); -#endif -//mexErrMsgTxt("Exit from Dynare"); -#ifdef PRINT_OUT - mexPrintf("prologue table_u\n"); -#endif - read_file_table_u(&prologue_table_u, &F_prologue_table_u, NULL, NULL, &nb_prologue_table_u, false, false, &nb_add_u_count, y_kmin, y_kmax, u_size); -#ifdef PRINT_OUT - mexPrintf("nb_prologue_table_u=%d\n", nb_prologue_table_u); -#endif - //mexErrMsgTxt("Exit from Dynare"); - SaveCode.read(reinterpret_cast(&middle_count_loop), sizeof(middle_count_loop)); -#ifdef PRINT_OUT - mexPrintf("middle_count_loop=%d\n",middle_count_loop); -#endif - //mexErrMsgTxt("Exit from Dynare"); -#ifdef PRINT_OUT - mexPrintf("middle table_u\n"); -#endif - read_file_table_u(&middle_table_u, &F_middle_table_u, &middle_i_table_u, &F_middle_i_table_u, &nb_middle_table_u, true, /*true*/false, &nb_add_u_count, y_kmin, y_kmax, u_size); -#ifdef PRINT_OUT - mexPrintf("nb_middle_table_u=%d\n",nb_middle_table_u); - //mexPrintf("last table_u\n"); -#endif - read_file_table_u(&last_table_u, &F_last_table_u, NULL, NULL, &nb_last_table_u, false, false, &nb_add_u_count, y_kmin, y_kmax, u_size); -#ifdef PRINT_OUT - mexPrintf("->nb_last_table_u=%d\n", nb_last_table_u); - mexPrintf("i=%d\n", i); - mexPrintf("going to read prologue_table_y\n"); -#endif - read_file_table_y(&prologue_table_y, NULL, &nb_prologue_table_y, false, false, y_kmin, y_kmax, u_size, y_size); -#ifdef PRINT_OUT - mexPrintf("nb_prologue_table_y=%d\n", nb_prologue_table_y); - mexPrintf("going to read first_table_y\n"); -#endif - read_file_table_y(&first_table_y, NULL, &nb_first_table_y, false, false, y_kmin, y_kmax, u_size, y_size); -#ifdef PRINT_OUT - mexPrintf("nb_first_table_y=%d\n", nb_first_table_y); - mexPrintf("going to read middle_table_y\n"); -#endif - read_file_table_y(&middle_table_y, &middle_i_table_y, &nb_middle_table_y, true, /*true*/false, y_kmin, y_kmax, u_size, y_size); -#ifdef PRINT_OUT - mexPrintf("nb_middle_table_y=%d\n", nb_middle_table_y); - mexPrintf("going to read last_table_y\n"); -#endif - read_file_table_y(&last_table_y, NULL, &nb_last_table_y, false, false, y_kmin, y_kmax, u_size, y_size); -#ifdef PRINT_OUT - mexPrintf("nb_last_table_y=%d\n", nb_last_table_y); - mexPrintf("->nb_last_table_y=%d\n", nb_last_table_y); - mexPrintf("max_u=%d\n",max_u); - mexPrintf("min_u=%d\n",min_u); -#endif - if (nb_last_table_u > 0) - { -#ifdef PRINT_OUT - mexPrintf("y_size=%d, periods=%d, y_kmin=%d, y_kmax=%d\n", y_size, periods, y_kmin, y_kmax); - mexPrintf("u=mxMalloc(%d)\n", max(u_count + 1,max_u+1)); -#endif - u = (double*)mxMalloc(max(u_count + 1,max_u+1) * sizeof(double)); - } - else - { -#ifdef PRINT_OUT - mexPrintf("u_size=%d, y_size=%d, periods=%d, y_kmin=%d, y_kmax=%d, u_count=%d, nb_add_u_count=%d\n", u_size, y_size, periods, y_kmin, y_kmax, u_count, nb_add_u_count); - mexPrintf("u=mxMalloc(%d)\n", u_count + (periods + y_kmin + y_kmax)* /*(u_count-u_size*(periods+y_kmin+y_kmax))*/nb_add_u_count); -#endif - u = (double*)mxMalloc((u_count + (periods + y_kmin + y_kmax)* /*(u_count-u_size*(periods+y_kmin+y_kmax)*/nb_add_u_count) * sizeof(double)); - memset(u, 0, (u_count + (periods + y_kmin + y_kmax)* /*(u_count-u_size*(periods+y_kmin+y_kmax)*/nb_add_u_count)*sizeof(double)); - } - if (u == NULL) - { - mexPrintf("memory exhausted\n"); - mexEvalString("st=fclose('all');clear all;"); - mexErrMsgTxt("Exit from Dynare"); - } - // mexErrMsgTxt("Exit from Dynare"); -} - - -void -SparseMatrix::close_SaveCode() -{ - if (SaveCode.is_open()) - SaveCode.close(); -} - - -/*void -SparseMatrix::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) -{ - periods=periods_arg; - nb_endo=nb_endo_arg; - y_kmin=y_kmin_arg; - y_kmax=y_kmax_arg; - y_size=y_size_arg; - u_count=u_count_arg; - u_count_init=u_count_init_arg; - u=u_arg; - y=y_arg; - ya=ya_arg; - slowc=slowc_arg; - slowc_save=slowc; - y_decal=y_decal_arg; - markowitz_c=markowitz_c_arg; - res1=res1_arg; - res2=res2_arg; - max_res=max_res_arg; -} -*/ diff --git a/mex/sources/simulate/simulate.cc b/mex/sources/simulate/simulate.cc index d34b3f6c6..0ea7a0087 100644 --- a/mex/sources/simulate/simulate.cc +++ b/mex/sources/simulate/simulate.cc @@ -42,28 +42,157 @@ max(int a, int b) +#ifdef DEBUG_EX +/*The Matlab c++ interface*/ +using namespace std; +#include +#include "mex_interface.hh" - - -/*class EvalException +int +main( int argc, const char* argv[] ) { -};*/ + FILE *fid; + printf("argc=%d\n",argc); + float f_tmp; + ostringstream tmp_out(""); + tmp_out << argv[1] << "_options.txt"; + cout << tmp_out.str().c_str() << "\n"; + int nb_params; + int i, row_y, col_y, row_x, col_x; + double *yd, *xd; + double *direction; + + string file_name(argv[1]); + //mexPrintf("file_name=%s\n",file_name.c_str()); + + fid = fopen(tmp_out.str().c_str(),"r"); + fscanf(fid,"%d",&periods); + fscanf(fid,"%d",&maxit_); + fscanf(fid,"%f",&f_tmp); + slowc = f_tmp; + //mexPrintf("slowc_save=%f\n",slowc_save); + fscanf(fid,"%f",&f_tmp); + markowitz_c = f_tmp; + fscanf(fid,"%f",&f_tmp); + solve_tolf = f_tmp; + fclose(fid); + + tmp_out.str(""); + tmp_out << argv[1] << "_M.txt"; + //printf("%s\n",tmp_out.str().c_str()); + fid = fopen(tmp_out.str().c_str(),"r"); + fscanf(fid,"%d",&y_kmin); + //printf("y_kmin=%d\n",y_kmin); + fscanf(fid,"%d",&y_kmax); + //printf("y_kmax=%d\n",y_kmax); + fscanf(fid,"%d",&y_decal); + //printf("y_decal=%d\n",y_decal); + fscanf(fid,"%d",&nb_params); + //printf("nb_params=%d\n",nb_params); + fscanf(fid,"%d",&row_x); + //printf("row_x=%d\n",row_x); + fscanf(fid,"%d",&col_x); + //printf("col_x=%d\n",col_x); + fscanf(fid,"%d",&row_y); + //printf("row_y=%d\n",row_y); + fscanf(fid,"%d",&col_y); + //printf("col_y=%d\n",col_y); + fscanf(fid,"%d",&nb_row_xd); + //printf("nb_row_xd=%d\n",nb_row_xd); + params = (double*)malloc(nb_params*sizeof(params[0])); + //printf("OK1\n"); + for(i=0; i < nb_params; i++) + { + fscanf(fid,"%f",&f_tmp); + params[i] = f_tmp; + //printf("param[%d]=%f\n",i,params[i]); + } + //printf("OK2\n"); + fclose(fid); + yd = (double*)malloc(row_y*col_y*sizeof(yd[0])); + xd = (double*)malloc(row_x*col_x*sizeof(xd[0])); + tmp_out.str(""); + tmp_out << argv[1] << "_oo.txt"; + fid = fopen(tmp_out.str().c_str(),"r"); + for(i=0; i < col_y*row_y; i++) + { + fscanf(fid,"%f",&f_tmp); + yd[i] = f_tmp; + } + for(i=0; i < col_x*row_x; i++) + { + fscanf(fid,"%f",&f_tmp); + xd[i] = f_tmp; + } + fclose(fid); + + size_of_direction=col_y*row_y*sizeof(double); + y=(double*)mxMalloc(size_of_direction); + ya=(double*)mxMalloc(size_of_direction); + direction=(double*)mxMalloc(size_of_direction); + memset(direction,0,size_of_direction); + x=(double*)mxMalloc(col_x*row_x*sizeof(double)); + for (i=0;i0) + { + plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL); + pind = mxGetPr(plhs[0]); + for (i=0;i0) - mexPrintf("\n%d %f ",i/row_x+1,x[i]); - else - mexPrintf("%f ",x[i]); - } - for(i=0;i0) { plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL); @@ -175,7 +283,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) for (i=0;i