From d90393f9df225b7b1430f23feab76d7fba7175ee Mon Sep 17 00:00:00 2001 From: ferhat Date: Tue, 20 Nov 2007 23:27:30 +0000 Subject: [PATCH] Adding the simulate files. git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1448 ac1d8469-bf42-47a9-8791-bf33cf982152 --- parser.src/simulate/Include/Interpreter.hh | 72 + parser.src/simulate/Include/Mem_Mngr.hh | 48 + parser.src/simulate/Include/SparseMatrix.hh | 149 + parser.src/simulate/Interpreter.cc | 1188 +++++++ parser.src/simulate/Makefile | 94 +- parser.src/simulate/Mem_Mngr.cc | 260 ++ parser.src/simulate/SparseMatrix.cc | 3468 +++++++++++++++++++ 7 files changed, 5236 insertions(+), 43 deletions(-) create mode 100644 parser.src/simulate/Include/Interpreter.hh create mode 100644 parser.src/simulate/Include/Mem_Mngr.hh create mode 100644 parser.src/simulate/Include/SparseMatrix.hh create mode 100644 parser.src/simulate/Interpreter.cc create mode 100644 parser.src/simulate/Mem_Mngr.cc create mode 100644 parser.src/simulate/SparseMatrix.cc diff --git a/parser.src/simulate/Include/Interpreter.hh b/parser.src/simulate/Include/Interpreter.hh new file mode 100644 index 000000000..7cfd39880 --- /dev/null +++ b/parser.src/simulate/Include/Interpreter.hh @@ -0,0 +1,72 @@ +#ifndef INTERPRETER_HH_INCLUDED +#define INTERPRETER_HH_INCLUDED + +#include +#include +#include +#include +#include "CodeInterpreter.hh" +#include "SymbolTableTypes.hh" +#include "ExprNode.hh" +//#include "Mem_Mngr.hh" +#include "SparseMatrix.hh" +#include "LinBCG.hh" +#include "mex.h" + +//#define DEBUGC + +using namespace std; + +typedef struct Block_contain_type +{ + int Equation, Variable, Own_Derivative; +}; + +typedef struct Block_type +{ + int begin, end, size, type; +}; + +#define pow_ pow +#define get_code_int *((int*)Code); Code+=sizeof(int); +#define get_code_double *((double*)Code); Code+=sizeof(double); +#define get_code_pdouble (double*)Code; Code+=sizeof(double); +#define get_code_bool *((bool*)(Code++)) +#define get_code_char *((char*)(Code++)) +#define get_code_pos (int)Code +#define get_code_pointer Code +#define set_code_pointer(pos) Code=pos + + +class Interpreter : SparseMatrix +{ + protected : + double pow1(double a, double b); + void compute_block_time(); + void simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination); + double *T; + vector Block_Contain; + vector Block; + char *Code; + stack Stack; + int Block_Count, Per_u_, Per_y_; + int it_, nb_row_x, nb_row_xd, maxit_, size_of_direction; + double *g1, *r; + //double max_res, res1, res2, slowc, markowitz_c; + double solve_tolf; + bool GaussSeidel; + double *x, *params; + //double *y, *ya, *x, *direction; + map ,int>, int> IM_i; + string filename; + public : + //ReadBinFile read_bin_file; + + Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *direction_arg, int y_size_arg, int nb_row_x_arg, + int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg, + double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg); + void compute_blocks(string file_name, string bin_basename); +}; + + +#endif diff --git a/parser.src/simulate/Include/Mem_Mngr.hh b/parser.src/simulate/Include/Mem_Mngr.hh new file mode 100644 index 000000000..bb14d2485 --- /dev/null +++ b/parser.src/simulate/Include/Mem_Mngr.hh @@ -0,0 +1,48 @@ +#ifndef MEM_MNGR_HH_INCLUDED +#define MEM_MNGR_HH_INCLUDED + +#include +#include +#include "mex.h" +using namespace std; + +typedef struct NonZeroElem + { + int u_index; + int r_index, c_index, lag_index; + NonZeroElem *NZE_R_N, *NZE_C_N; + }; + +typedef vector v_NonZeroElem; + +class Mem_Mngr +{ +public: + void write_swp_f(int *save_op_all,long int *nop_all); + bool read_swp_f(int **save_op_all,long int *nop_all); + void close_swp_f(); + void Print_heap(); + void init_Mem(); + void mxFree_NZE(void* pos); + NonZeroElem* mxMalloc_NZE(); + void init_CHUNK_BLCK_SIZE(int u_count); + void Free_All(); + Mem_Mngr(); + void fixe_file_name(string filename_arg); + int* malloc_std(long int nop); + int* realloc_std(int* save_op_o, long int &nopa); + void chk_avail_mem(int **save_op_all,long int *nop_all,long int *nopa_all,int add, int t); + bool swp_f; + //bool verbose; +private: + v_NonZeroElem Chunk_Stack; + int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK; + int CHUNK_heap_pos/*, CHUNK_heap_max_size*/; + NonZeroElem** NZE_Mem_add; + NonZeroElem* NZE_Mem; + int swp_f_b; + fstream SaveCode_swp; + string filename; +}; + +#endif diff --git a/parser.src/simulate/Include/SparseMatrix.hh b/parser.src/simulate/Include/SparseMatrix.hh new file mode 100644 index 000000000..e4632eeec --- /dev/null +++ b/parser.src/simulate/Include/SparseMatrix.hh @@ -0,0 +1,149 @@ +#ifndef SPARSEMATRIX_HH_INCLUDED +#define SPARSEMATRIX_HH_INCLUDED + +#include +#include +#include +#include +#include +#include "Mem_Mngr.hh" + +#define NEW_ALLOC +#define MARKOVITZ +//#define MEMORY_LEAKS + +using namespace std; + +typedef struct t_save_op_s +{ + short int lag, operat; + int first, second; +}; + +const int IFLD =0; +const int IFDIV =1; +const int IFLESS=2; +const int IFSUB =3; +const int IFLDZ =4; +const int IFMUL =5; +const int IFSTP =6; +const int IFADD =7; +const double eps=1e-7; +const double very_big=1e24; +const int alt_symbolic_count_max=1; + +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; + }; + + +class SparseMatrix + { + public: + SparseMatrix(); + int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter); + void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter); + void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1); + //void initialize(int periods_arg, int nb_endo_arg, int y_kmin_arg, int y_kmax_arg, int y_size_arg, int u_count_arg, int u_count_init_arg, double *u_arg, double *y_arg, double *ya_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg, double res1_arg, double res2_arg, double max_res_arg); + void Read_SparseMatrix(string file_name, int Size, int periods, int y_kmin, int y_kmax); + void Read_file(string file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax, int &nb_endo, int &u_count, int &u_count_init, double* u); + + private: + void Init(int periods, int y_kmin, int y_kmax, int Size, map ,int>, int> IM); + void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map ,int>, int> IM); + void End(int Size); + bool compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size +#ifdef PROFILER + , long int *ndiv, long int *nsub +#endif + ); + void Insert(int r, int c, int u_index, int lag_index); + void Delete(int r,int c, int Size, int *b); + int At_Row(int r, NonZeroElem **first); + int At_Pos(int r, int c, NonZeroElem **first); + int At_Col(int c, NonZeroElem **first); + int At_Col(int c, int lag, NonZeroElem **first); + int NRow(int r); + int NCol(int c); + int Union_Row(int row1, int row2); + void Print(int Size,int *b); + int Get_u(); + void Delete_u(int pos); + void Clear_u(); + void Print_u(); + int complete(int beg_t, int Size, int periods, int *b); + double bksub( int tbreak, int last_period, int Size, double slowc_l +#ifdef PROFILER + , long int *nmul +#endif + ); + void run_triangular(int nop_all,int *op_all); + void run_it(int nop_all,int *op_all); + void run_u_period1(int periods); + void close_swp_file(); + + 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, int y_kmin, int y_kmax, int u_size); + 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 y_kmin, int y_kmax, int u_size, int y_size); + void close_SaveCode(); + + stack Stack; + 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 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; + fstream SaveCode; + string filename; + int max_u, min_u; + clock_t time00; + + Mem_Mngr mem_mngr; + vector u_liste; + int *NbNZRow, *NbNZCol; + NonZeroElem **FNZE_R, **FNZE_C; + int nb_endo, u_count_init; + + int *pivot, *pivotk; + double *pivotv, *pivotva; + int *b; + bool *line_done; + bool symbolic, alt_symbolic; + int alt_symbolic_count; + int *g_save_op; + int first_count_loop; + int g_nop_all; + int u_count_alloc, u_count_alloc_save; + double markowitz_c_s; + double res1a; + long int nop_all, /*nopa_all,*/ nop1, nop2; + map ,int>, int> IM_i; +protected: + double *u, *y, *ya; + double res1, res2, max_res; + double slowc, slowc_save, markowitz_c; + int y_kmin, y_kmax, y_size, periods, y_decal; +//public: + int *index_vara, *index_equa; + int u_count, tbreak_g; + int iter; + double *direction; + + }; + + + + + +#endif diff --git a/parser.src/simulate/Interpreter.cc b/parser.src/simulate/Interpreter.cc new file mode 100644 index 000000000..5487f320c --- /dev/null +++ b/parser.src/simulate/Interpreter.cc @@ -0,0 +1,1188 @@ +#include "Interpreter.hh" + +Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *direction_arg, int y_size_arg, + int nb_row_x_arg, int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, + int maxit_arg_, double solve_tolf_arg, int size_of_direction_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg, + string &filename_arg) +{ + params=params_arg; + y=y_arg; + ya=ya_arg; + x=x_arg; + direction=direction_arg; + y_size=y_size_arg; + nb_row_x=nb_row_x_arg; + nb_row_xd=nb_row_xd_arg; + periods=periods_arg; + y_kmax=y_kmax_arg; + y_kmin=y_kmin_arg; + maxit_=maxit_arg_; + solve_tolf=solve_tolf_arg; + size_of_direction=size_of_direction_arg; + slowc=slowc_arg; + y_decal=y_decal_arg; + markowitz_c=markowitz_c_arg; + filename=filename_arg; + //GaussSeidel=true; +} + +double +Interpreter::pow1(double a, double b) +{ + double r=pow_(a,b); + if (isnan(r) || isinf(r)) + { + max_res=res1=res2=r; + return(r); + } + else + return r; +} + + +void +Interpreter::compute_block_time() /*throw(EvalException)*/ +{ + int var, lag, op; + double v1, v2; + char cc; + bool go_on=true; + double *ll; + while (go_on) + { + //mexPrintf("it_=%d",it_); + switch (cc=get_code_char) + { + case FLDV : + //load a variable in the processor +#ifdef DEBUGC + mexPrintf("FLDV"); + mexEvalString("drawnow;"); +#endif + switch (get_code_char) + { + case eParameter : + var=get_code_int +#ifdef DEBUGC + mexPrintf(" params[%d]=%f\n",var,params[var]); + mexEvalString("drawnow;"); +#endif + Stack.push(params[var]); + break; + case eEndogenous : + var=get_code_int + lag=get_code_int +#ifdef DEBUGC + mexPrintf(" y[%d]=%f\n",(it_+lag)*y_size+var,y[(it_+lag)*y_size+var]); + mexEvalString("drawnow;"); +#endif + Stack.push(y[(it_+lag)*y_size+var]); + break; + case eExogenous : + var=get_code_int + lag=get_code_int +#ifdef DEBUGC + mexPrintf(" x[%d]=%f\n",it_+lag+var*nb_row_x,x[it_+lag+var*nb_row_x]); + mexEvalString("drawnow;"); +#endif + Stack.push(x[it_+lag+var*nb_row_x]); + break; + case eExogenousDet : + 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;"); +#endif + Stack.push(x[it_+lag+var*nb_row_xd]); + break; + } + break; + case FLDT : + //load a temporary variable in the processor + var=get_code_int +#ifdef DEBUGC + mexPrintf("FLDT %d\n",var); + mexEvalString("drawnow;"); +#endif + Stack.push(T[var*(periods+y_kmin+y_kmax)+it_]); + break; + case FLDU : + //load u variable in the processor +#ifdef DEBUGC + mexPrintf("FLDU\n"); + mexEvalString("drawnow;"); +#endif + var=get_code_int + var+=Per_u_; + Stack.push(u[var]); + break; + case FLDR : + //load u variable in the processor +#ifdef DEBUGC + mexPrintf("FLDR\n"); + mexEvalString("drawnow;"); +#endif + var=get_code_int + Stack.push(r[var]); + break; + case FLDZ : + //load 0 in the processor +#ifdef DEBUGC + mexPrintf("FLDZ\n"); + mexEvalString("drawnow;"); +#endif + Stack.push(0); + break; + case FLDC : + //load a numerical constant in the processor + /*asm("fld\n\t" + "fstp %%st" : "=t" (ll) : "0" ((double)(*Code)));*/ + ll=get_code_pdouble +#ifdef DEBUGC + mexPrintf("FLDC %f\n",*ll); + mexEvalString("drawnow;"); +#endif + Stack.push(*ll); + break; + case FSTPV : + //load a variable in the processor +#ifdef DEBUGC + mexPrintf("FSTPV\n"); + mexEvalString("drawnow;"); +#endif + switch (get_code_char) + { + case eParameter : + var=get_code_int + params[var] = Stack.top(); + Stack.pop(); + break; + case eEndogenous : + 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()); + mexEvalString("drawnow;"); +#endif + y[(it_+lag)*y_size+var] = Stack.top(); +#ifdef DEBUGC + mexPrintf("%f\n",y[(it_+lag)*y_size+var]); + mexEvalString("drawnow;"); +#endif + Stack.pop(); + break; + case eExogenous : + var=get_code_int + lag=get_code_int + x[it_+lag+var*nb_row_x] = Stack.top(); + Stack.pop(); + break; + case eExogenousDet : + var=get_code_int + lag=get_code_int + x[it_+lag+var*nb_row_xd] = Stack.top(); + Stack.pop(); + break; + } + break; + case FSTPT : + //load a temporary variable in the processor + 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;"); +#endif + T[var*(periods+y_kmin+y_kmax)+it_] = Stack.top(); +#ifdef DEBUGC + mexPrintf("%f\n",T[var*(periods+y_kmin+y_kmax)+it_]); + mexEvalString("drawnow;"); +#endif + Stack.pop(); + break; + case FSTPU : + //load u variable in the processor + var=get_code_int + var+=Per_u_; +#ifdef DEBUGC + mexPrintf("FSTPU u[%d]",var); + mexEvalString("drawnow;"); +#endif + u[var] = Stack.top(); +#ifdef DEBUGC + mexPrintf("=%f\n",u[var]); + mexEvalString("drawnow;"); +#endif + Stack.pop(); + break; + case FSTPR : + //load u variable in the processor + var=get_code_int + r[var] = Stack.top(); +#ifdef DEBUGC + mexPrintf("FSTPR residual[%d]=%f\n",var,r[var]); + mexEvalString("drawnow;"); +#endif + Stack.pop(); + break; + case FSTPG : + //load u variable in the processor + var=get_code_int + g1[var] = Stack.top(); +#ifdef DEBUGC + mexPrintf("FSTPG g1[%d)=%f\n",var,g1[var]); + mexEvalString("drawnow;"); +#endif + Stack.pop(); + break; + case FBINARY : +#ifdef DEBUGC + mexPrintf("FBINARY\n"); + mexEvalString("drawnow;"); +#endif + op=get_code_int + v2=Stack.top(); + Stack.pop(); + v1=Stack.top(); + Stack.pop(); + switch (op) + { + case oPlus: + Stack.push(v1 + v2); +#ifdef DEBUGC + mexPrintf("+\n"); + mexEvalString("drawnow;"); +#endif + break; + case oMinus: + Stack.push(v1 - v2); +#ifdef DEBUGC + mexPrintf("-\n"); + mexEvalString("drawnow;"); +#endif + break; + case oTimes: + Stack.push(v1 * v2); +#ifdef DEBUGC + mexPrintf("*\n"); + mexEvalString("drawnow;"); +#endif + break; + case oDivide: + Stack.push(v1 / v2); +#ifdef DEBUGC + mexPrintf("/\n"); + mexEvalString("drawnow;"); +#endif + break; + case oLess: + Stack.push(double(v1v2)); +#ifdef DEBUGC + mexPrintf("%f > %f\n",v1,v2); + mexEvalString("drawnow;"); +#endif + break; + case oLessEqual: + Stack.push(double(v1<=v2)); +#ifdef DEBUGC + mexPrintf("%f <= %f\n",v1,v2); + mexEvalString("drawnow;"); +#endif + break; + case oGreaterEqual: + Stack.push(double(v1>=v2)); +#ifdef DEBUGC + mexPrintf("%f >= %f\n",v1,v2); + mexEvalString("drawnow;"); +#endif + break; + case oEqualEqual: + Stack.push(double(v1==v2)); +#ifdef DEBUGC + mexPrintf("%f == %f\n",v1,v2); + mexEvalString("drawnow;"); +#endif + break; + case oDifferent: + Stack.push(double(v1!=v2)); +#ifdef DEBUGC + mexPrintf("%f > %f\n",v1,v2); + mexEvalString("drawnow;"); +#endif + break; + case oPower: + Stack.push(pow1(v1, v2)); +#ifdef DEBUGC + mexPrintf("pow(%f, %f)\n",v1,v2); + mexEvalString("drawnow;"); +#endif + break; + case oMax: + Stack.push(max(v1, v2)); +#ifdef DEBUGC + mexPrintf("max(%f, %f)\n",v1,v2); + mexEvalString("drawnow;"); +#endif + break; + case oMin: + Stack.push(min(v1, v2)); +#ifdef DEBUGC + mexPrintf("min(%f, %f)\n",v1,v2); + mexEvalString("drawnow;"); +#endif + break; + case oEqual: + default: + /*throw EvalException();*/ + ; + } + break; + case FUNARY : +#ifdef DEBUGC + mexPrintf("FUNARY\n"); + mexEvalString("drawnow;"); +#endif + op=get_code_int + v1=Stack.top(); + Stack.pop(); + switch (op) + { + case oUminus: + Stack.push(-v1); +#ifdef DEBUGC + mexPrintf("-\n"); + mexEvalString("drawnow;"); +#endif + break; + case oExp: + Stack.push(exp(v1)); +#ifdef DEBUGC + mexPrintf("exp\n"); + mexEvalString("drawnow;"); +#endif + break; + case oLog: + Stack.push(log(v1)); +#ifdef DEBUGC + mexPrintf("log\n"); + mexEvalString("drawnow;"); +#endif + break; + case oLog10: + Stack.push(log10(v1)); +#ifdef DEBUGC + mexPrintf("log10\n"); + mexEvalString("drawnow;"); +#endif + break; + case oCos: + Stack.push(cos(v1)); +#ifdef DEBUGC + mexPrintf("cos\n"); + mexEvalString("drawnow;"); +#endif + break; + case oSin: + Stack.push(sin(v1)); +#ifdef DEBUGC + mexPrintf("sin\n"); + mexEvalString("drawnow;"); +#endif + break; + case oTan: + Stack.push(tan(v1)); +#ifdef DEBUGC + mexPrintf("tan\n"); + mexEvalString("drawnow;"); +#endif + break; + case oAcos: + Stack.push(acos(v1)); +#ifdef DEBUGC + mexPrintf("acos\n"); + mexEvalString("drawnow;"); +#endif + break; + case oAsin: + Stack.push(asin(v1)); +#ifdef DEBUGC + mexPrintf("asin\n"); + mexEvalString("drawnow;"); +#endif + break; + case oAtan: + Stack.push(atan(v1)); +#ifdef DEBUGC + mexPrintf("atan\n"); + mexEvalString("drawnow;"); +#endif + break; + case oCosh: + Stack.push(cosh(v1)); +#ifdef DEBUGC + mexPrintf("cosh\n"); + mexEvalString("drawnow;"); +#endif + break; + case oSinh: + Stack.push(sinh(v1)); +#ifdef DEBUGC + mexPrintf("sinh\n"); + mexEvalString("drawnow;"); +#endif + break; + case oTanh: + Stack.push(tanh(v1)); +#ifdef DEBUGC + mexPrintf("tanh\n"); + mexEvalString("drawnow;"); +#endif + break; + case oAcosh: + Stack.push(acosh(v1)); +#ifdef DEBUGC + mexPrintf("acosh\n"); + mexEvalString("drawnow;"); +#endif + break; + case oAsinh: + Stack.push(asinh(v1)); +#ifdef DEBUGC + mexPrintf("asinh\n"); + mexEvalString("drawnow;"); +#endif + break; + case oAtanh: + Stack.push(atanh(v1)); +#ifdef DEBUGC + mexPrintf("atanh\n"); + mexEvalString("drawnow;"); +#endif + break; + case oSqrt: + Stack.push(sqrt(v1)); +#ifdef DEBUGC + mexPrintf("sqrt\n"); + mexEvalString("drawnow;"); +#endif + break; + default: + /*throw EvalException();*/ + ; + } + break; + case FCUML : +#ifdef DEBUGC + mexPrintf("FCUML\n"); + mexEvalString("drawnow;"); +#endif + v1=Stack.top(); + Stack.pop(); + v2=Stack.top(); + Stack.pop(); + Stack.push(v1+v2); + break; + case FENDBLOCK : + //it's the block end +#ifdef DEBUGC + mexPrintf("FENDBLOCK\n"); + mexEvalString("drawnow;"); +#endif + //Block[Block_Count].end=get_code_pos; + go_on=false; + break; + case FENDEQU : +#ifdef DEBUGC + mexPrintf("FENDEQU\n"); + mexEvalString("drawnow;"); +#endif + /*if (GaussSeidel) + return;*/ + break; + case FOK : +#ifdef DEBUGC + mexPrintf("FOK\n"); + mexEvalString("drawnow;"); +#endif + op=get_code_int +#ifdef DEBUGC + mexPrintf("var=%d\n",op); +#endif + if (Stack.size()>0) + { + mexPrintf("error: Stack not empty!\n"); + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt("End of simulate"); + } + break; + default : + mexPrintf("Unknow opcode %d!! FENDEQU=%d\n",cc,FENDEQU); + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt("End of simulate"); + break; + } + } +} + +void +Interpreter::simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination) +{ + /*mexPrintf("simulate_a_block\n"); + mexEvalString("drawnow;");*/ + + char *begining; + int i; + bool is_linear, cvg; + int max_lag_plus_max_lead_plus_1; + int symbol_table_endo_nbr; + int Block_List_Max_Lag; + int Block_List_Max_Lead; + int giter; + int u_count_int; + double *y_save; + LinBCG linbcg; + Mat_DP a; + Vec_INT indx; + //SparseMatrix sparse_matrix; + + int nb_endo, u_count_init; + + + //mexPrintf("simulate_a_block\n"); + //mexEvalString("drawnow;"); + //mexPrintf("%d\n",debile); + + //GaussSeidel=false; + //slowc_save=slowc/2; + //mexPrintf("simulate_a_block size=%d type=%d\n",size,type); + switch (type) + { + case EVALUATE_FOREWARD : + case EVALUATE_FOREWARD_R : +#ifdef DEBUGC + mexPrintf("EVALUATE_FOREWARD\n"); +#endif + begining=get_code_pointer; + for (it_=y_kmin;it_y_kmin;it_--) + { + set_code_pointer(begining); + Per_y_=it_*y_size; + compute_block_time(); +#ifdef PRINT_OUT + for (j = 0; jmaxit_))) + { + set_code_pointer(begining); + Per_y_=it_*y_size; + compute_block_time(); + y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0]; + //mexPrintf("y[%d] += -r[0] (%f) / g1[0] (%f) = %f\n",Per_y_+Block_Contain[0].Variable,r[0],g1[0],y[Per_y_+Block_Contain[0].Variable]); + double rr; + rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]); + cvg=((rr*rr)y_kmin;it_--) + { + cvg=false; + iter=0; + Per_y_=it_*y_size; + while (!(cvg||(iter>maxit_))) + { + set_code_pointer(begining); + Per_y_=it_*y_size; + compute_block_time(); + y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0]; + double rr; + rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]); + cvg=((rr*rr)maxit_))) + { + res2=0; + res1=0; + max_res=0; + for (it_=y_kmin;it_maxit_))) + { + set_code_pointer(begining); + compute_block_time(); + Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + res2=0; + res1=0; + max_res=0; + for (i=0; iy_kmin;it_--) + { + cvg=false; + iter=0; + Per_y_=it_*y_size; + while (!(cvg||(iter>maxit_))) + { + set_code_pointer(begining); + compute_block_time(); + Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + res2=0; + res1=0; + max_res=0; + for (i=0; iy_kmin;it_--) + { + set_code_pointer(begining); + Per_y_=it_*y_size; + compute_block_time(); + Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + } + } + mxFree(g1); + mxFree(r); + mxFree(u); + break; + case SOLVE_TWO_BOUNDARIES_COMPLETE: +#ifdef DEBUGC + mexPrintf("SOLVE_TWO_BOUNDARIES_COMPLETE\n"); + mexEvalString("drawnow;"); +#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 +#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 +#ifdef DEBUGC + mexPrintf("symbol_table_endo_nbr=%d\n",symbol_table_endo_nbr); + mexEvalString("drawnow;"); +#endif + 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 +#ifdef DEBUGC + mexPrintf("Block_List_Max_Lead=%d\n",Block_List_Max_Lead); + mexEvalString("drawnow;"); +#endif + u_count_int=get_code_int +#ifdef DEBUGC + mexPrintf("u_count_int=%d\n",u_count_int); + mexPrintf("periods=%d\n",periods); + mexEvalString("drawnow;"); +#endif + + //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); + + fixe_u(&u, u_count_int, max_lag_plus_max_lead_plus_1); +#ifdef DEBUGC + mexPrintf("u=%x\n",u); +#endif + Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax); + + //mexPrintf("aft reading_sparse_matrix\n"); + //mexEvalString("drawnow;"); + u_count=u_count_int*(periods+y_kmax+y_kmin); + g1=(double*)mxMalloc(size*size*sizeof(double)); + r=(double*)mxMalloc(size*sizeof(double)); + y_save=(double*)mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); +#ifdef DEBUGC + mexPrintf("u_count=%d\n",u_count); + mexEvalString("drawnow;"); +#endif + begining=get_code_pointer; + if(!Gaussian_Elimination) + { + it_=y_kmin; + Per_u_=0; + Per_y_=it_*y_size; + set_code_pointer(begining); + compute_block_time(); + linbcg.Initialize(filename, res1, res2, max_res, slowc, ya, direction, iter); + linbcg.Preconditioner(periods, y_kmin, y_kmax, size, IM_i, index_vara, index_equa, y_size, y, true, 0, a, indx); + } + //GaussSeidel=false; + giter=0; + iter=0; + //mexPrintf("2 boudaries problem\n"); + //mexEvalString("drawnow;"); + //mexPrintf("GaussSeidel=%d\n",GaussSeidel); + if (!is_linear) + { + //double res1a=0; + cvg=false; + int u_count_saved=u_count; + while (!(cvg||(iter>maxit_))) + { + res2=0; + res1=0; + max_res=0; + memcpy(y_save, y, y_size*sizeof(double)*(periods+y_kmax+y_kmin)); + for (it_=y_kmin;it_solve_tolf)*/ + rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); + /*else + rr=r[i];*/ + if (max_res(Code), Code_Size); + CompiledCode.close(); + char *Init_Code=Code; + + //The big loop on intructions + Block_Count=-1; + bool go_on=true; + while (go_on) + { +#ifdef DEBUGC + mexPrintf("pos=%d\n",int(get_code_pos)-int(Init_Code)); + mexEvalString("drawnow;"); +#endif + char code=get_code_char; +#ifdef DEBUGC + int icode=(int)code; + mexPrintf("code=%d\n",icode); + mexEvalString("drawnow;"); +#endif + switch (code) + { + case FBEGINBLOCK : + //it's a new block + Block_Count++; + Block_type lBlock; + Block.clear(); + Block_Contain.clear(); + Block_contain_type lBlock_Contain; +#ifdef DEBUGC + mexPrintf("FBEGINBLOCK\n"); + mexEvalString("drawnow;"); +#endif + lBlock.begin=get_code_pos-int(Init_Code); +#ifdef DEBUGC + mexPrintf("Block[Block_Count].begin=%d\n",lBlock.begin); + mexEvalString("drawnow;"); +#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 + mexPrintf("Block[Block_Count].type=%d\n",lBlock.type); + mexEvalString("drawnow;"); +#endif + Block.push_back(lBlock); + for (int i=0;i=0)) + break; + } + /*if(verbose) + { + mexPrintf("push_back()\n"); + mexEvalString("drawnow;"); + }*/ + Chunk_Stack.push_back((NonZeroElem*)pos); + /*if(verbose) + { + mexPrintf("End\n"); + mexEvalString("drawnow;"); + }*/ +} + + +void +Mem_Mngr::write_swp_f(int *save_op_all,long int *nop_all) +{ + swp_f=true; + swp_f_b++; + mexPrintf("writing the block %d with size=%d\n",swp_f_b,*nop_all); + if (!SaveCode_swp.is_open()) + { + mexPrintf("open the swp file for writing\n"); +#ifdef PRINT_OUT + mexPrintf("file opened\n"); +#endif + SaveCode_swp.open((filename + ".swp").c_str(), std::ios::out | std::ios::binary); + if (!SaveCode_swp.is_open()) + { + mexPrintf("Error : Can't open file \"%s\" for writing\n", (filename + ".swp").c_str()); + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt("Exit from Dynare"); + } +#ifdef PRINT_OUT + mexPrintf("done\n"); +#endif + } + SaveCode_swp.write(reinterpret_cast(nop_all), sizeof(*nop_all)); + SaveCode_swp.write(reinterpret_cast(save_op_all), (*nop_all)*sizeof(int)); + (*nop_all)=0; +} + +bool +Mem_Mngr::read_swp_f(int **save_op_all,long int *nop_all) +{ + int j; + swp_f=true; + if (!SaveCode_swp.is_open()) + { +#ifdef PRINT_OUT + mexPrintf("file opened\n"); +#endif + mexPrintf("open the file %s\n",(filename + ".swp").c_str()); + SaveCode_swp.open((filename + ".swp").c_str(), std::ios::in | std::ios::binary); + j=SaveCode_swp.is_open(); + mexPrintf("is_open()=%d\n",j); + + if (!SaveCode_swp.is_open()) + { + mexPrintf("Error : Can't open file \"%s\" for reading\n", (filename + ".swp").c_str()); + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt("Exit from Dynare"); + } +#ifdef PRINT_OUT + mexPrintf("done\n"); +#endif + SaveCode_swp.seekg(0); + } + + j=SaveCode_swp.tellg(); + SaveCode_swp.read(reinterpret_cast(nop_all), sizeof(*nop_all)); + (*save_op_all)=(int*)mxMalloc((*nop_all)*sizeof(int)); + SaveCode_swp.read(reinterpret_cast(*save_op_all), (*nop_all)*sizeof(int)); + return(SaveCode_swp.good()); +} + + +void +Mem_Mngr::close_swp_f() +{ + if (SaveCode_swp.is_open()) + { + SaveCode_swp.close(); + mexPrintf("close the swp file\n"); + } +} + +int* +Mem_Mngr::malloc_std(long int nop) +{ + return((int*)malloc(nop*sizeof(int))); +} + +int* +Mem_Mngr::realloc_std(int* save_op_o, long int &nopa) +{ + int *save_op=(int*)realloc(save_op_o,nopa*sizeof(int)); + if (!save_op) + { + int nopag=int(nopa/3); + nopa=nopa-nopag; + while (!save_op && nopag>0) + { + nopag=int(nopag*0.66); + save_op=(int*)realloc(save_op_o,nopa*sizeof(int)); + } + if (!save_op) + { + mexPrintf("Memory exhausted\n"); + mexEvalString("drawnow;"); + mexEvalString("st=fclose('all');clear all;"); + filename+=" stopped"; + mexErrMsgTxt(filename.c_str()); + } + } + return(save_op); +} + +void +Mem_Mngr::chk_avail_mem(int **save_op_all,long int *nop_all,long int *nopa_all,int add, int t) + { + mexPrintf("Error: out of save_op_all[%d] nopa_all=%d t=%d\n",(*nop_all)+add,(*nopa_all),t); + int tmp_nopa_all=int(1.5*(*nopa_all)); + int *tmp_i; + if (tmp_nopa_all*sizeof(int)<1024*1024) + { + mexPrintf("allocate %d bites save_op_all=%x\n",tmp_nopa_all*sizeof(int),*save_op_all); + tmp_i=(int*)mxRealloc(*save_op_all,tmp_nopa_all*sizeof(int)); + mexPrintf("tmp_i="); + mexPrintf("%x\n",tmp_i); + } + else + tmp_i=NULL; + if (!tmp_i) + { + write_swp_f((*save_op_all),nop_all); + } + else + { + mexPrintf("allocated\n"); + (*save_op_all)=tmp_i; + (*nopa_all)=tmp_nopa_all; + } + mexPrintf("end of chk\n"); + } + +void +Mem_Mngr::Free_All() +{ + int i; + for (int i=0;ic_index==first2->c_index) + { + nb_elem++; + i1++; + i2++; + first1=first1->NZE_R_N; + first2=first2->NZE_R_N; + } + else if (first1->c_indexc_index) + { + nb_elem++; + i1++; + first1=first1->NZE_R_N; + } + else + { + nb_elem++; + i2++; + first2=first2->NZE_R_N; + } + } + return nb_elem; +} + +int +SparseMatrix::At_Pos(int r, int c, NonZeroElem **first) +{ + (*first)=FNZE_R[r]; + while ((*first)->c_index!=c /*&& (*first)->NZE_R_N*/) + { +#ifdef PRINT_OUT + mexPrintf("looking not CRS [%d, %d]\n",(*first)->r_index,(*first)->c_index); +#endif + (*first)=(*first)->NZE_R_N; + } + /*if ((*first)->c_index!=c) + mexPrintf("----------------------- cannot find M[%d, %d]\n",r,c);*/ + return NbNZRow[r]; +} + + +int SparseMatrix::At_Col(int c, NonZeroElem **first) +{ + (*first)=FNZE_C[c]; + return NbNZCol[c]; +} + +int SparseMatrix::At_Col(int c, int lag, NonZeroElem **first) +{ + (*first)=FNZE_C[c]; + int i=0; + while ((*first)->lag_index!=lag && (*first)) + { +#ifdef PRINT_OUT + mexPrintf("first->lag_index(%d) != %d\n",(*first)->lag_index,lag); +#endif + (*first)=(*first)->NZE_C_N; + } + if ((*first)) + { +#ifdef PRINT_OUT + mexPrintf("first=%x\n",(*first)); +#endif + NonZeroElem* firsta=(*first); + if (!firsta->NZE_C_N) + i++; + else + { + while (firsta->lag_index==lag && firsta->NZE_C_N) + { +#ifdef PRINT_OUT + mexPrintf("firsta->lag_index(%d) == %d, eq=%d, var=%d\n",firsta->lag_index,lag, firsta->r_index, firsta->c_index); +#endif + firsta=firsta->NZE_C_N; + i++; + } + if (firsta->lag_index==lag) i++; + } + } +#ifdef PRINT_OUT + mexPrintf("i=%d\n",i); +#endif + return i; +} + +#ifdef PROFILER +double tdelete1=0, tdelete2=0, tdelete21=0, tdelete22=0, tdelete221=0, tdelete222=0, tcompare=0; +#endif + +void SparseMatrix::Delete(int r,int c, int Size, int *b) +{ + NonZeroElem *first=FNZE_R[r], *firsta=NULL; +#ifdef PROFILER + clock_t td0, td1, td2; + td0=clock(); +#endif + while (first->c_index!=c /*&& first->NZE_R_N*/) + { +#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; + } +#ifdef PRINT_OUT + mexPrintf("CRS [%d, %d]=c(%d)\n",first->r_index,first->c_index,c); + mexEvalString("drawnow;"); +#endif + if (firsta!=NULL) + firsta->NZE_R_N=first->NZE_R_N; + if (first==FNZE_R[r]) + FNZE_R[r]=first->NZE_R_N; + NbNZRow[r]--; +#ifdef PROFILER + tdelete1+=clock()-td0; + td0=clock(); + td1=clock(); +#endif + first=FNZE_C[c]; + firsta=NULL; + while (first->r_index!=r /*&& first->NZE_C_N*/) + { +#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; + } +#ifdef PRINT_OUT + mexPrintf("CSS [%d, %d]=r(%d)\n",first->r_index,first->c_index,r); + mexEvalString("drawnow;"); +#endif +#ifdef PROFILER + tdelete21+=clock()-td1; + td1=clock(); +#endif + if (firsta!=NULL) + firsta->NZE_C_N=first->NZE_C_N; + if (first==FNZE_C[c]) + FNZE_C[c]=first->NZE_C_N; +#ifdef PROFILER + td2=clock(); +#endif + u_liste.push_back(first->u_index); +#ifdef PROFILER + tdelete221+=clock()-td2; + td2=clock(); +#endif +#ifdef NEW_ALLOC + mem_mngr.mxFree_NZE(first); +#else + mxFree(first); +#endif + NbNZCol[c]--; +#ifdef PROFILER + tdelete222+=clock()-td2; +#endif +#ifdef PROFILER + tdelete22+=clock()-td1; + tdelete2+=clock()-td0; +#endif +} + + +void SparseMatrix::Print(int Size, int *b) +{ + int a,i,j,k,l; + mexPrintf(" "); + for (k=0;kc_index-a);l++) + mexPrintf(" "); + mexPrintf("%-2d ",first->u_index); + a=first->c_index+1; + first=first->NZE_R_N; + } + for (k=a;kc_index-a);l++) + mexPrintf(" "); + mexPrintf("%8.4f",double(u[first->u_index])); + a=first->c_index+1; + first=first->NZE_R_N; + } + for (k=a;kc_index=%d, first->NZE_R_N=%x\n",first->c_index, first->NZE_R_N); +#endif + while (first->c_indexNZE_R_N) + { + firsta=first; +#ifdef PRINT_OUT + mexPrintf("drop first->c_index=%d c=%d\n",first->c_index,c); +#endif + first=first->NZE_R_N; + } + /*if (first->c_index!=c) + {*/ +#ifdef PRINT_OUT + mexPrintf("retain first->c_index=%d c=%d\n",first->c_index,c); +#endif +#ifdef NEW_ALLOC + firstn=mem_mngr.mxMalloc_NZE(); +#else + firstn=(NonZeroElem*)mxMalloc(sizeof(NonZeroElem)); +#endif + firstn->u_index=u_index; + firstn->r_index=r; + firstn->c_index=c; + firstn->lag_index=lag_index; + if (first->c_index>c) + { + if (first==FNZE_R[r]) + FNZE_R[r]=firstn; + if (firsta!=NULL) + firsta->NZE_R_N=firstn; + firstn->NZE_R_N=first; + } + else /*first.c_indexNZE_R_N=firstn; + firstn->NZE_R_N=NULL; + } + NbNZRow[r]++; + + /*} + else + mexPrintf("Error (in Insert): in CRS element r=%, c=%d already exists\n",r,c);*/ + } + else + { +#ifdef NEW_ALLOC + firstn=mem_mngr.mxMalloc_NZE(); +#else + firstn=(NonZeroElem*)mxMalloc(sizeof(NonZeroElem)); +#endif + firstn->u_index=u_index; + firstn->r_index=r; + firstn->c_index=c; + firstn->lag_index=lag_index; + FNZE_R[r]=firstn; + firstn->NZE_R_N=first; + NbNZRow[r]++; + } + first=FNZE_C[c]; + firsta=NULL; + while (first->r_indexNZE_C_N) + { + firsta=first; + first=first->NZE_C_N; + } + /*if (first->r_index!=r) + {*/ + if (first->r_index>r) + { + if (first==FNZE_C[c]) + FNZE_C[c]=firstn; + if (firsta!=NULL) + firsta->NZE_C_N=firstn; + firstn->NZE_C_N=first; + } + else /*first.r_indexNZE_C_N=firstn; + firstn->NZE_C_N=NULL; + } + NbNZCol[c]++; + /*} + else + mexPrintf("Error (in Insert): in CCS element r=%, c=%d already exists\n",r,c);*/ +} + +void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, int y_kmin, int y_kmax) +{ + int i,j,eq,var,lag; + /*mexPrintf("in Read_SparseMatrix\n"); + mexPrintf("file_name=%s SaveCode.is_open()=%d\n",file_name.c_str(), SaveCode.is_open()); + mexEvalString("drawnow;");*/ + filename=file_name; + mem_mngr.fixe_file_name(file_name); + 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 + } + IM_i.clear(); + /*mexPrintf("u_count_init=%d\n",u_count_init); + mexEvalString("drawnow;");*/ + for (i=0;i(&eq), sizeof(eq)); + SaveCode.read(reinterpret_cast(&var), sizeof(var)); + SaveCode.read(reinterpret_cast(&lag), sizeof(lag)); + SaveCode.read(reinterpret_cast(&j), sizeof(j)); + /*mexPrintf("0i=%d eq=%d var=%d lag=%d j=%d\n",i,eq,var,lag,j ); + mexEvalString("drawnow;");*/ + IM_i[std::make_pair(std::make_pair(eq, var), lag)] = j; + /*mexPrintf("1i=%d\n",i); + mexEvalString("drawnow;");*/ + } +#ifdef MEM_ALLOC_CHK + mexPrintf("index_vara=(int*)mxMalloc(%d*sizeof(int))\n",Size*(periods+y_kmin+y_kmax)); +#endif + /*mexPrintf("Size=%d periods=%d y_kmin=%d y_kmax=%d\n",Size,periods, y_kmin, y_kmax); + mexEvalString("drawnow;");*/ + index_vara=(int*)mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int)); +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + for (j=0;j(&index_vara[j]), sizeof(*index_vara)); + //mexPrintf("index_vara[%d]=%d\n",j,index_vara[j]); + } + for (i=1;i(&index_equa[j]), sizeof(*index_equa)); + //mexPrintf("index_equa[%d]=%d\n",j,index_equa[j]); + } +} + + +void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map ,int>, int> IM) +{ + int t,i, eq, var, lag; + double tmp_b=0.0; + std::map ,int>, int>::iterator it4; + NonZeroElem* first; + //mexPrintf("periods=%d, y_kmin=%d, y_kmax=%d, Size=%d, IM.size()=%d\n",periods, y_kmin, y_kmax, Size, IM.size()); +#ifdef MEM_ALLOC_CHK + mexPrintf("pivot=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods); +#endif + pivot=(int*)mxMalloc(Size*periods*sizeof(int)); +#ifdef MEM_ALLOC_CHK + mexPrintf("pivota=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods); +#endif + pivotk=(int*)mxMalloc(Size*periods*sizeof(int)); +#ifdef MEM_ALLOC_CHK + mexPrintf("pivotv=(double*)mxMalloc(%d*sizeof(double))\n",Size*periods); +#endif + pivotv=(double*)mxMalloc(Size*periods*sizeof(double)); +#ifdef MEM_ALLOC_CHK + mexPrintf("pivotva=(double*)mxMalloc(%d*sizeof(double))\n",Size*periods); +#endif + pivotva=(double*)mxMalloc(Size*periods*sizeof(double)); +#ifdef MEM_ALLOC_CHK + mexPrintf("b=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods); +#endif + b=(int*)mxMalloc(Size*periods*sizeof(int)); +#ifdef MEM_ALLOC_CHK + mexPrintf("line_done=(bool*)mxMalloc(%d*sizeof(bool))\n",Size*periods); +#endif + line_done=(bool*)mxMalloc(Size*periods*sizeof(bool)); + memset(line_done, 0, periods*Size*sizeof(*line_done)); + mem_mngr.init_CHUNK_BLCK_SIZE(u_count); + g_save_op=NULL; + g_nop_all=0; +#ifdef PRINT_OUT + mexPrintf("sizeof(NonZeroElem)=%d sizeof(NonZeroElem*)=%d\n",sizeof(NonZeroElem),sizeof(NonZeroElem*)); +#endif + i=(periods+y_kmax+1)*Size*sizeof(NonZeroElem*); +#ifdef MEM_ALLOC_CHK + mexPrintf("FNZE_R=(NonZeroElem**)mxMalloc(%d)\n",i); +#endif + FNZE_R=(NonZeroElem**)mxMalloc(i); +#ifdef MEM_ALLOC_CHK + mexPrintf("FNZE_C=(NonZeroElem**)mxMalloc(%d)\n",i); +#endif + FNZE_C=(NonZeroElem**)mxMalloc(i); + memset(FNZE_R, 0, i); + memset(FNZE_C, 0, i); +#ifdef MEM_ALLOC_CHK + mexPrintf("temp_NZE_R=(NonZeroElem**)(%d)\n",i); +#endif + NonZeroElem** temp_NZE_R=(NonZeroElem**)mxMalloc(i); +#ifdef MEM_ALLOC_CHK + mexPrintf("temp_NZE_R=(NonZeroElem**)(%d)\n",i); +#endif + NonZeroElem** temp_NZE_C=(NonZeroElem**)mxMalloc(i); + memset(temp_NZE_R, 0, i); + memset(temp_NZE_C, 0, i); + i=(periods+y_kmax+1)*Size*sizeof(int); +#ifdef MEM_ALLOC_CHK + mexPrintf("NbNZRow=(int*)mxMalloc(%d)\n",i); +#endif + NbNZRow=(int*)mxMalloc(i); +#ifdef MEM_ALLOC_CHK + mexPrintf("NbNZCol=(int*)mxMalloc(%d)\n",i); +#endif + NbNZCol=(int*)mxMalloc(i); +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + memset(NbNZRow, 0, i); + memset(NbNZCol, 0, i); + i=periods*Size*sizeof(*b); + memset(b,0,i); +#ifdef PRINT_OUT + mexPrintf("Now looping\n"); + mexEvalString("drawnow;"); +#endif + for (t=0;tfirst.first.second; + if (eq!=it4->first.first.first+Size*t) + tmp_b=0; + eq=it4->first.first.first+Size*t; +#ifdef PRINT_OUT + mexPrintf("eq=%d, var=%d",eq,var); + //mexEvalString("drawnow;"); +#endif + if (var<(periods+y_kmax)*Size) + { + lag=it4->first.second; +#ifdef PRINT_OUT + mexPrintf(", lag =%d, ti_y_kmin=%d, ti_y_kmax=%d ", lag, ti_y_kmin, ti_y_kmax); +#endif + if (lag<=ti_y_kmax && lag>=ti_y_kmin) + { + //mexPrintf("u_index=%d, eq=%d, var=%d, lag=%d ",it4->second+u_count_init*t, eq, var, lag); + var+=Size*t; + //mexPrintf("u_index=%d, eq=%d, var=%d, lag=%d ",it4->second+u_count_init*t, eq, var, lag); + NbNZRow[eq]++; + NbNZCol[var]++; +#ifdef NEW_ALLOC + first=mem_mngr.mxMalloc_NZE(); +#else + first=(NonZeroElem*)mxMalloc(sizeof(NonZeroElem)); +#endif + first->NZE_C_N=NULL; + first->NZE_R_N=NULL; + first->u_index=it4->second+u_count_init*t; + first->r_index=eq; + first->c_index=var; + first->lag_index=lag; + /*if(eq==0 && var==0) + mexPrintf("alloc FNZE_R[0]=%x\n",first);*/ + if (FNZE_R[eq]==NULL) + { + FNZE_R[eq]=first; + } + if (FNZE_C[var]==NULL) + FNZE_C[var]=first; + if (temp_NZE_R[eq]!=NULL) + temp_NZE_R[eq]->NZE_R_N=first; + if (temp_NZE_C[var]!=NULL) + temp_NZE_C[var]->NZE_C_N=first; + temp_NZE_R[eq]=first; + temp_NZE_C[var]=first; +#ifdef PRINT_OUT + mexPrintf("=> "); +#endif + } + else + { +#ifdef PRINT_OUT + mexPrintf("nn "); + mexPrintf("tmp_b+=u[%d]*y[index_var[%d]]\n",it4->second+u_count_init*t,var+Size*(y_kmin+t)); + mexPrintf("tmp_b+=u[%d](%f)*y[%d(%d)](%f)",it4->second+u_count_init*t,u[it4->second+u_count_init*t], index_vara[var+Size*(y_kmin+t)],var+Size*(y_kmin+t),y[index_vara[var+Size*(y_kmin+t)]]); + mexEvalString("drawnow;"); +#endif + tmp_b+=u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]]; + } + } + else + { +#ifdef PRINT_OUT + mexPrintf(""); +#endif + b[eq]=it4->second+u_count_init*t; + u[b[eq]]+=tmp_b; +#ifdef PRINT_OUT + mexPrintf("=> b[%d]=%f\n", eq, u[b[eq]]); + mexEvalString("drawnow;"); +#endif + } +#ifdef PRINT_OUT + mexPrintf(" u[%d] = %e\n",it4->second+u_count_init*t,double(u[it4->second+u_count_init*t])); + mexEvalString("drawnow;"); +#endif + it4++; + } + } +#ifdef PRINT_OUT + mexPrintf("end of Init\n"); + mexEvalString("drawnow;"); +#endif + mxFree(temp_NZE_R); + mxFree(temp_NZE_C); +} + +void SparseMatrix::ShortInit(int periods, int y_kmin, int y_kmax, int Size, std::map ,int>, int> IM) +{ + int t, eq, var, lag; + double tmp_b=0.0; + std::map ,int>, int>::iterator it4; + + for (t=0;tfirst.first.second; + if (eq!=it4->first.first.first+Size*t) + tmp_b=0; + eq=it4->first.first.first+Size*t; +#ifdef PRINT_OUT + mexPrintf("eq=%d, var=%d",eq,var); +#endif + if (var<(periods+y_kmax)*Size) + { + lag=it4->first.second; +#ifdef PRINT_OUT + mexPrintf(", lag =%d, ti_y_kmin=%d, ti_y_kmax=%d ", lag, ti_y_kmin, ti_y_kmax); +#endif + if (lag<=ti_y_kmax && lag>=ti_y_kmin) + { + var+=Size*t; + } + else + { +#ifdef PRINT_OUT + mexPrintf("nn "); + mexPrintf("tmp_b+=u[%d]*y[index_var[%d]]\n",it4->second+u_count_init*t,var+Size*(y_kmin+t)); + mexPrintf("tmp_b+=u[%d](%f)*y[%d(%d)](%f)",it4->second+u_count_init*t,u[it4->second+u_count_init*t], index_vara[var+Size*(y_kmin+t)],var+Size*(y_kmin+t),y[index_vara[var+Size*(y_kmin+t)]]); +#endif + tmp_b+=u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]]; + } + } + else + { +#ifdef PRINT_OUT + mexPrintf(""); +#endif + b[eq]=it4->second+u_count_init*t; + u[b[eq]]+=tmp_b; +#ifdef PRINT_OUT + mexPrintf("=> b[%d]=%f\n", eq, u[b[eq]]); +#endif + } +#ifdef PRINT_OUT + mexPrintf(" u[%d] = %e\n",it4->second+u_count_init*t,double(u[it4->second+u_count_init*t])); +#endif + it4++; + } + } +} + + + +int SparseMatrix::Get_u() +{ + if (!u_liste.empty()) + { + int i=u_liste.back(); + u_liste.pop_back(); +#ifdef PRINT_OUT + mexPrintf("Get_u=%d\n",i); +#endif + return i; + } + else + { + if (u_countNZE_R_N; + mxFree(first); + first=firsta; + } + } +#endif + mxFree(FNZE_R); + mxFree(FNZE_C); + mxFree(NbNZRow); + mxFree(NbNZCol); + mxFree(b); + mxFree(line_done); + mxFree(pivot); + mxFree(pivotk); + mxFree(pivotv); + mxFree(pivotva); +} + +bool +SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size +#ifdef PROFILER +, long int *ndiv, long int *nsub +#endif +) +{ + long int i,j,/*nop=nop4/4*/nop=nop4/2, t, index_d, k; + double r=0.0; + bool OK=true; + t_save_op_s *save_op_s, *save_opa_s, *save_opaa_s; + int *diff1, *diff2; +#ifdef MEM_ALLOC_CHK + mexPrintf("diff1=(int*)mxMalloc(%d*sizeof(int))\n",nop); +#endif + diff1=(int*)mxMalloc(nop*sizeof(int)); +#ifdef MEM_ALLOC_CHK + mexPrintf("diff1=(int*)mxMalloc(%d*sizeof(int))\n",nop); +#endif + diff2=(int*)mxMalloc(nop*sizeof(int)); +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + j=k=i=0; + while (ifirst-save_opa_s->first; + switch (save_op_s->operat) + { + case IFLD: + case IFDIV: + OK=(save_op_s->operat==save_opa_s->operat && save_opa_s->operat==save_opaa_s->operat + && diff1[j]==(save_opa_s->first-save_opaa_s->first)); + i+=2; + break; + case IFLESS: + case IFSUB: + diff2[j]=save_op_s->second-save_opa_s->second; + OK=(save_op_s->operat==save_opa_s->operat && save_opa_s->operat==save_opaa_s->operat + && diff1[j]==(save_opa_s->first-save_opaa_s->first) + && diff2[j]==(save_opa_s->second-save_opaa_s->second)); + i+=3; + break; + default: + mexPrintf("unknown operator = %d ",save_op_s->operat); + mexEvalString("st=fclose('all');clear all;"); + filename+=" stopped"; + mexErrMsgTxt(filename.c_str()); + break; + } + j++; + } +#ifdef PROFILER + if(OK) + mexPrintf("at %d same construction\n",beg_t); + else + mexPrintf("at %d different construction\n",beg_t); +#endif + // the same pivot for all remaining periods + if (OK) + for (i=beg_t;ifirst+t*diff1[j]; + if (index_d>=u_count_alloc) + { + u_count_alloc+=5*u_count_alloc_save; +#ifdef MEM_ALLOC_CHK + mexPrintf("u=(double*)mxRealloc(u,u_count_alloc*sizeof(double))\n",u_count_alloc); +#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)); + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt("Exit from Dynare"); + } + } + 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; + } + j++; + } + } + int t1=t; + for (t=t1;tlag<((periods-beg_t)-t)) + { + index_d=save_op_s->first+t*diff1[j]; + if (index_d>=u_count_alloc) + { + u_count_alloc+=5*u_count_alloc_save; +#ifdef MEM_ALLOC_CHK + mexPrintf("u=(double*)mxRealloc(u,u_count_alloc*sizeof(double))\n",u_count_alloc); +#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)); + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt("Exit from Dynare"); + } + } + 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; + } + } + else + { + switch (save_op_s->operat) + { + case IFLD : + case IFDIV : + i+=2; + break; + case IFSUB : + case IFLESS : + i+=3; + break; + } + } + j++; + } + } +#ifdef WRITE_u + toto.close(); + filename+=" stopped"; + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt(filename.c_str()); +#endif + } + mxFree(diff1); + mxFree(diff2); + return OK; +} + +void +SparseMatrix::run_u_period1(int periods) +{ + double r=0; + int index_d, t; + for (t=0;t=u_count_alloc) + { + u_count_alloc+=5*u_count_alloc_save; +#ifdef MEM_ALLOC_CHK + mexPrintf("u=(double*)mxRealloc(u,u_count_alloc*sizeof(double))\n",u_count_alloc); +#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)); + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt("Exit from Dynare"); + } + } + switch (g_save_op[i]) + { + case IFLD : + r=u[index_d]; +#ifdef PRINT_u + mexPrintf("FLD u[%d] (%f)\n",index_d,u[index_d]); +#endif + 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 + break; + case IFSUB : + u[index_d]-=u[g_save_op[i+3]+t*g_save_op[i+4]]*r; +#ifdef PRINT_u + mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f) index1=%d index2=%d\n",index_d,g_save_op[i+3]+t*g_save_op[i+4],u[g_save_op[i+3]+t*g_save_op[i+4]],r,u[index_d],g_save_op[i+1],g_save_op[i+2] ); +#endif + break; + case IFLESS: + u[index_d]=-u[g_save_op[i+3]+t*g_save_op[i+4]]*r; +#ifdef PRINT_u + mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f) index1=%d index2=%d\n",index_d,g_save_op[i+3]+t*g_save_op[i+4],u[g_save_op[i+3]+t*g_save_op[i+4]],r,u[index_d],g_save_op[i+1],g_save_op[i+2] ); +#endif + break; + } + i+=5; + } + } +} + + + + +void +SparseMatrix::run_it(int nop_all,int *op_all) +{ + double r=0; + int index_d; + for (long int i=0;i=u_count_alloc) + { + u_count_alloc+=5*u_count_alloc_save; +#ifdef MEM_ALLOC_CHK + mexPrintf("u=(double*)mxRealloc(u,u_count_alloc*sizeof(double))\n",u_count_alloc); +#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)); + mexEvalString("st=fclose('all');clear all;"); + mexErrMsgTxt("Exit from Dynare"); + } + } + switch (op_all[i]) + { + 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[op_all[i+2]]*r; +#ifdef PRINT_u + mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,op_all[i+2],u[op_all[i+2]],r,u[index_d]); +#endif + i+=3; + break; + case IFLESS: + u[index_d]=-u[op_all[i+2]]*r; +#ifdef PRINT_u + mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,op_all[i+2],u[op_all[i+2]],r,u[index_d]); +#endif + i+=3; + break; + } + } +} + + +void +SparseMatrix::run_triangular(int nop_all,int *op_all) +{ + int j=0; + mexPrintf("begining of run_triangular nop_all=%d\n",nop_all); + if (mem_mngr.swp_f) + { + bool OK=true; + int* save_op; + long int nop; + while (OK) + { + mexPrintf("reading blck%d\n",j++); + OK=mem_mngr.read_swp_f(&save_op,&nop); + if (OK) + { + run_it(nop,save_op); + mxFree(save_op); + } + } + } + run_it(nop_all,op_all); +} + +int +SparseMatrix::complete(int beg_t, int Size, int periods, int *b) +{ + long int i, j, k, nop, nopa, nop1, cal_y, nb_var, pos, t, ti, max_var, min_var; + NonZeroElem *first; + int *save_code; + int *diff; + double yy=0.0, err; + + int size_of_save_code=(1+y_kmax)*Size*(Size+1+4)/2*4; +#ifdef MEM_ALLOC_CHK + mexPrintf("save_code=(int*)mxMalloc(%d*sizeof(int))\n",size_of_save_code); +#endif + save_code=(int*)mxMalloc(size_of_save_code*sizeof(int)); + int size_of_diff=(1+y_kmax)*Size*(Size+1+4); +#ifdef MEM_ALLOC_CHK + mexPrintf("diff=(int*)mxMalloc(%d*sizeof(int))\n",size_of_diff); +#endif + diff=(int*)mxMalloc(size_of_diff*sizeof(int)); +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + cal_y=y_size*y_kmin; + + i=(beg_t+1)*Size-1; + nop=0; + for (j=i;j>i-Size;j--) + { + pos=pivot[j]; + nb_var=At_Row(pos,&first); + first=first->NZE_R_N; + nb_var--; + save_code[nop]=IFLDZ; + save_code[nop+1]=0; + save_code[nop+2]=0; + save_code[nop+3]=0; + if ((nop+3)>=size_of_save_code) + mexPrintf("out of save_code[%d] (bound=%d)\n",nop+2,size_of_save_code); + nop+=4; + for (k=0;kc_index]+cal_y; + save_code[nop+2]=first->u_index; + save_code[nop+3]=first->lag_index; + if ((nop+3)>=size_of_save_code) + mexPrintf("out of save_code[%d] (bound=%d)\n",nop+2,size_of_save_code); + nop+=4; + first=first->NZE_R_N; + } + //yy=-(yy+u[b[pos]]); + //mexPrintf("|u[%d]|\n",b[pos]); + save_code[nop]=IFADD; + save_code[nop+1]=b[pos]; + save_code[nop+2]=0; + save_code[nop+3]=0; + if ((nop+3)>=size_of_save_code) + mexPrintf("out of save_code[%d] (bound=%d)\n",nop+2,size_of_save_code); + nop+=4; + save_code[nop]=IFSTP; + save_code[nop+1]=index_vara[j]+y_size*y_kmin; + save_code[nop+2]=0; + save_code[nop+3]=0; + if ((nop+2)>=size_of_save_code) + mexPrintf("out of save_code[%d] (bound=%d)\n",nop+2,size_of_save_code); + nop+=4; + } + i=beg_t*Size-1; + nop1=nopa=0; + for (j=i;j>i-Size;j--) + { + pos=pivot[j]; + nb_var=At_Row(pos,&first); + first=first->NZE_R_N; + nb_var--; + diff[nopa]=0; + diff[nopa+1]=0; + nopa+=2; + nop1+=4; + for (k=0;kc_index]+cal_y); + diff[nopa+1]=save_code[nop1+2]-(first->u_index); + if ((nop1+2)>=size_of_save_code) + mexPrintf("out of save_code[%d] (bound=%d)\n",nop1+2,size_of_save_code); + if ((nopa+1)>=size_of_diff) + mexPrintf("out of diff[%d] (bound=%d)\n",nopa+2,size_of_diff); + nopa+=2; + nop1+=4; + first=first->NZE_R_N; + } + diff[nopa]=save_code[nop1+1]-(b[pos]); + diff[nopa+1]=0; + if ((nop1+3)>=size_of_save_code) + mexPrintf("out of save_code[%d] (bound=%d)\n",nop1+2,size_of_save_code); + if ((nopa+1)>=size_of_diff) + mexPrintf("out of diff[%d] (bound=%d)\n",nopa+2,size_of_diff); + nopa+=2; + nop1+=4; + diff[nopa]=save_code[nop1+1]-(index_vara[j]+y_size*y_kmin); + diff[nopa+1]=0; + if ((nop1+4)>=size_of_save_code) + mexPrintf("out of save_code[%d] (bound=%d)\n",nop1+2,size_of_save_code); + if ((nopa+1)>=size_of_diff) + mexPrintf("out of diff[%d] (bound=%d)\n",nopa+2,size_of_diff); + nopa+=2; + nop1+=4; + } + max_var=(periods+y_kmin)*y_size; + min_var=y_kmin*y_size; + int k1=0; + for (t=periods+y_kmin-1;t>=beg_t+y_kmin;t--) + { + j=0; + ti=t-y_kmin-beg_t; + for (i=0;imin_var) + { + yy+=y[k]*u[save_code[i+2]+ti*diff[j+1]]; +#ifdef PRINT_OUT_y1 + mexPrintf("y[%d](%f)*u[%d](%f)+",k, double(y[k]), save_code[i+2]+ti*diff[j+1], double(u[save_code[i+2]+ti*diff[j+1]])); +#endif + } + break; + case IFADD : + yy=-(yy+u[save_code[i+1]+ti*diff[j]]); +#ifdef PRINT_OUT_y1 + mexPrintf("|u[%d](%f)|",save_code[i+1]+ti*diff[j],double(u[save_code[i+1]+ti*diff[j]])); +#endif + break; + case IFSTP : + k=save_code[i+1]+ti*diff[j]; + k1=k; + err = yy - y[k]; + y[k] += slowc*(err); +#ifdef PRINT_OUT_y1 + mexPrintf("=y[%d]=%f diff[%d]=%d save_code[%d]=%d ti=%d\n",save_code[i+1]+ti*diff[j],y[k],j,diff[j],i+1,save_code[i+1],ti); +#endif + break; + } + j+=2; + } + } + mxFree(save_code); + mxFree(diff); + return(beg_t); +} + +double +SparseMatrix::bksub( int tbreak, int last_period, int Size, double slowc_l +#ifdef PROFILER +, /*NonZeroElem *first,*/ long int *nmul +#endif +) +{ + NonZeroElem *first; + int i, j, k; + double yy; + res1 = res2 = max_res = 0; + for (i=0;i=y_kmin;t--) + { + int ti=(t-y_kmin)*Size; +#ifdef PRINT_OUT + mexPrintf("t=%d ti=%d\n",t,ti); +#endif + int cal=y_kmin*Size; + int cal_y=y_size*y_kmin; + for (i=ti-1;i>=ti-Size;i--) + { + j=i+cal; +#ifdef PRINT_OUT_y + mexPrintf("t=%d, ti=%d j=%d i+Size=%d\n",t,ti,j,i+Size); +#endif + int pos=pivot[/*j*/i+Size]; +#ifdef PRINT_OUT_y + mexPrintf("i-ti+Size=%d pos=%d j=%d\n",i-ti+Size,pos,j); +#endif + int nb_var=At_Row(pos,&first); + first=first->NZE_R_N; + nb_var--; + int eq=index_vara[j]+y_size; +#ifdef PRINT_OUT_y1 + mexPrintf("y[index_vara[%d]=%d]=",j,index_vara[j]+y_size); +#endif + yy=0; + for (k=0;kc_index]+cal_y]*u[first->u_index]; +#ifdef PROFILER + (*nmul)++; +#endif + first=first->NZE_R_N; + } +#ifdef PRINT_OUT_y1 + mexPrintf("|u[%d](%f)|",b[pos],double(u[b[pos]])); +#endif + yy=-(yy+y[eq]+u[b[pos]]); + direction[eq]=yy; + y[eq] += slowc_l*yy; +#ifdef PRINT_OUT_y1 + mexPrintf("=%f (%f)\n",double(yy),double(y[eq])); +#endif + } + } + return res1; +} + + +void +SparseMatrix::close_swp_file() +{ + mem_mngr.close_swp_f(); +} + + + + + +void +SparseMatrix::Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter) +{ + int i, j, k, l, m=0, nop; + double err; +#ifdef PRINT_OUT_y + int m1; +#endif + int period = it_ * y_size, s_middle_count_loop = 0 ; + if (periods>0) + 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); +#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;"); + } +} + + + + + + + + + + + +int +SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter) +{ + /*Triangularisation at each period of a block using a simple gaussian Elimination*/ + t_save_op_s *save_op_s; + int start_compare=y_kmin; + bool record=false; + int *save_op=NULL, *save_opa=NULL, *save_opaa=NULL; + long int nop=0, nopa=0; + int tbreak=0, last_period=periods; + int i,j,k; + int pivj=0, pivk=0; + NonZeroElem *first, *firsta, *first_sub, *first_piv, *first_suba; + double piv_abs, first_elem; + //SparseMatrix sparse_matrix; + //mexPrintf("->u_count=%d &u_count=%x\n",u_count,&u_count); +#ifdef RECORD_ALL + int save_u_count=u_count; +#endif + u_count_alloc_save=u_count_alloc; +#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; +#endif + //pctimer_t t01; + clock_t t01; + //pctimer_t t1 = pctimer(); + clock_t t1 = clock(); +#ifdef PROFILER + tdelete1=0; tdelete2=0; tdelete21=0; tdelete22=0; tdelete221=0; tdelete222=0; +#endif + if (iter>0) + mexPrintf("Sim : %f ms\n",(1000.0*(double(clock())-double(time00)))/double(CLOCKS_PER_SEC)); +#ifdef MEMORY_LEAKS + mexEvalString("feature('memstats');"); +#endif + if (isnan(res1) || isinf(res1)) + { + if (slowc_save<1e-8) + { + mexPrintf("Dynare cannot improve the simulation\n"); + mexEvalString("drawnow;"); + mexEvalString("st=fclose('all');clear all;"); + filename+=" stopped"; + mexErrMsgTxt(filename.c_str()); + } + slowc_save/=2; + mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n",slowc_save); + for (i=0;i-0.3) && symbolic) + { + if (start_compare==y_kmin) + { + mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied for a longer period.\n"); + start_compare=min(tbreak_g,periods); + } + else + { + mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied to all periods.\n"); + symbolic=false; + alt_symbolic=true; + markowitz_c_s=markowitz_c; + markowitz_c=0; + } + } + res1a=res1; + 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"); + if (cvg) + return(0); +#ifdef PRINT_OUT + mexPrintf("Size=%d y_size=%d y_kmin=%d y_kmax=%d u_count=%d u_count_alloc=%d periods=%d\n",Size,y_size,y_kmin,y_kmax,u_count,u_count_alloc,periods); + mexEvalString("drawnow;"); +#endif +#ifdef PROFILER + clock_t t00 = clock(); +#endif +#ifdef WRITE_u + fstream toto; + int i_toto=0; + if (!symbolic) + { + toto.open("compare_good.txt", std::ios::out); + } +#endif +#ifdef PROFILER + t01=clock(); +#endif +#ifdef PROFILER + mexPrintf("initialization time=%f ms\n",1000.0*(double(t01)-double(t00))/double(CLOCKS_PER_SEC)); +#endif + +#ifdef PRINT_OUT + mexPrintf("sizeof(NonZeroElem)=%d\n",sizeof(NonZeroElem)); + /*for (i=0;i0) + { +#ifdef PRINT_OUT + mexPrintf("run_triangular\n"); + mexEvalString("drawnow;"); +#endif + run_u_period1(periods); +#ifdef PRINT_OUT + mexPrintf("done\n"); + mexEvalString("drawnow;"); +#endif + } + else + { +#ifdef PRINT_OUT + mexPrintf("Init\n"); + mexEvalString("drawnow;"); +#endif + Init(periods, y_kmin, y_kmax, Size, IM_i); +#ifdef PRINT_OUT + mexPrintf("done\n"); + mexEvalString("drawnow;"); +#endif + for (int t=0;tr_index,line_done[first->r_index]); +#endif + if (!line_done[first->r_index]) + { + k=first->u_index; +#ifdef PRINT_OUT + mexPrintf("u[%d]=%f fabs(u[%d])=%f\n",k,double(u[k]),k,double(fabs(u[k]))); +#endif + int jj=first->r_index; + int NRow_jj=NRow(jj); +#ifdef PROFILER + nbpivot++; + nbpivot_it++; +#endif + +#ifdef MARKOVITZ + piv_v[l]=u[k]; + double piv_fabs=fabs(u[k]); + pivj_v[l]=jj; + pivk_v[l]=k; + NR[l]=NRow_jj; + if (NRow_jj==1 && !one) + { + one=true; + piv_abs=piv_fabs; + N_max=NRow_jj; + } + if (!one) + { + if (piv_fabs>piv_abs) + piv_abs=piv_fabs; + if (NRow_jj>N_max) + N_max=NRow_jj; + } + else + { + if (NRow_jj==1) + { + if (piv_fabs>piv_abs) + piv_abs=piv_fabs; + if (NRow_jj>N_max) + N_max=NRow_jj; + } + } + l++; +#else + if (piv_absNZE_C_N; + } +#ifdef MARKOVITZ + double markovitz=0, markovitz_max=-9e70; + if (!one) + { + for (j=0;jmarkovitz_max) + { + piv=piv_v[j]; + pivj=pivj_v[j]; //Line number + pivk=pivk_v[j]; //positi + markovitz_max=markovitz; + } + } + } + else + { + for (j=0;jmarkovitz_max && NR[j]==1) + { + piv=piv_v[j]; + pivj=pivj_v[j]; //Line number + pivk=pivk_v[j]; //positi + markovitz_max=markovitz; + } + } + } +#endif +#ifdef PROFILER + tpivot+=clock()-td0; +#endif + +#ifdef PRINT_OUT + mexPrintf("Thats the pivot: %d with value %f in u[%d] \n",pivj,double(piv),pivk); + mexPrintf("______________________________________________\n"); + mexPrintf("pivot[%d]=%d\n",i,pivj); + mexEvalString("drawnow;"); +#endif + + pivot[i]=pivj; + pivotk[i]=pivk; + pivotv[i]=piv; + } + else + { + pivj=pivot[i-Size]+Size; + pivot[i]=pivj; + At_Pos(pivj, i, &first); + pivk=first->u_index; + piv=u[pivk]; + piv_abs=fabs(piv); + } + line_done[pivj]=true; +#ifdef PRINT_u + mexPrintf("FLD u[%d] (%f=%f) |",pivk,u[pivk],piv); + Print_u();mexPrintf("\n"); + mexEvalString("drawnow;"); +#endif + if (symbolic) + { + if (record) + { + if (nop+1>=nopa) + { + nopa=int(1.5*nopa); +#ifdef MEM_ALLOC_CHK + mexPrintf("save_op=(int*)mxRealloc(save_op,%d*sizeof(int))\n",nopa); +#endif +#ifdef PROFILER + nbRealloc++; +#endif +#ifdef N_MX_ALLOC + save_op=realloc_std(save_op, nopa); +#else + save_op=(int*)mxRealloc(save_op,nopa*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + } + save_op_s=(t_save_op_s*)(&(save_op[nop])); + save_op_s->operat=IFLD; + save_op_s->first=pivk; + save_op_s->lag=0; + } + nop+=2; + } +#ifdef RECORD_ALL + else if (record_all) + { + if (nop_all+1>=nopa_all) + chk_avail_mem(&save_op_all,&nop_all,&nopa_all,1,t); + save_op_all[nop_all]=IFLD; + save_op_all[nop_all+1]=pivk; + nop_all+=2; + } +#endif + if (piv_abslag_index, first->r_index, first->c_index, first->u_index); +#endif + u[first->u_index]/=piv; +#ifdef PROFILER + nbdiv++; +#endif +#ifdef WRITE_u + if ((periods-t)<=y_kmax) + { + toto << i_toto << " u[" /*<< first->u_index*/ << "]/=" << piv << "=" << u[first->u_index] << endl; + i_toto++; + } +#endif +#ifdef PRINT_u + mexPrintf("FDIV u[%d](%f)/=piv(%f)=(%f) |",first->u_index,u[first->u_index],piv,u[first->u_index]); + Print_u();mexPrintf("\n"); +#endif + if (symbolic) + { + if (record) + { + if (nop+1>=nopa) + { + nopa=int(1.5*nopa); +#ifdef MEM_ALLOC_CHK + mexPrintf("save_op=(int*)mxRealloc(save_op,%d*sizeof(int))\n",nopa); +#endif +#ifdef PROFILER + nbRealloc++; +#endif +#ifdef N_MX_ALLOC + save_op=realloc_std(save_op, nopa); +#else + save_op=(int*)mxRealloc(save_op,nopa*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + } + save_op_s=(t_save_op_s*)(&(save_op[nop])); + save_op_s->operat=IFDIV; + save_op_s->first=first->u_index; + save_op_s->lag=first->lag_index; + } + nop+=2; + } +#ifdef RECORD_ALL + else if (record_all) + { + if (nop_all+1>=nopa_all) + chk_avail_mem(&save_op_all,&nop_all,&nopa_all,1,t); + save_op_all[nop_all]=IFDIV; + save_op_all[nop_all+1]=first->u_index; + nop_all+=2; + } +#endif + first=first->NZE_R_N; + } +#ifdef PRINT_OUT + mexPrintf("dividing at u[%d]\n",b[pivj]); +#endif + u[b[pivj]]/=piv; +#ifdef WRITE_u + if ((periods-t)<=y_kmax) + { + toto << i_toto << " u[" /*<< b[pivj]*/ << "]/=" << piv << "=" << u[b[pivj]] << endl; + i_toto++; + } +#endif +#ifdef PRINT_OUT + mexPrintf("ok3\n"); + mexEvalString("drawnow;"); +#endif +#ifdef PRINT_u + mexPrintf("FDIV u[%d](%f)/=piv(%f)=(%f) |",b[pivj],u[b[pivj]],piv,u[b[pivj]]); + Print_u();mexPrintf("\n"); +#endif + if (symbolic) + { + if (record) + { + if (nop+1>=nopa) + { + nopa=int(1.5*nopa); +#ifdef MEM_ALLOC_CHK + mexPrintf("save_op=(int*)mxRealloc(save_op,%d*sizeof(int))\n",nopa); +#endif +#ifdef PROFILER + nbRealloc++; +#endif +#ifdef N_MX_ALLOC + save_op=realloc_std(save_op, nopa); +#else + save_op=(int*)mxRealloc(save_op,nopa*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + } + save_op_s=(t_save_op_s*)(&(save_op[nop])); + save_op_s->operat=IFDIV; + save_op_s->first=b[pivj]; + save_op_s->lag=0; + } + nop+=2; + } +#ifdef RECORD_ALL + else if (record_all) + { + if (nop_all+1>=nopa_all) + chk_avail_mem(&save_op_all,&nop_all,&nopa_all,1,t); + save_op_all[nop_all]=IFDIV; + save_op_all[nop_all+1]=b[pivj]; + nop_all+=2; + } +#endif + /*substract the elements on the non treated lines*/ + nb_eq=At_Col(i,&first); + NonZeroElem *first_piva; + int nb_var_piva=At_Row(pivj,&first_piva); +#ifdef PRINT_OUT + if(iter>0) + { + mexPrintf("ok4 nb_eq=%d iter=%d\n",nb_eq,iter); + mexEvalString("drawnow;"); + } +#endif + for (j=0;jr_index; +#ifdef PRINT_OUT + mexPrintf("line_done[%d]=%d\n", row, int(line_done[row])); +#endif + if (!line_done[row]) + { +#ifdef PRINT_OUT + mexPrintf("Substracting from line %d lag %d\n",row,first->lag_index); +#endif + first_elem=u[first->u_index]; +#ifdef PRINT_u + mexPrintf("FLD u[%d] (%f) |",first->u_index,u[first->u_index]); + Print_u();mexPrintf("\n"); +#endif + if (symbolic) + { + if (record) + { + if (nop+1>=nopa) + { + nopa=int(1.5*nopa); +#ifdef MEM_ALLOC_CHK + mexPrintf("save_op=(int*)mxRealloc(save_op,%d*sizeof(int))\n",nopa); +#endif +#ifdef PROFILER + nbRealloc++; +#endif +#ifdef N_MX_ALLOC + save_op=realloc_std(save_op, nopa); +#else + save_op=(int*)mxRealloc(save_op,nopa*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + } + save_op_s=(t_save_op_s*)(&(save_op[nop])); + save_op_s->operat=IFLD; + save_op_s->first=first->u_index; + save_op_s->lag=abs(first->lag_index); + } + nop+=2; + } +#ifdef RECORD_ALL + else if (record_all) + { + if (nop_all+1>=nopa_all) + chk_avail_mem(&save_op_all,&nop_all,&nopa_all,1,t); + save_op_all[nop_all]=IFLD; + save_op_all[nop_all+1]=first->u_index; + nop_all+=2; + } +#endif + int nb_var_piv=nb_var_piva; + first_piv=first_piva; + int nb_var_sub=At_Row(row,&first_sub); + int l_sub=0, l_piv=0; + int sub_c_index=first_sub->c_index, piv_c_index=first_piv->c_index; + int tmp_lag=first_sub->lag_index; +#ifdef PROFILER + td1 = clock(); +#endif + while (l_subr_index, first_piv->lag_index, first_piv->c_index,l_piv); + mexPrintf("l_sub(%d)r_index, first_sub->lag_index, first_sub->c_index,l_sub); +#endif + if (l_sub=nb_var_piv)) + { +#ifdef PRINT_OUT + if(iter>0) + { + mexPrintf(" u[%d]=u[%d]\n",first_sub->u_index,first_sub->u_index); + mexEvalString("drawnow;"); + } +#endif + first_sub=first_sub->NZE_R_N; + if (first_sub) + sub_c_index=first_sub->c_index; + else + sub_c_index=Size*periods; + l_sub++; + } + else if (sub_c_index>piv_c_index || l_sub>=nb_var_sub) + { + int tmp_u_count=Get_u(); +#ifdef PROFILER + clock_t td0=clock(); +#endif + int lag=first_piv->c_index/Size-row/Size; + Insert(row,first_piv->c_index,tmp_u_count,lag); +#ifdef PROFILER + tinsert+=clock()-td0; +#endif + u[tmp_u_count]=-u[first_piv->u_index]*first_elem; +#ifdef PROFILER + nbless++; +#endif +#ifdef WRITE_u + if ((periods-t)<=y_kmax) + { + toto << i_toto << " u[" /*<< tmp_u_count*/ << "]=-u[" /*<< first_piv->u_index*/ << "]*" << first_elem << "=" << u[tmp_u_count] << endl; + i_toto++; + } +#endif +#ifdef PRINT_u + if(iter>0) + { + mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f) |",tmp_u_count,first_piv->u_index,u[first_piv->u_index],first_elem,u[tmp_u_count]); + /*Print_u();*/mexPrintf("\n"); + mexEvalString("drawnow;"); + } + +#endif + if (symbolic) + { + if (record) + { + if (nop+2>=nopa) + { + nopa=int(1.5*nopa); +#ifdef MEM_ALLOC_CHK + mexPrintf("save_op=(int*)mxRealloc(save_op,%d*sizeof(int))\n",nopa); +#endif +#ifdef PROFILER + nbRealloc++; +#endif +#ifdef N_MX_ALLOC + save_op=realloc_std(save_op, nopa); +#else + save_op=(int*)mxRealloc(save_op,nopa*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + } + save_op_s=(t_save_op_s*)(&(save_op[nop])); + save_op_s->operat=IFLESS; + save_op_s->first=tmp_u_count; + save_op_s->second=first_piv->u_index; + save_op_s->lag=max(first_piv->lag_index,abs(tmp_lag)); + } + nop+=3; + } +#ifdef RECORD_ALL + else if (record_all) + { + if (nop_all+2>=nopa_all) + chk_avail_mem(&save_op_all,&nop_all,&nopa_all,2,t); + save_op_all[nop_all]=IFLESS; + save_op_all[nop_all+1]=tmp_u_count; + save_op_all[nop_all+2]=first_piv->u_index; + nop_all+=3; + } +#endif +#ifdef PRINT_OUT + if(iter>0) + { + mexPrintf(" u[%d at (%d, %d) lag %d]=-u[%d]*%f\n",tmp_u_count,row , first_piv->c_index, first_piv->c_index/Size-row/Size, first_piv->u_index ,double(first_elem)); + mexEvalString("drawnow;"); + } + +#endif + first_piv=first_piv->NZE_R_N; + if (first_piv) + piv_c_index=first_piv->c_index; + else + piv_c_index=Size*periods; + l_piv++; + } + else /*first_sub->c_index==first_piv->c_index*/ + { + if (i==sub_c_index) + { +#ifdef PRINT_OUT + /*if(iter>0) + { + mexPrintf(" delete element [%d, %d] lag %d u[%d]\n",first_sub->r_index,first_sub->c_index,first_sub->lag_index,first_sub->u_index); + mexEvalString("drawnow;"); + }*/ +#endif + firsta=first; + first_suba=first_sub->NZE_R_N; + /*if(iter>0) + { + mexPrintf("just after\n"); + mexEvalString("drawnow;"); + }*/ +#ifdef PROFILER + clock_t td0=clock(); +#endif + Delete(first_sub->r_index,first_sub->c_index, Size, b); +#ifdef PROFILER + tdelete+=clock()-td0; +#endif + first=firsta->NZE_C_N; + first_sub=first_suba; + if (first_sub) + sub_c_index=first_sub->c_index; + else + sub_c_index=Size*periods; + l_sub++; + first_piv=first_piv->NZE_R_N; + if (first_piv) + piv_c_index=first_piv->c_index; + else + piv_c_index=Size*periods; + l_piv++; +#ifdef PRINT_OUT + Print(Size,b); +#endif + } + else + { +#ifdef PRINT_OUT + mexPrintf(" u[%d]-=u[%d]*%f\n",first_sub->u_index,first_piv->u_index,double(first_elem)); +#endif + + + u[first_sub->u_index]-=u[first_piv->u_index]*first_elem; +#ifdef PROFILER + nbless++; +#endif +#ifdef WRITE_u + if ((periods-t)<=y_kmax) + { + toto << i_toto << " u[" /*<< first_sub->u_index*/ << "]-=u[" /*<< first_piv->u_index*/ << "]*" << first_elem << "=" << u[first_sub->u_index] << endl; + i_toto++; + } +#endif +#ifdef PRINT_u + if(iter>0) + { + mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f) |",first_sub->u_index,first_piv->u_index,u[first_piv->u_index],first_elem,u[first_sub->u_index]); + /*Print_u();*/mexPrintf("\n"); + mexEvalString("drawnow;"); + } +#endif + if (symbolic) + { + if (record) + { + if (nop+2>=nopa) + { + nopa=int(1.5*nopa); +#ifdef MEM_ALLOC_CHK + mexPrintf("save_op=(int*)mxRealloc(save_op,%d*sizeof(int))\n",nopa); +#endif +#ifdef PROFILER + nbRealloc++; +#endif +#ifdef N_MX_ALLOC + save_op=realloc_std(save_op, nopa); +#else + save_op=(int*)mxRealloc(save_op,nopa*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + } + save_op_s=(t_save_op_s*)(&(save_op[nop])); + save_op_s->operat=IFSUB; + save_op_s->first=first_sub->u_index; + save_op_s->second=first_piv->u_index; + save_op_s->lag=max(abs(tmp_lag),first_piv->lag_index); + } + nop+=3; + } +#ifdef RECORD_ALL + else if (record_all) + { + if (nop_all+2>=nopa_all) + chk_avail_mem(&save_op_all,&nop_all,&nopa_all,2,t); + save_op_all[nop_all]=IFSUB; + save_op_all[nop_all+1]=first_sub->u_index; + save_op_all[nop_all+2]=first_piv->u_index; + nop_all+=3; + } +#endif + first_sub=first_sub->NZE_R_N; + if (first_sub) + sub_c_index=first_sub->c_index; + else + sub_c_index=Size*periods; + l_sub++; + first_piv=first_piv->NZE_R_N; + if (first_piv) + piv_c_index=first_piv->c_index; + else + piv_c_index=Size*periods; + l_piv++; + } + } + } +#ifdef PRINT_OUT + mexPrintf("row=%d pivj=%d\n",row,pivj); + mexPrintf(" u[%d](%f)-=u[%d](%f)*%f\n",b[row],double(u[b[row]]), b[pivj],double(u[b[pivj]]),double(first_elem)); +#endif + + u[b[row]]-=u[b[pivj]]*first_elem; +#ifdef PROFILER + nbless++; + tbigloop += clock() - td1; +#endif +#ifdef WRITE_u + if ((periods-t)<=y_kmax) + { + toto << i_toto << " u[" /*<< b[row]*/ << "]-=u[" /*<< b[pivj]*/ << "]*" << first_elem << "=" << u[b[row]] << endl; + i_toto++; + } +#endif +#ifdef PRINT_u + mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f) |",b[row],b[pivj],u[b[pivj]],first_elem,u[b[row]]); + Print_u();mexPrintf("\n"); +#endif + /*if(iter>0) + { + mexPrintf("ok4f j=%d\n",j); + mexEvalString("drawnow;"); + }*/ + + if (symbolic) + { + if (record) + { + if (nop+2>=nopa) + { + nopa=int(1.5*nopa); +#ifdef MEM_ALLOC_CHK + mexPrintf("save_op=(int*)mxRealloc(save_op,%d*sizeof(int))\n",nopa); +#endif +#ifdef PROFILER + nbRealloc++; +#endif +#ifdef N_MX_ALLOC + save_op=realloc_std(save_op, nopa); +#else + save_op=(int*)mxRealloc(save_op,nopa*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + } + save_op_s=(t_save_op_s*)(&(save_op[nop])); + save_op_s->operat=IFSUB; + save_op_s->first=b[row]; + save_op_s->second=b[pivj]; + save_op_s->lag=abs(tmp_lag); + } + nop+=3; + } +#ifdef RECORD_ALL + else if (record_all) + { + if (nop_all+2>=nopa_all) + chk_avail_mem(&save_op_all,&nop_all,&nopa_all,2,t); + save_op_all[nop_all]=IFSUB; + save_op_all[nop_all+1]=b[row]; + save_op_all[nop_all+2]=b[pivj]; + nop_all+=3; + } + if(iter>0) + { + mexPrintf("ok4g j=%d\n",j); + mexEvalString("drawnow;"); + } + +#endif + } + else + first=first->NZE_C_N; + +#ifdef PRINT_OUT + mexPrintf(" bef first=%x\n",first); +#endif + } + } + if (symbolic) + { +#ifdef PROFILER + td1=clock(); + mexPrintf("at %d nop=%d ?= nop1=%d record=%d save_opa=%x save_opaa=%x \n", t, nop, nop1, record, save_opa, save_opaa); +#endif + if (record && (nop==nop1)) + { +#ifdef PRINT_OUT + mexPrintf("nop=%d, nop1=%d, record=%d save_opa=%x save_opaa=%x\n",nop, nop1, record, save_opa, save_opaa); +#endif + if (save_opa && save_opaa) + { +#ifdef PROFILER + clock_t ta=clock(); +#endif + if (compare( save_op, save_opa, save_opaa, t, periods, nop, Size +#ifdef PROFILER + , &ndiv, &nsub +#endif + )) + { +#ifdef PROFILER + mexPrintf("t=%d over periods=%d t0=%f t1=%f y_kmin=%d y_kmax=%d\n",t,periods,1000*(ta-t00), 1000.0*(double(clock())-double(ta))/double(CLOCKS_PER_SEC),y_kmin, y_kmax); + mexPrintf("compare time %f ms\n",1000.0*(double(clock())-double(ta))/double(CLOCKS_PER_SEC)); +#endif + tbreak=t; + tbreak_g=tbreak; + break; + } + } + if (save_opa) + { + if (save_opaa) + { +#ifdef N_MX_ALLOC + free(save_opaa); +#else + mxFree(save_opaa); +#endif + save_opaa=NULL; + } +#ifdef MEM_ALLOC_CHK + mexPrintf("save_opaa=(int*)mxMalloc(%d*sizeof(int))\n",nop1); +#endif +#ifdef N_MX_ALLOC + save_opaa=malloc_std(nop1); +#else + save_opaa=(int*)mxMalloc(nop1*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + memcpy(save_opaa, save_opa, nop1*sizeof(int)); + } + if (save_opa) + { +#ifdef N_MX_ALLOC + free(save_opa); +#else + mxFree(save_opa); +#endif + save_opa=NULL; + } +#ifdef MEM_ALLOC_CHK + mexPrintf("save_opa=(int*)mxMalloc(%d*sizeof(int))\n",nop); +#endif +#ifdef N_MX_ALLOC + save_opa=malloc_std(nop); +#else + save_opa=(int*)mxMalloc(nop*sizeof(int)); +#endif +#ifdef MEM_ALLOC_CHK + mexPrintf("ok\n"); +#endif + memcpy(save_opa, save_op, nop*sizeof(int)); + } + else + { + if (nop==nop1) + record=true; + else + { + record=false; + if (save_opa) + { +#ifdef N_MX_ALLOC + free(save_opa); +#else + mxFree(save_opa); +#endif + save_opa=NULL; + } + if (save_opaa) + { +#ifdef N_MX_ALLOC + free(save_opaa); +#else + mxFree(save_opaa); +#endif + save_opaa=NULL; + } + } + } + nop2=nop1; + nop1=nop; +#ifdef PROFILER + tcompare+=clock()-td1; +#endif + } + record=true; + } + } + nop_all+=nop; +#ifdef PROFILER + t01=clock(); + mexPrintf("resolution time=%f ms ndiv=%d nsub=%d ncomp=%d \n",1000.0*(double(t01)-double(t00))/double(CLOCKS_PER_SEC),ndiv,nsub,ncomp); + mexPrintf(" tinsert=%f tdelete=%f tpivot=%f tbigloop=%f tdelete1=%f tdelete2=%f \n",double(1000*tinsert), double(1000*tdelete), double(1000*tpivot), double(1000*tbigloop), double(1000*tdelete1), double(1000*tdelete2)); + mexPrintf(" tdelete21=%f tdelete22=%f \n",double(1000*tdelete21),double(1000*tdelete22)); + mexPrintf(" tdelete221=%f tdelete222=%f \n",double(1000*tdelete221),double(1000*tdelete222)); + mexPrintf(" tcompare=%f \n",double(1000*tcompare)); + mexPrintf("nbpivot=%d, nbdiv=%d, nbless=%d, nop=%d nbpivot_it=%d nbRealloc=%d\n", nbpivot, nbdiv, nbless, nbpivot + nbdiv + nbless, nbpivot_it, nbRealloc); +#endif + if (symbolic) + { +#ifdef N_MX_ALLOC + if (save_op) + free(save_op); + if (save_opa) + free(save_opa); + if (save_opaa) + free(save_opaa); +#else + if (save_op) + mxFree(save_op); + if (save_opa) + mxFree(save_opa); + if (save_opaa) + mxFree(save_opaa); +#endif + } + close_swp_file(); + + /*The backward substitution*/ +#ifdef PRINT_OUT + Print(Size,b); +#endif + double slowc_lbx=slowc, res1bx; +#ifdef PROFILER + t00 = clock(); +#endif + for (i=0;i0)) + { + u_count=save_u_count; + End(Size); + } +#endif + if (print_it) + { + //pctimer_t t2 = pctimer(); + clock_t t2 = clock(); + mexPrintf("(** %f milliseconds **)\n", 1000.0*(double(t2) - double(t1))/double(CLOCKS_PER_SEC)); + mexEvalString("drawnow;"); + } + +#ifdef WRITE_u + if (!symbolic) + { + toto.close(); + mexEvalString("st=fclose('all');clear all;"); + filename+=" stopped"; + mexErrMsgTxt(filename.c_str()); + } +#endif + close_swp_file(); + time00=clock(); + return(0); +} + +void +SparseMatrix::fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1) +{ + u_count=u_count_int * periods; + u_count_alloc = 2*u_count; + (*u)=(double*)mxMalloc(u_count_alloc*sizeof(double)); + memset((*u), 0, u_count_alloc*sizeof(double)); + u_count_init=max_lag_plus_max_lead_plus_1; +} + + +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; +} +*/