diff --git a/tests/ferhat/bidon.m b/tests/ferhat/bidon.m
new file mode 100644
index 0000000000000000000000000000000000000000..88751b114ee02c769def22868bb2e677a05e2a35
--- /dev/null
+++ b/tests/ferhat/bidon.m
@@ -0,0 +1,176 @@
+%
+% Status : main Dynare file
+%
+% Warning : this file is generated automatically by Dynare
+% from model file (.mod)
+
+clear all
+tic;
+global M_ oo_ exedet_ exdet_ recur_ recurs_
+global options_ endval_
+global ys0_ recurs0_ ex0_ ct_
+options_ = [];
+M_.fname = 'bidon';
+%
+% Some global variables initialisation
+%
+global_initialization;
+diary off;
+warning off;
+
+delete bidon.log;
+warning on;
+warning backtrace;
+logname_ = 'bidon.log';
+diary 'bidon.log';
+erase_compiled_function('bidon_static');
+erase_compiled_function('bidon_dynamic');
+M_.exo_names = 'g_bar';
+M_.exo_names_tex = '';
+M_.endo_names = 'c';
+M_.endo_names_tex = '';
+M_.endo_names = strvcat(M_.endo_names,'i');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'g');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'infl');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'y');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'k');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'l');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'r');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'p1');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'q1');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'p2');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'q2');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.endo_names = strvcat(M_.endo_names,'r0');
+M_.endo_names_tex = strvcat(M_.endo_names_tex,'');
+M_.param_names = 'a';
+M_.param_names_tex = '';
+M_.param_names = strvcat(M_.param_names,'b');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'d');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'e');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'f');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'h');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'j');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'m');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'n');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'o');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'p');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.param_names = strvcat(M_.param_names,'q');
+M_.param_names_tex = strvcat(M_.param_names_tex,'');
+M_.exo_det_nbr = 0;
+M_.exo_nbr = 1;
+M_.Sigma_e = zeros(1, 1);
+M_.endo_nbr = 13;
+M_.recur_nbr = 0;
+M_.param_nbr = 12;
+M_.lead_lag_incidence = [
+ 0 2 5 18 0;
+ 0 0 6 0 0;
+ 0 0 7 0 0;
+ 0 0 8 0 0;
+ 1 3 9 19 20;
+ 0 4 10 0 0;
+ 0 0 11 0 0;
+ 0 0 12 0 0;
+ 0 0 13 0 0;
+ 0 0 14 0 0;
+ 0 0 15 0 0;
+ 0 0 16 0 0;
+ 0 0 17 0 0;]';
+M_.exo_names_orig_ord = [1:1];
+M_.maximum_lag = 10;
+M_.maximum_lead = 2;
+M_.maximum_endo_lag = 2;
+M_.maximum_endo_lead = 2;
+oo_.steady_state = zeros(13, 1);
+M_.maximum_exo_lag = 10;
+M_.maximum_exo_lead = 0;
+oo_.exo_steady_state = zeros(1, 1);
+M_.params = zeros(12, 1);
+M_.params( 1 ) = (0.4*0.6);
+a = M_.params( 1 );
+M_.params( 2 ) = (0.3*0.6);
+b = M_.params( 2 );
+M_.params( 3 ) = 0.1;
+d = M_.params( 3 );
+M_.params( 4 ) = 0.15;
+e = M_.params( 4 );
+M_.params( 5 ) = 1;
+f = M_.params( 5 );
+M_.params( 6 ) = 0.15;
+h = M_.params( 6 );
+M_.params( 7 ) = 1;
+j = M_.params( 7 );
+M_.params( 8 ) = 1;
+m = M_.params( 8 );
+M_.params( 9 ) = 1;
+n = M_.params( 9 );
+M_.params( 10 ) = 1;
+o = M_.params( 10 );
+%
+% INITVAL instructions
+%
+options_.initval_file = 0;
+endval_=0;
+oo_.exo_steady_state( 1 ) = 0.15;
+oo_.steady_state( 1 ) = 0.7;
+oo_.steady_state( 2 ) = 0.15;
+oo_.steady_state( 3 ) = 0.15;
+oo_.steady_state( 5 ) = 1;
+oo_.steady_state( 6 ) = 1;
+oo_.steady_state( 7 ) = 1;
+oo_.steady_state( 4 ) = 0.02;
+oo_.steady_state( 8 ) = 0;
+oo_.steady_state( 13 ) = oo_.steady_state(8);
+oo_.steady_state( 9 ) = (2/3);
+oo_.steady_state( 10 ) = (3.1/3);
+oo_.steady_state( 12 ) = (4/9);
+oo_.steady_state( 11 ) = (8/9);
+oo_.y_simul=[oo_.steady_state*ones(1,M_.maximum_lag)];
+if M_.exo_nbr > 0;
+ oo_.exo_simul = [ones(M_.maximum_lag,1)*oo_.exo_steady_state'];
+end;
+if M_.exo_det_nbr > 0;
+ oo_.exo_det_simul = [ones(M_.maximum_lag,1)*oo_.exo_det_steady_state'];
+end;
+options_.solve_algo = 2;
+steady;
+options_.periods = 100;
+options_.simul = 1;
+if ~ options_.initval_file
+ make_y_;
+ make_ex_;
+end
+disp('compiling...');
+mex bidon_dynamic.cc;
+oo_.endo_simul=bidon_dynamic;
+var_list_=[];
+var_list_ = strvcat(var_list_,'c');
+rplot(var_list_);
+var_list_=[];
+var_list_ = strvcat(var_list_,'y');
+rplot(var_list_);
+save('bidon_results', 'oo_');
+diary off
+
+disp(['Total computing time : ' sec2hms(round(toc)) ]);
diff --git a/tests/ferhat/bidon.mod b/tests/ferhat/bidon.mod
new file mode 100644
index 0000000000000000000000000000000000000000..f38e75afdd6944cfe0d083e419f6907b19110d52
--- /dev/null
+++ b/tests/ferhat/bidon.mod
@@ -0,0 +1,66 @@
+var c i g infl y k l r p1 q1 p2 q2 r0;
+varexo g_bar;
+
+parameters a b d e f h j m n o p q;
+a=0.4*0.6;
+b=0.3*0.6;
+d=0.1;
+e=0.15;
+f=1;
+h=0.15;
+j=1;
+m=1;
+n=1;
+o=1;
+
+model(SPARSE_DLL,gcc_compiler);
+/*0*/ k=(1-h)*k(-1)+i; /*k:0*/
+/*1*/ y=l^j*k^m; /*l:1*/
+/*2*/ c=y*a+b+0.3*c(-1)+0.1*c(+1)+0.*g_bar(-10); /*c:2*/
+/*3*/ infl=0.02*y+0.5*r; /*infl:3*/
+/*4*/ i=d*(y-y(-1))+e/**r*/; /*i4*/
+/*5*/ g=f*g_bar; /*g:5*/
+/*6*/ y=0.6*(c+i+g)+0.1*y(-2)+0.1*y(+2)+0.1*y(-1)+0.1*y(+1); /*y:6*/
+/*7*/ r=y-1+infl-0.02; /*r7*/
+/*8*/ p1=i+0.5*q1;
+/*9*/ q1=0.5*p1+c;
+/*10*/ q2=0.5*p2+r0;
+/*11*/ p2=0.5*q2+p1;
+/*12*/ r0=r;
+end;
+
+initval;
+g_bar=0.15;
+c=0.7;
+i=0.15;
+g=0.15;
+y=1;
+k=1;
+l=1;
+infl=0.02;
+r=0;
+r0=r;
+p1=2/3;
+q1=3.1/3;
+q2=4/9;
+p2=8/9;
+
+end;
+
+steady(solve_algo=2);
+//check;
+
+/*shocks;
+var g_bar;
+periods 1;
+values 0.16;
+end;*/
+
+
+
+
+simul(periods=100);
+
+
+rplot c;
+rplot y;
diff --git a/tests/ferhat/pctimer_h.hh b/tests/ferhat/pctimer_h.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a3db34be17c38928b982135ebf9f235f127fd8e4
--- /dev/null
+++ b/tests/ferhat/pctimer_h.hh
@@ -0,0 +1,57 @@
+/*
+ * pctimer.h 1.1 02/03/05
+ *
+ * Uses Win32 performance counter functions to get a high-resolution timer
+ *
+ * By Wu Yongwei
+ *
+ */
+
+#ifndef _PCTIMER_H
+
+typedef double pctimer_t;
+
+#if defined(_WIN32) || defined(__CYGWIN__)
+
+#ifndef _WIN32
+#define PCTIMER_NO_WIN32
+#endif /* _WIN32 */
+
+#define WIN32_LEAN_AND_MEAN
+#include <windows.h>
+
+#ifdef PCTIMER_NO_WIN32
+#undef PCTIMER_NO_WIN32
+#undef _WIN32
+#endif /* PCTIMER_NO_WIN32 */
+
+__inline pctimer_t pctimer()
+{
+ static LARGE_INTEGER pcount, pcfreq;
+ static int initflag;
+
+ if (!initflag)
+ {
+ QueryPerformanceFrequency(&pcfreq);
+ initflag++;
+ }
+
+ QueryPerformanceCounter(&pcount);
+ return (double)pcount.QuadPart / (double)pcfreq.QuadPart;
+}
+
+#else /* Win32/Cygwin */
+
+#include <sys/time.h>
+
+__inline pctimer_t pctimer()
+{
+ struct timeval tv;
+ gettimeofday(&tv, NULL);
+ return (double)tv.tv_sec + (double)tv.tv_usec / 1000000;
+}
+
+#endif /* Win32/Cygwin */
+
+#endif /* _PCTIMER_H */
+
diff --git a/tests/ferhat/ramst.mod b/tests/ferhat/ramst.mod
new file mode 100644
index 0000000000000000000000000000000000000000..6d7aa8bf915b705cb698385c6bf8268eaa024836
--- /dev/null
+++ b/tests/ferhat/ramst.mod
@@ -0,0 +1,36 @@
+var c k;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+
+model/*(SPARSE_DLL,GCC_COMPILER)*/;
+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;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k;
+end;
+
+steady(solve_algo=2);
+
+//check;
+
+shocks;
+var x;
+periods 1;
+values 1.2;
+end;
+
+simul(periods=200);
+
+rplot c;
+rplot k;
diff --git a/tests/ferhat/ramst_a.mod b/tests/ferhat/ramst_a.mod
new file mode 100644
index 0000000000000000000000000000000000000000..4a81935c3a5283dc7ededbbce917696c8608bffe
--- /dev/null
+++ b/tests/ferhat/ramst_a.mod
@@ -0,0 +1,37 @@
+// check shocks on several periods
+var c k;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+
+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;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k +1 ;
+end;
+
+steady;
+
+check;
+
+shocks;
+var x;
+periods 1 2 3 4;
+values 1.1 1.2 1.3 1.4;
+end;
+
+simul(periods=200);
+
+rplot c;
+rplot k;
diff --git a/tests/ferhat/simulate.c b/tests/ferhat/simulate.c
new file mode 100644
index 0000000000000000000000000000000000000000..5a824caada2add9b7e35f18682682e2f3ce1ced2
--- /dev/null
+++ b/tests/ferhat/simulate.c
@@ -0,0 +1,788 @@
+ ////////////////////////////////////////////////////////////////////////
+ // simulate.c //
+ // simulate file designed for Matlab LCC C compiler //
+ // use LCC_COMPILER option in MODEL command [defailt option] //
+ ////////////////////////////////////////////////////////////////////////
+
+
+#define PRINT_OUT
+//#define FIXE
+#define SIZE_OF_INT sizeof(int)
+
+typedef struct t_table_y
+{
+ int index, nb;
+ int *u_index, *y_index;
+} t_table_y;
+
+typedef struct t_table_u
+{
+ struct t_table_u* pNext;
+ unsigned char type;
+ int index;
+ int op1, op2;
+} t_table_u;
+
+FILE *SaveCode = NULL, *fopen();
+
+
+void
+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)
+{
+ char type;
+ int i;
+ fread(nb_table_u, SIZE_OF_INT, 1, SaveCode);
+#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;
+ for(i = 0;i < *nb_table_u;i++)
+ {
+ fread(&type, sizeof(type), 1, SaveCode);
+ 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;
+ fread(&(*table_u)->index, SIZE_OF_INT, 1, SaveCode);
+ fread(&(*table_u)->op1, SIZE_OF_INT, 1, SaveCode);
+ if(shifting)
+ {
+ (*table_u)->index -= y_kmin * u_size;
+ (*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
+ 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;
+ fread(&(*table_u)->index, SIZE_OF_INT, 1, SaveCode);
+ fread(&(*table_u)->op1, SIZE_OF_INT, 1, SaveCode);
+ fread(&(*table_u)->op2, SIZE_OF_INT, 1, SaveCode);
+ if(shifting)
+ {
+ (*table_u)->index -= y_kmin * u_size;
+ (*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
+ 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;
+ fread(&(*table_u)->index, SIZE_OF_INT, 1, SaveCode);
+ if(shifting)
+ (*table_u)->index -= y_kmin * u_size;
+#ifdef PRINT_OUT
+ mexPrintf("push(u[%d])\n", (*table_u)->index);
+#endif
+ break;
+ }
+ }
+ 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++)
+ {
+ fread(&type, sizeof(type), 1, SaveCode);
+ switch (type)
+ {
+ case 3:
+ case 7:
+ (*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;
+ fread(&(*i_table_u)->index, SIZE_OF_INT, 1, SaveCode);
+ fread(&(*i_table_u)->op1, SIZE_OF_INT, 1, SaveCode);
+#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("u[%d]=1/(1-u[%d])\n", (*i_table_u)->index, (*i_table_u)->op1);
+ else
+ mexPrintf("u[%d]*=u[%d]\n", (*i_table_u)->index, (*i_table_u)->op1);
+#endif
+ break;
+ case 1:
+ case 2:
+ case 6:
+ (*i_table_u)->pNext = (t_table_u*)mxMalloc(sizeof(t_table_u));
+ (*i_table_u) = (*i_table_u)->pNext;
+ (*i_table_u)->type = type;
+ fread(&(*i_table_u)->index, SIZE_OF_INT, 1, SaveCode);
+ fread(&(*i_table_u)->op1, SIZE_OF_INT, 1, SaveCode);
+ fread(&(*i_table_u)->op2, SIZE_OF_INT, 1, SaveCode);
+#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("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("u[%d]+=u[%d]*u[%d]\n", (*i_table_u)->index, (*i_table_u)->op1, (*i_table_u)->op2);
+ else
+ mexPrintf("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:
+ (*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;
+ fread(&(*i_table_u)->index, SIZE_OF_INT, 1, SaveCode);
+#ifdef FIXE
+ (*i_table_u)->index = u_size;
+#endif
+#ifdef PRINT_OUT
+ mexPrintf("push(u[%d])\n", (*i_table_u)->index);
+#endif
+ break;
+ }
+ }
+ }
+}
+
+void
+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 i, k;
+ fread(nb_table_y, SIZE_OF_INT, 1, SaveCode);
+
+ mexPrintf("nb_table_y=%d\n", *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++)
+ {
+ fread(&((*table_y)[i].index), SIZE_OF_INT, 1, SaveCode);
+ fread(&((*table_y)[i].nb), SIZE_OF_INT, 1, SaveCode);
+ 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++)
+ {
+ fread(&((*table_y)[i].u_index[k]), SIZE_OF_INT, 1, SaveCode);
+ fread(&((*table_y)[i].y_index[k]), SIZE_OF_INT, 1, SaveCode);
+ /*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
+ //mexPrintf("sizeof((*i_table_y)[i].y_index[k])=%d\n",sizeof((*i_table_y)[i].y_index[k]));
+ 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);
+ mexPrintf("i_to_do=%d\n", i_to_do);
+#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++)
+ {
+ fread(&((*i_table_y)[i].index), SIZE_OF_INT, 1, SaveCode);
+ fread(&((*i_table_y)[i].nb), SIZE_OF_INT, 1, SaveCode);
+#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++)
+ {
+ fread(&((*i_table_y)[i].u_index[k]), SIZE_OF_INT, 1, SaveCode);
+ fread(&((*i_table_y)[i].y_index[k]), SIZE_OF_INT, 1, SaveCode);
+#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
+ }
+ }
+ }
+}
+
+
+int i, j, k, nb_endo, u_count, u_count_init;
+int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
+int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
+int first_count_loop, middle_count_loop;
+char type;
+t_table_u *prologue_table_u, *first_table_u, *first_i_table_u, *middle_table_u, *middle_i_table_u, *last_table_u;
+t_table_y *prologue_table_y, *first_table_y, *middle_table_y, *middle_i_table_y, *last_table_y;
+t_table_u *F_prologue_table_u, *F_first_table_u, *F_first_i_table_u, *F_middle_table_u, *F_middle_i_table_u, *F_last_table_u;
+
+
+void
+Read_file(char* a_file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax)
+{
+ int nb_add_u_count = 0;
+ char *file_name=(char*)malloc(strlen(a_file_name)+5);
+ file_name=strcat(a_file_name,".bin");
+ u_size = u_size1;
+#ifdef PRINT_OUT
+ mexPrintf("file_name=%s\n", file_name);
+#endif
+ if(!SaveCode)
+ {
+#ifdef PRINT_OUT
+ mexPrintf("file opened\n");
+#endif
+ if(!(SaveCode=fopen(file_name,"r")))
+ {
+ mexPrintf("Error : Can't open file \"%s\" for reading\n", file_name);
+ mexErrMsgTxt("Exit from Dynare");
+ }
+#ifdef PRINT_OUT
+ mexPrintf("done\n");
+#endif
+ }
+ fread(&nb_endo, SIZE_OF_INT, 1, SaveCode);
+ fread(&u_count, SIZE_OF_INT, 1, SaveCode);
+ fread(&u_count_init, SIZE_OF_INT, 1, SaveCode);
+
+#ifdef PRINT_OUT
+ mexPrintf("nb_endo=%d, u_count=%d, u_count_init=%d\n",nb_endo,u_count, u_count_init);
+ 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);
+#ifdef PRINT_OUT
+ mexPrintf("nb_prologue_table_u=%d\n", nb_prologue_table_u);
+ 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);
+#ifdef PRINT_OUT
+ mexPrintf("nb_first_table_u=%d\n", nb_first_table_u);
+ mexPrintf("nb_endo=%d\n", nb_endo);
+ mexPrintf("nb_first_table_u=%d\n", nb_first_table_u);
+ mexPrintf("u_count=%d\n", u_count);
+#endif
+ fread(&middle_count_loop, SIZE_OF_INT, 1, SaveCode);
+#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);
+#ifdef PRINT_OUT
+ 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);
+#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);
+#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);
+#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);
+#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);
+#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);
+#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", u_count + 1);
+#endif
+ u = (double*)mxMalloc((u_count + 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");
+ mexErrMsgTxt("Exit from Dynare");
+ }
+}
+
+
+typedef struct t_Stack
+{
+ double uu;
+ struct t_Stack *pPrev;
+} t_Stack;
+
+t_Stack *Stack=NULL, *tmp_Stack;
+
+void push(double uu)
+{
+ if(!Stack)
+ {
+ Stack=(t_Stack*)malloc(sizeof(t_Stack));
+ Stack->pPrev=NULL;
+ }
+ else
+ {
+ tmp_Stack=Stack;
+ Stack=(t_Stack*)malloc(sizeof(t_Stack));
+ Stack->pPrev=tmp_Stack;
+ }
+ Stack->uu=uu;
+}
+
+double pop()
+{
+ double uu;
+ uu=Stack->uu;
+ Stack=Stack->pPrev;
+ return(uu);
+}
+
+void
+simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax)
+{
+ int i, j, k, l, m, m1, nop;
+ int period = it_ * y_size, s_middle_count_loop = 0 ;
+ pctimer_t t1 = pctimer(), t2;
+ double uu;
+#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, u[j*u_size + i]);
+ }
+ mexPrintf("\n");
+ }
+#endif
+ if(nb_first_table_u > 0)
+ {
+ first_count_loop = it_;
+ s_middle_count_loop = it_ -y_kmin - middle_count_loop + 1;
+//#ifdef PRINT_OUT
+ mexPrintf("----------------------------------------------------------------------\n");
+ mexPrintf("// Simulate \\\n");
+ mexPrintf("----------------------------------------------------------------------\n");
+//#endif
+ }
+ nop = 0;
+ for(j = 0 ;j < first_count_loop - y_kmin;j++)
+ {
+ 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
+ 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, 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
+ 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, 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
+ 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, u[first_table_u->index + j*first_i_table_u->index]);
+#endif
+ break;
+ case 5:
+ push(u[first_table_u->index + j*first_i_table_u->index]);
+#ifdef PRINT_OUT
+ 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
+ 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, 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
+ 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, u[first_table_u->index + j*first_i_table_u->index]);
+#endif
+ break;
+ }
+ first_table_u = first_table_u->pNext;
+ first_i_table_u = first_i_table_u->pNext;
+ nop++;
+ }
+ }
+#ifdef PRINT_OUT
+ mexPrintf("prologue\n");
+#endif
+ 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
+ mexPrintf("u[%d]=u[%d]*u[%d]=%f\n", prologue_table_u->index , prologue_table_u->op1 , prologue_table_u->op2 , 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
+ mexPrintf("u[%d]+=u[%d]*u[%d]=%f\n" , prologue_table_u->index , prologue_table_u->op1 , prologue_table_u->op2 , u[prologue_table_u->index ]);
+#endif
+ break;
+ case 3:
+ u[prologue_table_u->index ] = 1 / ( -u[prologue_table_u->op1 ]);
+#ifdef PRINT_OUT
+ mexPrintf("u[%d]=1/(-u[%d])=%f\n", prologue_table_u->index, prologue_table_u->op1, u[prologue_table_u->index]);
+#endif
+ break;
+ case 5:
+ push(u[prologue_table_u->index]);
+#ifdef PRINT_OUT
+ 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
+ mexPrintf("u[%d]=1/(1-u[%d]*u[%d])=%f\n", prologue_table_u->index, prologue_table_u->op1, prologue_table_u->op2, u[prologue_table_u->index]);
+#endif
+ break;
+ case 7:
+ u[prologue_table_u->index] *= u[prologue_table_u->op1];
+#ifdef PRINT_OUT
+ mexPrintf("u[%d]*=u[%d]=%f\n", prologue_table_u->index, prologue_table_u->op1, u[prologue_table_u->index]);
+#endif
+ break;
+ }
+ prologue_table_u = prologue_table_u->pNext;
+ nop++;
+ }
+#ifdef PRINT_OUT
+ mexPrintf("middle_u (s_middle_count_loop=%d\n", s_middle_count_loop);
+#endif
+ for(j = 0;j < s_middle_count_loop - y_kmin;j++)
+ {
+ //cout << "j=" << j << "\n";
+#ifdef PRINT_OUT
+ 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
+ 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, 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
+ 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, 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
+ 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, u[middle_table_u->index + j*middle_i_table_u->index]);
+#endif
+ break;
+ case 5:
+#ifdef PRINT_OUT
+ 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
+ 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
+ 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, 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
+ 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, u[middle_table_u->index + j*middle_i_table_u->index]);
+#endif
+ break;
+ }
+ 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
+ mexPrintf("u[%d]=u[%d]*u[%d]=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, 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
+ mexPrintf("u[%d]+=u[%d]*u[%d]=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, u[last_table_u->index]);
+#endif
+ break;
+ case 3:
+ u[last_table_u->index] = 1 / ( -u[last_table_u->op1]);
+#ifdef PRINT_OUT
+ mexPrintf("u[%d]=1/(-u[%d])=%f\n", last_table_u->index, last_table_u->op1, u[last_table_u->index]);
+#endif
+ break;
+ case 5:
+ push(u[last_table_u->index]);
+#ifdef PRINT_OUT
+ 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
+ mexPrintf("u[%d]=1/(1-u[%d]*u[%d])=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, u[last_table_u->index]);
+#endif
+ break;
+ case 7:
+ u[last_table_u->index] *= u[last_table_u->op1];
+#ifdef PRINT_OUT
+ mexPrintf("u[%d]*=u[%d]=%f\n", last_table_u->index, last_table_u->op1, u[last_table_u->index]);
+#endif
+ break;
+ }
+ last_table_u = last_table_u->pNext;
+ nop++;
+ }
+#ifdef PRINT_OUT
+ mexPrintf("nb_last_table_y=%d\n",nb_last_table_y);
+#endif
+ for(i = nb_last_table_y - 1;i >= 0;i--)
+ {
+ k = last_table_y[i].index;
+ y[period + k] = 0;
+#ifdef PRINT_OUT
+ mexPrintf("it_=%d\n", it_);
+ mexPrintf("->y[it_*y_size+%d]=y[%d]=", k, it_*y_size + k);
+#endif
+ for(j = last_table_y[i].nb - 1;j >= 0;j--)
+ {
+ uu = pop();
+ m = last_table_y[i].y_index[j];
+#ifdef PRINT_OUT
+ if(j > 0)
+ {
+ if(m >= 0)
+ mexPrintf("u[%d](%f)*y[%d]+", last_table_y[i].u_index[j], uu, last_table_y[i].y_index[j]);
+ else
+ mexPrintf("u[%d](%f)+", last_table_y[i].u_index[j], uu);
+ }
+ else
+ {
+ if(m >= 0)
+ mexPrintf("u[%d](%f)*y[%d]", last_table_y[i].u_index[j], uu, last_table_y[i].y_index[j]);
+ else
+ mexPrintf("u[%d](%f)", last_table_y[i].u_index[j], uu);
+ }
+#endif
+ if(m >= 0)
+ y[period + k] += uu * y[period + m];
+ else
+ y[period + k] += uu;
+ }
+#ifdef PRINT_OUT
+ mexPrintf("=%f\n", y[period + k]);
+#endif
+ nop++;
+ }
+#ifdef PRINT_OUT
+ mexPrintf("s_middle_count_loop=%d\n",s_middle_count_loop);
+ mexPrintf("nb_middle_table_y=%d\n",nb_middle_table_y);
+#endif
+ for(j = s_middle_count_loop - y_kmin - 1;j >= 0;j--)
+ {
+ for(i = nb_middle_table_y - 1;i >= 0;i--)
+ {
+ k = middle_table_y[i].index + j * middle_i_table_y[i].index;
+ y[k] = 0;
+#ifdef PRINT_OUT
+ mexPrintf("y[%d]=", k );
+#endif
+ for(l = middle_table_y[i].nb - 1;l >= 0;l--)
+ {
+ uu = pop();
+ m = middle_table_y[i].y_index[l] + j * middle_i_table_y[i].y_index[l];
+#ifdef PRINT_OUT
+ if(m >= 0)
+ {
+ m1 = middle_table_y[i].u_index[l] + j * middle_i_table_y[i].u_index[l];
+ if(l > 0)
+ mexPrintf("u[%d](%f)*y[%d](%f)+", m1, uu, m, y[m]);
+ else
+ mexPrintf("u[%d](%f)*y[%d](%f)", m1, uu, m, y[m]);
+ }
+ else
+ {
+ m1 = middle_table_y[i].u_index[l] + j * middle_i_table_y[i].u_index[l];
+ if(l > 0)
+ mexPrintf("u[%d](%f)*y[%d](%f)+", m1, uu, m, 1.0);
+ else
+ mexPrintf("u[%d](%f)*y[%d](%f)", m1, uu, m, 1.0);
+ }
+#endif
+ if(m >= 0)
+ y[k] += uu * y[m];
+ else
+ y[k] += uu;
+ }
+#ifdef PRINT_OUT
+ mexPrintf("=%f\n", y[k]);
+#endif
+ nop++;
+ }
+ }
+ for(i = nb_prologue_table_y - 1;i >= 0;i--)
+ {
+ k = prologue_table_y[i].index;
+ y[k] = 0;
+#ifdef PRINT_OUT
+ mexPrintf("y[%d]=", k);
+#endif
+ for(j = prologue_table_y[i].nb - 1;j >= 0;j--)
+ {
+ uu = pop();
+#ifdef PRINT_OUT
+ 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], uu, prologue_table_y[i].y_index[j], y[prologue_table_y[i].y_index[j]]);
+ else
+ mexPrintf("u[%d](%f)*y[%d](%f)", prologue_table_y[i].u_index[j], uu, prologue_table_y[i].y_index[j], 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], 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], uu, prologue_table_y[i].y_index[j], 1.0);
+ }
+#endif
+ if(prologue_table_y[i].y_index[j] >= 0)
+ y[k] += uu * y[prologue_table_y[i].y_index[j]];
+ else
+ y[k] += uu;
+ }
+#ifdef PRINT_OUT
+ mexPrintf("=%f\n", y[k]);
+#endif
+ nop++;
+ }
+ for(i = nb_first_table_y - 1;i >= 0;i--)
+ {
+ k = first_table_y[i].index;
+ y[k] = 0;
+#ifdef PRINT_OUT
+ mexPrintf("y[%d]=", k);
+#endif
+ for(j = first_table_y[i].nb - 1;j >= 0;j--)
+ {
+ uu = pop();
+#ifdef PRINT_OUT
+ if(j > 0)
+ mexPrintf("u[%d](%f)*y[%d](%f)+", first_table_y[i].u_index[j], uu, first_table_y[i].y_index[j], y[first_table_y[i].y_index[j]]);
+ else
+ mexPrintf("u[%d](%f)*y[%d](%f)", first_table_y[i].u_index[j], uu, first_table_y[i].y_index[j], y[first_table_y[i].y_index[j]]);
+#endif
+ if(m >= 0)
+ y[k] += uu * y[first_table_y[i].y_index[j]];
+ else
+ y[k] += uu;
+ }
+#ifdef PRINT_OUT
+ mexPrintf("=%f\n", y[k]);
+#endif
+ nop++;
+ }
+ t2 = pctimer();
+ if(nb_first_table_u > 0)
+ mexPrintf("(**%f milliseconds u_count : %d nop : %d **)\n", 1000*(t2 - t1), u_count, nop);
+}
+
diff --git a/tests/ferhat/simulate.cc b/tests/ferhat/simulate.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b6fcd46d3a3c88226ae279e02ed62ec42703e3c1
--- /dev/null
+++ b/tests/ferhat/simulate.cc
@@ -0,0 +1,850 @@
+ ////////////////////////////////////////////////////////////////////////
+ // simulate.c //
+ // simulate file designed for GNU GCC C++ compiler //
+ // use GCC_COMPILER option in MODEL command //
+ ////////////////////////////////////////////////////////////////////////
+
+//#define PRINT_OUT
+//#define PRINT_OUT_p
+#include <stack>
+#include <math.h>
+#include <iostream>
+#include <fstream>
+#include "pctimer_h.hh"
+#include "mex.h"
+
+
+/*#include "simulate.hh"*/
+
+int Per_y_, Per_u_, it_, nb_row_x, u_size, y_size, x_size, y_kmin, y_kmax, y_decal;
+int periods, maxit_;
+double *params, *u, *y, *x, *r, *g1, *g2;
+double slowc, solve_tolf, max_res, res1, res2;
+bool cvg;
+pctimer_t t0, t1;
+
+typedef struct t_table_y
+{
+ int index, nb;
+ int *u_index, *y_index;
+};
+
+typedef struct t_table_u
+{
+ t_table_u* pNext;
+ unsigned char type;
+ int index;
+ int op1, op2;
+};
+
+
+std::ifstream SaveCode;
+
+
+
+
+void
+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)
+{
+ char type;
+ int i;
+ SaveCode.read(reinterpret_cast<char *>(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<char *>(&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<char *>(&(*table_u)->index), sizeof((*table_u)->index));
+ SaveCode.read(reinterpret_cast<char *>(&(*table_u)->op1), sizeof((*table_u)->op1));
+ if(shifting)
+ {
+ (*table_u)->index -= y_kmin * u_size;
+ (*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<char *>(&(*i_table_u)->index), sizeof((*i_table_u)->index));
+ SaveCode.read(reinterpret_cast<char *>(&(*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<char *>(&(*table_u)->index), sizeof((*table_u)->index));
+ SaveCode.read(reinterpret_cast<char *>(&(*table_u)->op1), sizeof((*table_u)->op1));
+ SaveCode.read(reinterpret_cast<char *>(&(*table_u)->op2), sizeof((*table_u)->op2));
+ if(shifting)
+ {
+ (*table_u)->index -= y_kmin * u_size;
+ (*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<char *>(&(*i_table_u)->index), sizeof((*i_table_u)->index));
+ SaveCode.read(reinterpret_cast<char *>(&(*i_table_u)->op1), sizeof((*i_table_u)->op1));
+ SaveCode.read(reinterpret_cast<char *>(&(*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<char *>(&(*table_u)->index), sizeof((*table_u)->index));
+ if(shifting)
+ (*table_u)->index -= y_kmin * u_size;
+#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<char *>(&(*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
+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 i, k;
+ SaveCode.read(reinterpret_cast<char *>(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<char *>(&((*table_y)[i].index)), sizeof((*table_y)[i].index));
+ SaveCode.read(reinterpret_cast<char *>(&((*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<char *>(&((*table_y)[i].u_index[k])), sizeof((*table_y)[i].u_index[k]));
+ SaveCode.read(reinterpret_cast<char *>(&((*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<char *>(&((*i_table_y)[i].index)), sizeof((*i_table_y)[i].index));
+ SaveCode.read(reinterpret_cast<char *>(&((*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<char *>(&((*i_table_y)[i].u_index[k])), sizeof((*i_table_y)[i].u_index[k]));
+ SaveCode.read(reinterpret_cast<char *>(&((*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
+ }
+ }
+ }
+}
+
+
+int i, j, k, nb_endo, u_count, u_count_init, iter;
+int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
+int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
+int first_count_loop, middle_count_loop;
+char type;
+t_table_u *prologue_table_u, *first_table_u, *first_i_table_u, *middle_table_u, *middle_i_table_u, *last_table_u;
+t_table_y *prologue_table_y, *first_table_y, *middle_table_y, *middle_i_table_y, *last_table_y;
+t_table_u *F_prologue_table_u, *F_first_table_u, *F_first_i_table_u, *F_middle_table_u, *F_middle_i_table_u, *F_last_table_u;
+std::string filename;
+
+
+void
+Read_file(std::string file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax)
+{
+ int nb_add_u_count = 0;
+ u_size = u_size1;
+ filename=file_name;
+#ifdef PRINT_OUT
+ 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());
+ mexErrMsgTxt("Exit from Dynare");
+ }
+#ifdef PRINT_OUT
+ mexPrintf("done\n");
+#endif
+ }
+ SaveCode.read(reinterpret_cast<char *>(&nb_endo), sizeof(nb_endo));
+ SaveCode.read(reinterpret_cast<char *>(&u_count), sizeof(u_count));
+ SaveCode.read(reinterpret_cast<char *>(&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);
+#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);
+#ifdef PRINT_OUT
+ mexPrintf("nb_prologue_table_u=%d\n", nb_prologue_table_u);
+#endif
+ //mexErrMsgTxt("Exit from Dynare");
+ SaveCode.read(reinterpret_cast<char *>(&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);
+#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);
+#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);
+#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);
+#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);
+#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);
+#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);
+#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", u_count + 1);
+#endif
+ u = (double*)mxMalloc((u_count + 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");
+ mexErrMsgTxt("Exit from Dynare");
+ }
+ // mexErrMsgTxt("Exit from Dynare");
+}
+
+
+std::stack <double> Stack;
+
+void
+simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax)
+{
+ int i, j, k, l, m, m1, nop;
+ int period = it_ * y_size, s_middle_count_loop = 0 ;
+ pctimer_t t1 = pctimer();
+ double uu, yy;
+ char tmp_s[150];
+#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, u[j*u_size + i]);
+ }
+ mexPrintf("\n");
+ }
+#endif
+ if(nb_first_table_u > 0)
+ {
+ first_count_loop = it_;
+ s_middle_count_loop = it_ -y_kmin - middle_count_loop + 1;
+//#ifdef PRINT_OUT
+ mexPrintf("----------------------------------------------------------------------\n");
+ mexPrintf(" Simulate it�ration� %d \n",iter+1);
+ mexPrintf(" max. error=%.10e \n",max_res);
+ mexPrintf(" sqr. error=%.10e \n",res2);
+ mexPrintf(" abs. error=%.10e \n",res1);
+ mexPrintf("----------------------------------------------------------------------\n");
+//#endif
+ }
+ nop = 0;
+ for(j = 0 ;j < first_count_loop - y_kmin;j++)
+ {
+ 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
+ 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, 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
+ 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, 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
+ 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, 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
+ 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
+ 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, 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
+ 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, 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));
+ 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_p
+ 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_p
+ mexPrintf("u[%d]=u[%d]*u[%d]=%f\n", prologue_table_u->index , prologue_table_u->op1 , prologue_table_u->op2 , 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_p
+ mexPrintf("u[%d]+=u[%d]*u[%d]=%f\n" , prologue_table_u->index , prologue_table_u->op1 , prologue_table_u->op2 , u[prologue_table_u->index ]);
+#endif
+ break;
+ case 3:
+ u[prologue_table_u->index ] = 1 / ( -u[prologue_table_u->op1 ]);
+#ifdef PRINT_OUT_p
+ mexPrintf("u[%d]=1/(-u[%d])=%f\n", prologue_table_u->index, prologue_table_u->op1, u[prologue_table_u->index]);
+#endif
+ break;
+ case 5:
+ //nb_prologue_push++;
+ Stack.push(u[prologue_table_u->index]);
+#ifdef PRINT_OUT_p
+ 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_p
+ mexPrintf("u[%d]=1/(1-u[%d]*u[%d])=%f\n", prologue_table_u->index, prologue_table_u->op1, prologue_table_u->op2, u[prologue_table_u->index]);
+#endif
+ break;
+ case 7:
+ u[prologue_table_u->index] *= u[prologue_table_u->op1];
+#ifdef PRINT_OUT_p
+ mexPrintf("u[%d]*=u[%d]=%f\n", prologue_table_u->index, prologue_table_u->op1, 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] (in prologue_table_u)",prologue_table_u->index);
+ filename+=" stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
+ prologue_table_u = prologue_table_u->pNext;
+ nop++;
+ }
+#ifdef PRINT_OUT
+ 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
+ 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
+ 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, 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
+ 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, 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
+ 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, u[middle_table_u->index + j*middle_i_table_u->index]);
+#endif
+ break;
+ case 5:
+#ifdef PRINT_OUT
+ 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
+ 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, 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
+ 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, 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 (in middle_table_u)",middle_table_u->index,j);
+ 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
+ mexPrintf("u[%d]=u[%d]*u[%d]=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, 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
+ mexPrintf("u[%d]+=u[%d]*u[%d]=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, u[last_table_u->index]);
+#endif
+ break;
+ case 3:
+ u[last_table_u->index] = 1 / ( -u[last_table_u->op1]);
+#ifdef PRINT_OUT
+ mexPrintf("u[%d]=1/(-u[%d])=%f\n", last_table_u->index, last_table_u->op1, u[last_table_u->index]);
+#endif
+ break;
+ case 5:
+ Stack.push(u[last_table_u->index]);
+#ifdef PRINT_OUT
+ 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
+ mexPrintf("u[%d]=1/(1-u[%d]*u[%d])=%f\n", last_table_u->index, last_table_u->op1, last_table_u->op2, u[last_table_u->index]);
+#endif
+ break;
+ case 7:
+ u[last_table_u->index] *= u[last_table_u->op1];
+#ifdef PRINT_OUT
+ mexPrintf("u[%d]*=u[%d]=%f\n", last_table_u->index, last_table_u->op1, 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] (in last_table_u)",last_table_u->index);
+ filename+=" stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
+ last_table_u = last_table_u->pNext;
+ nop++;
+ }
+ 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
+ mexPrintf("it_=%d\n", it_);
+ mexPrintf("->y[it_*y_size+%d]=y[%d]=", k, it_*y_size + 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
+ if(j > 0)
+ {
+ if(m >= 0)
+ mexPrintf("u[%d](%f)*y[%d]+", last_table_y[i].u_index[j], uu, last_table_y[i].y_index[j]);
+ else
+ mexPrintf("u[%d](%f)+", last_table_y[i].u_index[j], uu);
+ }
+ else
+ {
+ if(m >= 0)
+ mexPrintf("u[%d](%f)*y[%d]", last_table_y[i].u_index[j], uu, last_table_y[i].y_index[j]);
+ else
+ mexPrintf("u[%d](%f)", last_table_y[i].u_index[j], 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);
+ filename+=" stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
+ /*if(((k-73) % y_size)==0)
+ mexPrintf("y[it_*y_size +73]=%f \n",yy);*/
+ y[period + k] += slowc*(yy - y[period + k]);
+//#ifdef PRINT_OUT
+ mexPrintf("=%f\n", y[period + k]);
+//#endif
+ nop++;
+ }
+ int nb_middle_pop=0;
+ for(j = s_middle_count_loop - y_kmin - 1+y_decal;j >= y_decal;j--)
+ {
+ for(i = nb_middle_table_y - 1;i >= 0;i--)
+ {
+ k = middle_table_y[i].index + j * middle_i_table_y[i].index;
+ yy = 0;
+//#ifdef PRINT_OUT
+ mexPrintf("(0)y[%d]=", k );
+//#endif
+ for(l = middle_table_y[i].nb - 1;l >= 0;l--)
+ {
+ uu = Stack.top();
+ //mexPrintf("{");
+ Stack.pop();
+ //mexPrintf("}");
+ //nb_middle_pop++;
+ m = middle_table_y[i].y_index[l] + j * middle_i_table_y[i].y_index[l];
+//#ifdef PRINT_OUT
+ if(m >= 0)
+ {
+ m1 = middle_table_y[i].u_index[l] + j * middle_i_table_y[i].u_index[l];
+ if(l > 0)
+ mexPrintf("u[%d](%f)*y[%d](%f)+", m1, uu, m, y[m]);
+ else
+ mexPrintf("u[%d](%f)*y[%d](%f)", m1, uu, m, y[m]);
+ }
+ else
+ {
+ m1 = middle_table_y[i].u_index[l] + j * middle_i_table_y[i].u_index[l];
+ if(l > 0)
+ mexPrintf("u[%d](%f)*y[%d](%f)+", m1, uu, m, 1.0);
+ else
+ mexPrintf("u[%d](%f)*y[%d](%f)", m1, uu, m, 1.0);
+ }
+//#endif
+ if(m >= 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);
+ filename+=" stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
+ /*if(((k-73) % y_size)==0)
+ mexPrintf("y[it_*y_size +73]=%f \n",yy);*/
+ y[k] += slowc*(yy - y[k]);
+//#ifdef PRINT_OUT
+ mexPrintf("=%f\n", y[k]);
+//#endif
+ nop++;
+ }
+ }
+ //mexPrintf("nb_middle_push=%d et nb_middle_pop=%d\n",nb_middle_push,nb_middle_pop);
+ //int nb_prologue_pop=0;
+ for(i = nb_prologue_table_y - 1;i >= 0;i--)
+ {
+ k = prologue_table_y[i].index;
+ yy = 0;
+//#ifdef PRINT_OUT_p
+ 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_p
+ 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], uu, prologue_table_y[i].y_index[j], y[prologue_table_y[i].y_index[j]]);
+ else
+ mexPrintf("u[%d](%f)*y[%d](%f)", prologue_table_y[i].u_index[j], uu, prologue_table_y[i].y_index[j], 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], 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], 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);
+ filename+=" stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
+ /*if(((k-73) % y_size)==0)
+ mexPrintf("y[it_*y_size +73]=%f \n",yy);*/
+ y[k+y_decal*y_size] += slowc*(yy - y[k+y_decal*y_size]);
+//#ifdef PRINT_OUT_p
+ mexPrintf("=%f\n", y[k+y_decal*y_size]);
+//#endif
+ nop++;
+ }
+ //mexPrintf("nb_prologue_push=%d et nb_prologue_pop=%d\n",nb_prologue_push,nb_prologue_pop);
+ for(i = nb_first_table_y - 1;i >= 0;i--)
+ {
+ k = first_table_y[i].index;
+ yy = 0;
+//#ifdef PRINT_OUT_p
+ 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_p
+ if(j > 0)
+ mexPrintf("u[%d](%f)*y[%d](%f)+", first_table_y[i].u_index[j], uu, first_table_y[i].y_index[j], y[first_table_y[i].y_index[j]]);
+ else
+ mexPrintf("u[%d](%f)*y[%d](%f)", first_table_y[i].u_index[j], uu, first_table_y[i].y_index[j], 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);
+ filename+=" stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
+ /*if(((k-73) % y_size)==0)
+ mexPrintf("y[it_*y_size +73]=%f \n",yy);*/
+ y[k] += slowc*(yy - y[k]);
+//#ifdef PRINT_OUT_p
+ mexPrintf("=%f\n", y[k]);
+//#endif
+ nop++;
+ }
+ pctimer_t t2 = pctimer();
+ if(nb_first_table_u > 0)
+ {
+ mexPrintf("(**%f milliseconds u_count : %d nop : %d **)\n", 1000*(t2 - t1), u_count, nop);
+ mexEvalString("drawnow;");
+ }
+ /*if((nb_last_table_u>0)&&(it_>y_kmin))
+ mexErrMsgTxt("Exit from Dynare");*/
+}
+