diff --git a/mex/sources/build_matlab.m b/mex/sources/build_matlab.m index 561e4b634..81c53639b 100644 --- a/mex/sources/build_matlab.m +++ b/mex/sources/build_matlab.m @@ -36,9 +36,9 @@ if strcmpi('GLNX86', computer) || strcmpi('GLNXA64', computer) ... elseif strcmpi('PCWIN', computer) || strcmpi('PCWIN64', computer) % Windows (x86-32 or x86-64) with Microsoft or gcc compiler if strcmpi('PCWIN', computer) - LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win32/microsoft/']; + LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win32/microsoft/']; else - LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win64/microsoft/']; + LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win64/microsoft/']; end LAPACK_PATH = ['"' LIBRARY_PATH 'libmwlapack.lib"']; if matlab_ver_less_than('7.5') @@ -59,7 +59,7 @@ COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -DMATLAB_MEX_FILE -DMATLAB_VERSION=0x' spr % Large array dims for 64 bits platforms appeared in Matlab 7.3 if (strcmpi('GLNXA64', computer) || strcmpi('PCWIN64', computer)) ... - && ~matlab_ver_less_than('7.3') + && ~matlab_ver_less_than('7.3') COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -largeArrayDims' ]; end diff --git a/mex/sources/build_matlab_multithread.m b/mex/sources/build_matlab_multithread.m index aaf4886ba..057bb4b08 100644 --- a/mex/sources/build_matlab_multithread.m +++ b/mex/sources/build_matlab_multithread.m @@ -36,9 +36,9 @@ if strcmpi('GLNX86', computer) || strcmpi('GLNXA64', computer) ... elseif strcmpi('PCWIN', computer) || strcmpi('PCWIN64', computer) % Windows (x86-32 or x86-64) with Microsoft or gcc compiler if strcmpi('PCWIN', computer) - LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win32/microsoft/']; + LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win32/microsoft/']; else - LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win64/microsoft/']; + LIBRARY_PATH = [MATLAB_PATH '/extern/lib/win64/microsoft/']; end LAPACK_PATH = ['"' LIBRARY_PATH 'libmwlapack.lib"']; if matlab_ver_less_than('7.5') @@ -58,7 +58,7 @@ COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -DMATLAB_MEX_FILE -DMATLAB_VERSION=0x' spr % Large array dims for 64 bits platforms appeared in Matlab 7.3 if (strcmpi('GLNXA64', computer) || strcmpi('PCWIN64', computer)) ... - && ~matlab_ver_less_than('7.3') + && ~matlab_ver_less_than('7.3') COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -largeArrayDims' ]; end @@ -122,7 +122,7 @@ eval([ COMPILE_COMMAND ' -I. mjdgges/mjdgges.c ' LAPACK_PATH ]); disp('Compiling sparse_hessian_times_B_kronecker_C...') eval([ COMPILE_COMMAND_OMP ' -I. kronecker/sparse_hessian_times_B_kronecker_C.cc' ]); - + disp('Compiling A_times_B_kronecker_C...') eval([ COMPILE_COMMAND_OMP ' -I. kronecker/A_times_B_kronecker_C.cc ']); diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index 821a631c4..9b13454a4 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -18,7 +18,7 @@ */ #include #include -#include "Interpreter.hh" +#include "Interpreter.hh" #define BIG 1.0e+8; #define SMALL 1.0e-5; @@ -28,28 +28,28 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub 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, int minimal_solving_periods_arg) { - params=params_arg; - y=y_arg; - ya=ya_arg; - x=x_arg; + params = params_arg; + y = y_arg; + ya = ya_arg; + x = x_arg; steady_y = steady_y_arg; steady_x = steady_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; + 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; slowc_save = slowc; - y_decal=y_decal_arg; - markowitz_c=markowitz_c_arg; - filename=filename_arg; - T=NULL; + y_decal = y_decal_arg; + markowitz_c = markowitz_c_arg; + filename = filename_arg; + T = NULL; error_not_printed = true; minimal_solving_periods = minimal_solving_periods_arg; } @@ -57,17 +57,17 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub double Interpreter::pow1(double a, double b) { - double r = pow_(a, b); + double r = pow_(a, b); if (isnan(r) || isinf(r)) { - if(a<0 && error_not_printed) - { - mexPrintf("Error: X^a with X=%5.15f\n",a); - error_not_printed = false; - r = 0.0000000000000000000000001; - } - res1=NAN; - return(r); + if (a < 0 && error_not_printed) + { + mexPrintf("Error: X^a with X=%5.15f\n", a); + error_not_printed = false; + r = 0.0000000000000000000000001; + } + res1 = NAN; + return (r); } else return r; @@ -76,16 +76,16 @@ Interpreter::pow1(double a, double b) double Interpreter::log1(double a) { - double r = log(a); + double r = log(a); if (isnan(r) || isinf(r)) { - if(a<=0 && error_not_printed) - { - mexPrintf("Error: log(X) with X=%5.15f\n",a); - error_not_printed = false; - } - res1=NAN; - return(r); + if (a <= 0 && error_not_printed) + { + mexPrintf("Error: log(X) with X=%5.15f\n", a); + error_not_printed = false; + } + res1 = NAN; + return (r); } else return r; @@ -97,575 +97,575 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) int var, lag = 0, op; ostringstream tmp_out; double v1, v2; - bool go_on=true; + bool go_on = true; double ll; while (go_on) { switch (it_code->first) { - case FLDV : - //load a variable in the processor - switch (((FLDV_*)it_code->second)->get_type()) - { - case eParameter : - var = ((FLDV_*)it_code->second)->get_pos(); - Stack.push(params[var]); + case FLDV: + //load a variable in the processor + switch (((FLDV_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FLDV_ *) it_code->second)->get_pos(); + Stack.push(params[var]); #ifdef DEBUG - tmp_out << " params[" << var << "](" << params[var] << ")"; + tmp_out << " params[" << var << "](" << params[var] << ")"; #endif - break; - case eEndogenous : - var = ((FLDV_*)it_code->second)->get_pos(); - lag = ((FLDV_*)it_code->second)->get_lead_lag(); - if(evaluate) - Stack.push(ya[(it_+lag)*y_size+var]); - else - { - Stack.push(y[(it_+lag)*y_size+var]); - } + break; + case eEndogenous: + var = ((FLDV_ *) it_code->second)->get_pos(); + lag = ((FLDV_ *) it_code->second)->get_lead_lag(); + if (evaluate) + Stack.push(ya[(it_+lag)*y_size+var]); + else + { + Stack.push(y[(it_+lag)*y_size+var]); + } #ifdef DEBUG - tmp_out << " y[" << it_+lag << ", " << var << "](" << y[(it_+lag)*y_size+var] << ")"; + tmp_out << " y[" << it_+lag << ", " << var << "](" << y[(it_+lag)*y_size+var] << ")"; #endif - break; - case eExogenous : - var = ((FLDV_*)it_code->second)->get_pos(); - lag = ((FLDV_*)it_code->second)->get_lead_lag(); - Stack.push(x[it_+lag+var*nb_row_x]); + break; + case eExogenous: + var = ((FLDV_ *) it_code->second)->get_pos(); + lag = ((FLDV_ *) it_code->second)->get_lead_lag(); + Stack.push(x[it_+lag+var*nb_row_x]); #ifdef DEBUG - tmp_out << " x[" << it_+lag << ", " << var << "](" << x[it_+lag+var*nb_row_x] << ")"; + tmp_out << " x[" << it_+lag << ", " << var << "](" << x[it_+lag+var*nb_row_x] << ")"; #endif - break; - case eExogenousDet : - var = ((FLDV_*)it_code->second)->get_pos(); - lag = ((FLDV_*)it_code->second)->get_lead_lag(); - Stack.push(x[it_+lag+var*nb_row_xd]); - break; - case eModelLocalVariable: + break; + case eExogenousDet: + var = ((FLDV_ *) it_code->second)->get_pos(); + lag = ((FLDV_ *) it_code->second)->get_lead_lag(); + Stack.push(x[it_+lag+var*nb_row_xd]); + break; + case eModelLocalVariable: #ifdef DEBUG - mexPrintf("FLDV a local variable in Block %d Stack.size()=%d",block_num, Stack.size()); - mexPrintf(" value=%f\n",Stack.top()); + mexPrintf("FLDV a local variable in Block %d Stack.size()=%d", block_num, Stack.size()); + mexPrintf(" value=%f\n", Stack.top()); #endif - break; - default: - mexPrintf("FLDV: Unknown variable type\n"); - } - break; - case FLDSV : - //load a variable in the processor - switch (((FLDSV_*)it_code->second)->get_type()) - { - case eParameter : - var = ((FLDSV_*)it_code->second)->get_pos(); - Stack.push(params[var]); + break; + default: + mexPrintf("FLDV: Unknown variable type\n"); + } + break; + case FLDSV: + //load a variable in the processor + switch (((FLDSV_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FLDSV_ *) it_code->second)->get_pos(); + Stack.push(params[var]); #ifdef DEBUG - tmp_out << " params[" << var << "](" << params[var] << ")"; + tmp_out << " params[" << var << "](" << params[var] << ")"; #endif - break; - case eEndogenous : - var = ((FLDSV_*)it_code->second)->get_pos(); - if(evaluate) - Stack.push(ya[var]); - else - Stack.push(y[var]); + break; + case eEndogenous: + var = ((FLDSV_ *) it_code->second)->get_pos(); + if (evaluate) + Stack.push(ya[var]); + else + Stack.push(y[var]); #ifdef DEBUG - tmp_out << " y[" << var << "](" << y[var] << ")"; + tmp_out << " y[" << var << "](" << y[var] << ")"; #endif - break; - case eExogenous : - var = ((FLDSV_*)it_code->second)->get_pos(); - Stack.push(x[var]); + break; + case eExogenous: + var = ((FLDSV_ *) it_code->second)->get_pos(); + Stack.push(x[var]); #ifdef DEBUG - tmp_out << " x[" << var << "](" << x[var] << ")"; + tmp_out << " x[" << var << "](" << x[var] << ")"; #endif - break; - case eExogenousDet : - var = ((FLDSV_*)it_code->second)->get_pos(); - Stack.push(x[var]); - break; - case eModelLocalVariable: + break; + case eExogenousDet: + var = ((FLDSV_ *) it_code->second)->get_pos(); + Stack.push(x[var]); + break; + case eModelLocalVariable: #ifdef DEBUG - mexPrintf("FLDSV a local variable in Block %d Stack.size()=%d",block_num, Stack.size()); - mexPrintf(" value=%f\n",Stack.top()); + mexPrintf("FLDSV a local variable in Block %d Stack.size()=%d", block_num, Stack.size()); + mexPrintf(" value=%f\n", Stack.top()); #endif - break; - default: - mexPrintf("FLDSV: Unknown variable type\n"); - } - break; - case FLDVS : - //load a variable in the processor - switch (((FLDVS_*)it_code->second)->get_type()) - { - case eParameter : - var = ((FLDVS_*)it_code->second)->get_pos(); + break; + default: + mexPrintf("FLDSV: Unknown variable type\n"); + } + break; + case FLDVS: + //load a variable in the processor + switch (((FLDVS_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FLDVS_ *) it_code->second)->get_pos(); #ifdef DEBUG - mexPrintf("params[%d]=%f\n", var, params[var]); + mexPrintf("params[%d]=%f\n", var, params[var]); #endif - Stack.push(params[var]); - break; - case eEndogenous : - var = ((FLDVS_*)it_code->second)->get_pos(); + Stack.push(params[var]); + break; + case eEndogenous: + var = ((FLDVS_ *) it_code->second)->get_pos(); #ifdef DEBUG - mexPrintf(" steady_y[%d]=%f\n", var, steady_y[var]); + mexPrintf(" steady_y[%d]=%f\n", var, steady_y[var]); #endif - Stack.push(steady_y[var]); - break; - case eExogenous : - var = ((FLDVS_*)it_code->second)->get_pos(); + Stack.push(steady_y[var]); + break; + case eExogenous: + var = ((FLDVS_ *) it_code->second)->get_pos(); #ifdef DEBUG - mexPrintf(" x[%d] = %f\n", var, x[var]); + mexPrintf(" x[%d] = %f\n", var, x[var]); #endif - Stack.push(x[var]); - break; - case eExogenousDet : - var = ((FLDVS_*)it_code->second)->get_pos(); + Stack.push(x[var]); + break; + case eExogenousDet: + var = ((FLDVS_ *) it_code->second)->get_pos(); #ifdef DEBUG - mexPrintf(" xd[%d] = %f\n", var, x[var]); + mexPrintf(" xd[%d] = %f\n", var, x[var]); #endif - Stack.push(x[var]); - break; - case eModelLocalVariable: + Stack.push(x[var]); + break; + case eModelLocalVariable: #ifdef DEBUG - mexPrintf("FLDVS a local variable in Block %d Stack.size()=%d",block_num, Stack.size()); - mexPrintf(" value=%f\n",Stack.top()); + mexPrintf("FLDVS a local variable in Block %d Stack.size()=%d", block_num, Stack.size()); + mexPrintf(" value=%f\n", Stack.top()); #endif - break; - default: - mexPrintf("FLDVS: Unknown variable type\n"); - } - break; - case FLDT : - //load a temporary variable in the processor - var = ((FLDT_*)it_code->second)->get_pos(); + break; + default: + mexPrintf("FLDVS: Unknown variable type\n"); + } + break; + case FLDT: + //load a temporary variable in the processor + var = ((FLDT_ *) it_code->second)->get_pos(); #ifdef DEBUG - tmp_out << " T[" << it_ << ", " << var << "](" << T[var*(periods+y_kmin+y_kmax)+it_] << ")"; + tmp_out << " T[" << it_ << ", " << var << "](" << T[var*(periods+y_kmin+y_kmax)+it_] << ")"; #endif - Stack.push(T[var*(periods+y_kmin+y_kmax)+it_]); - break; - case FLDST : - //load a temporary variable in the processor - var = ((FLDST_*)it_code->second)->get_pos(); + Stack.push(T[var*(periods+y_kmin+y_kmax)+it_]); + break; + case FLDST: + //load a temporary variable in the processor + var = ((FLDST_ *) it_code->second)->get_pos(); #ifdef DEBUG - tmp_out << " T[" << var << "](" << T[var] << ")"; + tmp_out << " T[" << var << "](" << T[var] << ")"; #endif - Stack.push(T[var]); - break; - case FLDU : - //load u variable in the processor - var = ((FLDU_*)it_code->second)->get_pos(); - var+=Per_u_; + Stack.push(T[var]); + break; + case FLDU: + //load u variable in the processor + var = ((FLDU_ *) it_code->second)->get_pos(); + var += Per_u_; #ifdef DEBUG - tmp_out << " u[" << var << "](" << u[var] << ")"; + tmp_out << " u[" << var << "](" << u[var] << ")"; #endif - Stack.push(u[var]); - break; - case FLDSU : - //load u variable in the processor - var = ((FLDSU_*)it_code->second)->get_pos(); + Stack.push(u[var]); + break; + case FLDSU: + //load u variable in the processor + var = ((FLDSU_ *) it_code->second)->get_pos(); #ifdef DEBUG - tmp_out << " u[" << var << "](" << u[var] << ")"; + tmp_out << " u[" << var << "](" << u[var] << ")"; #endif - Stack.push(u[var]); - break; - case FLDR : - //load u variable in the processor - var = ((FLDR_*)it_code->second)->get_pos(); - Stack.push(r[var]); - break; - case FLDZ : - //load 0 in the processor - Stack.push(0); + Stack.push(u[var]); + break; + case FLDR: + //load u variable in the processor + var = ((FLDR_ *) it_code->second)->get_pos(); + Stack.push(r[var]); + break; + case FLDZ: + //load 0 in the processor + Stack.push(0); #ifdef DEBUG - tmp_out << " 0"; + tmp_out << " 0"; #endif - break; - case FLDC : - //load a numerical constant in the processor - ll = ((FLDC_*)it_code->second)->get_value(); + break; + case FLDC: + //load a numerical constant in the processor + ll = ((FLDC_ *) it_code->second)->get_value(); #ifdef DEBUG - tmp_out << " " << ll; + tmp_out << " " << ll; #endif - Stack.push(ll); - break; - case FSTPV : - //load a variable in the processor - switch (((FSTPV_*)it_code->second)->get_type()) - { - case eParameter : - var = ((FSTPV_*)it_code->second)->get_pos(); - params[var] = Stack.top(); - Stack.pop(); - break; - case eEndogenous : - var = ((FSTPV_*)it_code->second)->get_pos(); - lag = ((FSTPV_*)it_code->second)->get_lead_lag(); - y[(it_+lag)*y_size+var] = Stack.top(); + Stack.push(ll); + break; + case FSTPV: + //load a variable in the processor + switch (((FSTPV_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FSTPV_ *) it_code->second)->get_pos(); + params[var] = Stack.top(); + Stack.pop(); + break; + case eEndogenous: + var = ((FSTPV_ *) it_code->second)->get_pos(); + lag = ((FSTPV_ *) it_code->second)->get_lead_lag(); + y[(it_+lag)*y_size+var] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" y[%d, %d](%f)=%s\n", it_+lag, var, y[(it_+lag)*y_size+var], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" y[%d, %d](%f)=%s\n", it_+lag, var, y[(it_+lag)*y_size+var], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case eExogenous : - var = ((FSTPV_*)it_code->second)->get_pos(); - lag = ((FSTPV_*)it_code->second)->get_lead_lag(); - x[it_+lag+var*nb_row_x] = Stack.top(); + Stack.pop(); + break; + case eExogenous: + var = ((FSTPV_ *) it_code->second)->get_pos(); + lag = ((FSTPV_ *) it_code->second)->get_lead_lag(); + x[it_+lag+var*nb_row_x] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" x[%d, %d](%f)=%s\n", it_+lag, var, x[it_+lag+var*nb_row_x], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" x[%d, %d](%f)=%s\n", it_+lag, var, x[it_+lag+var*nb_row_x], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case eExogenousDet : - var = ((FSTPV_*)it_code->second)->get_pos(); - lag = ((FSTPV_*)it_code->second)->get_lead_lag(); - x[it_+lag+var*nb_row_xd] = Stack.top(); + Stack.pop(); + break; + case eExogenousDet: + var = ((FSTPV_ *) it_code->second)->get_pos(); + lag = ((FSTPV_ *) it_code->second)->get_lead_lag(); + x[it_+lag+var*nb_row_xd] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" x[%d, %d](%f)=%s\n", it_+lag, var, x[it_+lag+var*nb_row_xd], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" x[%d, %d](%f)=%s\n", it_+lag, var, x[it_+lag+var*nb_row_xd], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - default: - mexPrintf("FSTPV: Unknown variable type\n"); - } - break; - case FSTPSV : - //load a variable in the processor - switch (((FSTPSV_*)it_code->second)->get_type()) - { - case eParameter : - var = ((FSTPSV_*)it_code->second)->get_pos(); - params[var] = Stack.top(); - Stack.pop(); - break; - case eEndogenous : - var = ((FSTPSV_*)it_code->second)->get_pos(); - y[var] = Stack.top(); + Stack.pop(); + break; + default: + mexPrintf("FSTPV: Unknown variable type\n"); + } + break; + case FSTPSV: + //load a variable in the processor + switch (((FSTPSV_ *) it_code->second)->get_type()) + { + case eParameter: + var = ((FSTPSV_ *) it_code->second)->get_pos(); + params[var] = Stack.top(); + Stack.pop(); + break; + case eEndogenous: + var = ((FSTPSV_ *) it_code->second)->get_pos(); + y[var] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" y[%d](%f)=%s\n", var, y[var], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" y[%d](%f)=%s\n", var, y[var], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case eExogenous : - case eExogenousDet : - var = ((FSTPSV_*)it_code->second)->get_pos(); - x[var] = Stack.top(); + Stack.pop(); + break; + case eExogenous: + case eExogenousDet: + var = ((FSTPSV_ *) it_code->second)->get_pos(); + x[var] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" x[%d, %d](%f)=%s\n", it_+lag, var, x[var], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" x[%d, %d](%f)=%s\n", it_+lag, var, x[var], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - default: - mexPrintf("FSTPSV: Unknown variable type\n"); - } - break; - case FSTPT : - //store in a temporary variable from the processor - var = ((FSTPT_*)it_code->second)->get_pos(); - T[var*(periods+y_kmin+y_kmax)+it_] = Stack.top(); + Stack.pop(); + break; + default: + mexPrintf("FSTPSV: Unknown variable type\n"); + } + break; + case FSTPT: + //store in a temporary variable from the processor + var = ((FSTPT_ *) it_code->second)->get_pos(); + T[var*(periods+y_kmin+y_kmax)+it_] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" T[%d, %d](%f)=%s\n", it_, var, T[var*(periods+y_kmin+y_kmax)+it_], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" T[%d, %d](%f)=%s\n", it_, var, T[var*(periods+y_kmin+y_kmax)+it_], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case FSTPST : - //store in a temporary variable from the processor - var = ((FSTPST_*)it_code->second)->get_pos(); - T[var] = Stack.top(); + Stack.pop(); + break; + case FSTPST: + //store in a temporary variable from the processor + var = ((FSTPST_ *) it_code->second)->get_pos(); + T[var] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" T[%d](%f)=%s\n", var, T[var], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" T[%d](%f)=%s\n", var, T[var], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case FSTPU : - //store in u variable from the processor - var = ((FSTPU_*)it_code->second)->get_pos(); - var+=Per_u_; - u[var] = Stack.top(); + Stack.pop(); + break; + case FSTPU: + //store in u variable from the processor + var = ((FSTPU_ *) it_code->second)->get_pos(); + var += Per_u_; + u[var] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" u[%d](%f)=%s\n", var, u[var], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" u[%d](%f)=%s\n", var, u[var], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case FSTPSU : - //store in u variable from the processor - var = ((FSTPSU_*)it_code->second)->get_pos(); - u[var] = Stack.top(); + Stack.pop(); + break; + case FSTPSU: + //store in u variable from the processor + var = ((FSTPSU_ *) it_code->second)->get_pos(); + u[var] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" u[%d](%f)=%s\n", var, u[var], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" u[%d](%f)=%s\n", var, u[var], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case FSTPR : - //store in residual variable from the processor - var = ((FSTPR_*)it_code->second)->get_pos(); - r[var] = Stack.top(); + Stack.pop(); + break; + case FSTPR: + //store in residual variable from the processor + var = ((FSTPR_ *) it_code->second)->get_pos(); + r[var] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" r[%d](%f)=%s\n", var, r[var], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" r[%d](%f)=%s\n", var, r[var], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case FSTPG : - //store in derivative (g) variable from the processor - var = ((FSTPG_*)it_code->second)->get_pos(); - g1[var] = Stack.top(); + Stack.pop(); + break; + case FSTPG: + //store in derivative (g) variable from the processor + var = ((FSTPG_ *) it_code->second)->get_pos(); + g1[var] = Stack.top(); #ifdef DEBUG - tmp_out << "=>"; - mexPrintf(" g1[%d](%f)=%s\n", var, g1[var], tmp_out.str().c_str()); - tmp_out.str(""); + tmp_out << "=>"; + mexPrintf(" g1[%d](%f)=%s\n", var, g1[var], tmp_out.str().c_str()); + tmp_out.str(""); #endif - Stack.pop(); - break; - case FBINARY : - op = ((FBINARY_*)it_code->second)->get_op_type(); - v2=Stack.top(); - Stack.pop(); - v1=Stack.top(); - Stack.pop(); - switch (op) - { - case oPlus: - Stack.push(v1 + v2); + Stack.pop(); + break; + case FBINARY: + op = ((FBINARY_ *) it_code->second)->get_op_type(); + v2 = Stack.top(); + Stack.pop(); + v1 = Stack.top(); + Stack.pop(); + switch (op) + { + case oPlus: + Stack.push(v1 + v2); #ifdef DEBUG - tmp_out << " |" << v1 << "+" << v2 << "|"; + tmp_out << " |" << v1 << "+" << v2 << "|"; #endif - break; - case oMinus: - Stack.push(v1 - v2); + break; + case oMinus: + Stack.push(v1 - v2); #ifdef DEBUG - tmp_out << " |" << v1 << "-" << v2 << "|"; + tmp_out << " |" << v1 << "-" << v2 << "|"; #endif - break; - case oTimes: - Stack.push(v1 * v2); + break; + case oTimes: + Stack.push(v1 * v2); #ifdef DEBUG - tmp_out << " |" << v1 << "*" << v2 << "|"; + tmp_out << " |" << v1 << "*" << v2 << "|"; #endif - break; - case oDivide: - Stack.push(v1 / v2); + break; + case oDivide: + Stack.push(v1 / v2); #ifdef DEBUG - tmp_out << " |" << v1 << "/" << v2 << "|"; + tmp_out << " |" << v1 << "/" << v2 << "|"; #endif - break; - case oLess: - Stack.push(double(v1v2)); + break; + case oGreater: + Stack.push(double (v1 > v2)); #ifdef DEBUG - tmp_out << " |" << v1 << ">" << v2 << "|"; + tmp_out << " |" << v1 << ">" << v2 << "|"; #endif - break; - case oLessEqual: - Stack.push(double(v1<=v2)); + break; + case oLessEqual: + Stack.push(double (v1 <= v2)); #ifdef DEBUG - tmp_out << " |" << v1 << "<=" << v2 << "|"; + tmp_out << " |" << v1 << "<=" << v2 << "|"; #endif - break; - case oGreaterEqual: - Stack.push(double(v1>=v2)); + break; + case oGreaterEqual: + Stack.push(double (v1 >= v2)); #ifdef DEBUG - tmp_out << " |" << v1 << ">=" << v2 << "|"; + tmp_out << " |" << v1 << ">=" << v2 << "|"; #endif - break; - case oEqualEqual: - Stack.push(double(v1==v2)); + break; + case oEqualEqual: + Stack.push(double (v1 == v2)); #ifdef DEBUG - tmp_out << " |" << v1 << "==" << v2 << "|"; + tmp_out << " |" << v1 << "==" << v2 << "|"; #endif - break; - case oDifferent: - Stack.push(double(v1!=v2)); + break; + case oDifferent: + Stack.push(double (v1 != v2)); #ifdef DEBUG - tmp_out << " |" << v1 << "!=" << v2 << "|"; + tmp_out << " |" << v1 << "!=" << v2 << "|"; #endif - break; - case oPower: - Stack.push(pow1(v1, v2)); + break; + case oPower: + Stack.push(pow1(v1, v2)); #ifdef DEBUG - tmp_out << " |" << v1 << "^" << v2 << "|"; + tmp_out << " |" << v1 << "^" << v2 << "|"; #endif - break; - case oMax: - Stack.push(max(v1, v2)); + break; + case oMax: + Stack.push(max(v1, v2)); #ifdef DEBUG - tmp_out << " |max(" << v1 << "," << v2 << ")|"; + tmp_out << " |max(" << v1 << "," << v2 << ")|"; #endif - break; - case oMin: - Stack.push(min(v1, v2)); + break; + case oMin: + Stack.push(min(v1, v2)); #ifdef DEBUG - tmp_out << " |min(" << v1 << "," << v2 << ")|"; + tmp_out << " |min(" << v1 << "," << v2 << ")|"; #endif - break; - case oEqual: - default: - /*throw EvalException();*/ - ; - } - break; - case FUNARY : - op = ((FUNARY_*)it_code->second)->get_op_type(); - v1=Stack.top(); - Stack.pop(); - switch (op) - { - case oUminus: - Stack.push(-v1); + break; + case oEqual: + default: + /*throw EvalException();*/ + ; + } + break; + case FUNARY: + op = ((FUNARY_ *) it_code->second)->get_op_type(); + v1 = Stack.top(); + Stack.pop(); + switch (op) + { + case oUminus: + Stack.push(-v1); #ifdef DEBUG - tmp_out << " |-(" << v1 << ")|"; + tmp_out << " |-(" << v1 << ")|"; #endif - break; - case oExp: - Stack.push(exp(v1)); + break; + case oExp: + Stack.push(exp(v1)); #ifdef DEBUG - tmp_out << " |exp(" << v1 << ")|"; + tmp_out << " |exp(" << v1 << ")|"; #endif - break; - case oLog: - Stack.push(log1(v1)); + break; + case oLog: + Stack.push(log1(v1)); #ifdef DEBUG - tmp_out << " |log(" << v1 << ")|"; + tmp_out << " |log(" << v1 << ")|"; #endif - break; - case oLog10: - Stack.push(log10(v1)); + break; + case oLog10: + Stack.push(log10(v1)); #ifdef DEBUG - tmp_out << " |log10(" << v1 << ")|"; + tmp_out << " |log10(" << v1 << ")|"; #endif - break; - case oCos: - Stack.push(cos(v1)); + break; + case oCos: + Stack.push(cos(v1)); #ifdef DEBUG - tmp_out << " |cos(" << v1 << ")|"; + tmp_out << " |cos(" << v1 << ")|"; #endif - break; - case oSin: - Stack.push(sin(v1)); + break; + case oSin: + Stack.push(sin(v1)); #ifdef DEBUG - tmp_out << " |sin(" << v1 << ")|"; + tmp_out << " |sin(" << v1 << ")|"; #endif - break; - case oTan: - Stack.push(tan(v1)); + break; + case oTan: + Stack.push(tan(v1)); #ifdef DEBUG - tmp_out << " |tan(" << v1 << ")|"; + tmp_out << " |tan(" << v1 << ")|"; #endif - break; - case oAcos: - Stack.push(acos(v1)); + break; + case oAcos: + Stack.push(acos(v1)); #ifdef DEBUG - tmp_out << " |acos(" << v1 << ")|"; + tmp_out << " |acos(" << v1 << ")|"; #endif - break; - case oAsin: - Stack.push(asin(v1)); + break; + case oAsin: + Stack.push(asin(v1)); #ifdef DEBUG - tmp_out << " |asin(" << v1 << ")|"; + tmp_out << " |asin(" << v1 << ")|"; #endif - break; - case oAtan: - Stack.push(atan(v1)); + break; + case oAtan: + Stack.push(atan(v1)); #ifdef DEBUG - tmp_out << " |atan(" << v1 << ")|"; + tmp_out << " |atan(" << v1 << ")|"; #endif - break; - case oCosh: - Stack.push(cosh(v1)); + break; + case oCosh: + Stack.push(cosh(v1)); #ifdef DEBUG - tmp_out << " |cosh(" << v1 << ")|"; + tmp_out << " |cosh(" << v1 << ")|"; #endif - break; - case oSinh: - Stack.push(sinh(v1)); + break; + case oSinh: + Stack.push(sinh(v1)); #ifdef DEBUG - tmp_out << " |sinh(" << v1 << ")|"; + tmp_out << " |sinh(" << v1 << ")|"; #endif - break; - case oTanh: - Stack.push(tanh(v1)); + break; + case oTanh: + Stack.push(tanh(v1)); #ifdef DEBUG - tmp_out << " |tanh(" << v1 << ")|"; + tmp_out << " |tanh(" << v1 << ")|"; #endif - break; - case oAcosh: - Stack.push(acosh(v1)); + break; + case oAcosh: + Stack.push(acosh(v1)); #ifdef DEBUG - tmp_out << " |acosh(" << v1 << ")|"; + tmp_out << " |acosh(" << v1 << ")|"; #endif - break; - case oAsinh: - Stack.push(asinh(v1)); + break; + case oAsinh: + Stack.push(asinh(v1)); #ifdef DEBUG - tmp_out << " |asinh(" << v1 << ")|"; + tmp_out << " |asinh(" << v1 << ")|"; #endif - break; - case oAtanh: - Stack.push(atanh(v1)); + break; + case oAtanh: + Stack.push(atanh(v1)); #ifdef DEBUG - tmp_out << " |atanh(" << v1 << ")|"; + tmp_out << " |atanh(" << v1 << ")|"; #endif - break; - case oSqrt: - Stack.push(sqrt(v1)); + break; + case oSqrt: + Stack.push(sqrt(v1)); #ifdef DEBUG - tmp_out << " |sqrt(" << v1 << ")|"; + tmp_out << " |sqrt(" << v1 << ")|"; #endif - break; - default: - ; - } - break; - case FCUML : - v1=Stack.top(); - Stack.pop(); - v2=Stack.top(); - Stack.pop(); - Stack.push(v1+v2); - break; - case FENDBLOCK : - //it's the block end - go_on=false; - break; - case FENDEQU : - break; - case FOK : - op = ((FOK_*)it_code->second)->get_arg(); - if (Stack.size()>0) - { - mexPrintf("error: Stack not empty!\n"); - mexErrMsgTxt("End of bytecode"); - } - break; - default : - mexPrintf("Unknown opcode %d!! FENDEQU=%d\n",it_code->first,FENDEQU); - mexErrMsgTxt("End of bytecode"); - break; + break; + default: + ; + } + break; + case FCUML: + v1 = Stack.top(); + Stack.pop(); + v2 = Stack.top(); + Stack.pop(); + Stack.push(v1+v2); + break; + case FENDBLOCK: + //it's the block end + go_on = false; + break; + case FENDEQU: + break; + case FOK: + op = ((FOK_ *) it_code->second)->get_arg(); + if (Stack.size() > 0) + { + mexPrintf("error: Stack not empty!\n"); + mexErrMsgTxt("End of bytecode"); + } + break; + default: + mexPrintf("Unknown opcode %d!! FENDEQU=%d\n", it_code->first, FENDEQU); + mexErrMsgTxt("End of bytecode"); + break; } it_code++; } @@ -678,149 +678,149 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam it_code_type begining; switch (type) { - case EVALUATE_FORWARD : - if(steady_state) - compute_block_time(0, true, block_num); - else - { - begining=it_code; - for (it_=y_kmin;it_=y_kmin;it_--) - { - it_code = begining; - Per_y_=it_*y_size; - compute_block_time(0, true, block_num); - } - } - break; - case SOLVE_BACKWARD_SIMPLE : - g1=(double*)mxMalloc(size*size*sizeof(double)); - r=(double*)mxMalloc(size*sizeof(double)); - if(steady_state) - { - compute_block_time(0, true, block_num); - for(int j=0; j=y_kmin;it_--) - { - it_code = begining; - Per_y_=it_*y_size; - compute_block_time(0,true, block_num); - for(int j=0; j=y_kmin;it_--) - { - it_code = begining; - Per_y_=it_*y_size; - compute_block_time(0,true, block_num); - for(int j=0; j= y_kmin; it_--) + { + it_code = begining; + Per_y_ = it_*y_size; + compute_block_time(0, true, block_num); + } + } + break; + case SOLVE_BACKWARD_SIMPLE: + g1 = (double *) mxMalloc(size*size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + if (steady_state) + { + compute_block_time(0, true, block_num); + for (int j = 0; j < size; j++) + y[Block_Contain[j].Variable] += r[j]; + } + else + { + begining = it_code; + for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--) + { + it_code = begining; + Per_y_ = it_*y_size; + compute_block_time(0, true, block_num); + for (int j = 0; j < size; j++) + y[it_*y_size+Block_Contain[j].Variable] += r[j]; + } + } + mxFree(g1); + mxFree(r); + break; + case SOLVE_BACKWARD_COMPLETE: + fixe_u(&u, u_count_int, u_count_int); + Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false); + r = (double *) mxMalloc(size*sizeof(double)); + if (steady_state) + { + compute_block_time(0, true, block_num); + for (int j = 0; j < size; j++) + y[Block_Contain[j].Variable] += r[j]; + } + else + { + begining = it_code; + for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--) + { + it_code = begining; + Per_y_ = it_*y_size; + compute_block_time(0, true, block_num); + for (int j = 0; j < size; j++) + y[it_*y_size+Block_Contain[j].Variable] += r[j]; + } + } + mxFree(r); + break; + case SOLVE_TWO_BOUNDARIES_SIMPLE: + case SOLVE_TWO_BOUNDARIES_COMPLETE: + fixe_u(&u, u_count_int, u_count_int); + Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true); + u_count = u_count_int*(periods+y_kmax+y_kmin); + r = (double *) mxMalloc(size*sizeof(double)); + begining = it_code; + for (it_ = y_kmin; it_ < periods+y_kmin; it_++) + { + Per_u_ = (it_-y_kmin)*u_count_int; + Per_y_ = it_*y_size; + it_code = begining; + compute_block_time(Per_u_, true, block_num); + for (int j = 0; j < size; j++) + y[it_*y_size+Block_Contain[j].Variable] += r[j]; + } + mxFree(r); + break; } } @@ -836,532 +836,532 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, double *y_save; switch (type) { - case EVALUATE_FORWARD : - if(steady_state) - compute_block_time(0, false, block_num); - else - { - begining = it_code; - for (it_=y_kmin;it_=y_kmin;it_--) - { - it_code = begining; - Per_y_=it_*y_size; - compute_block_time(0, false, block_num); - } - } - break; - case SOLVE_FORWARD_SIMPLE : - g1=(double*)mxMalloc(size*size*sizeof(double)); - r=(double*)mxMalloc(size*sizeof(double)); - begining = it_code; - if(steady_state) - { - cvg=false; - iter=0; - while (!(cvg||(iter>maxit_))) - { - it_code = begining; - Per_y_=it_*y_size; - compute_block_time(0, false, block_num); - double rr; - rr=r[0]; - cvg=(fabs(rr)maxit_))) - { - it_code = begining; - Per_y_=it_*y_size; - compute_block_time(0, false, block_num); - double rr; - if(fabs(1+y[Per_y_+Block_Contain[0].Variable])>eps) - rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]); - else - rr=r[0]; - cvg=(fabs(rr)maxit_))) - { - it_code = begining; - Per_y_=it_*y_size; - compute_block_time(0, false, block_num); - double rr; - rr=r[0]; - cvg=(fabs(rr)y_kmin;it_--) - { - cvg=false; - iter=0; - Per_y_=it_*y_size; - while (!(cvg||(iter>maxit_))) - { - it_code = begining; - Per_y_=it_*y_size; - compute_block_time(0, false, block_num); - double rr; - if(fabs(1+y[Per_y_+Block_Contain[0].Variable])>eps) - rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]); - else - rr=r[0]; - cvg=(fabs(rr)maxit_))) - { - it_code= begining; - error_not_printed = true; - res2=0; - res1=0; - max_res=0; - compute_block_time(0, false, block_num); - if (!(isnan(res1)||isinf(res1))) - { - for (i=0; imaxit_))) - { - it_code = begining; - error_not_printed = true; - res2=0; - res1=0; - max_res=0; - compute_block_time(0, false, block_num); - if (!(isnan(res1)||isinf(res1))) - { - for (i=0; ieps) - rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); - else - rr=r[i]; - if (max_resmaxit_))) - { - it_code = begining; - error_not_printed = true; - res2=0; - res1=0; - max_res=0; - compute_block_time(0, false, block_num); - if (!(isnan(res1)||isinf(res1))) - { - for (i=0; iy_kmin;it_--) - { - cvg=false; - iter=0; - Per_y_=it_*y_size; - while (!(cvg||(iter>maxit_))) - { - it_code = begining; - error_not_printed = true; - res2=0; - res1=0; - max_res=0; - compute_block_time(0, false, block_num); - if (!(isnan(res1)||isinf(res1))) - { - for (i=0; ieps) - rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); - else - rr=r[i]; - if (max_resy_kmin;it_--) - { - it_code = begining; - Per_y_=it_*y_size; - error_not_printed = true; - compute_block_time(0, false, block_num); - cvg=false; - result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false); - } - } - } - mxFree(index_equa); - mxFree(index_vara); - memset(direction,0,size_of_direction); - mxFree(g1); - mxFree(r); - mxFree(u); - break; - case SOLVE_TWO_BOUNDARIES_SIMPLE : - case SOLVE_TWO_BOUNDARIES_COMPLETE: - if(steady_state) - { - mexPrintf("SOLVE_TWO_BOUNDARIES in a steady state model: impossible case\n"); - return false; - } - fixe_u(&u, u_count_int, u_count_int); - Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true); - u_count=u_count_int*(periods+y_kmax+y_kmin); - r=(double*)mxMalloc(size*sizeof(double)); - y_save=(double*)mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); - begining = it_code; - if(!Gaussian_Elimination) - { - } - giter=0; - iter=0; - if (!is_linear) - { - cvg=false; - int u_count_saved=u_count; - while (!(cvg||(iter>maxit_))) - { - res2=0; - res1=0; - max_res=0; - max_res_idx=0; - memcpy(y_save, y, y_size*sizeof(double)*(periods+y_kmax+y_kmin)); - for (it_=y_kmin;it_eps) - rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); - else - rr=r[i]; - if (max_res= y_kmin; it_--) + { + it_code = begining; + Per_y_ = it_*y_size; + compute_block_time(0, false, block_num); + } + } + break; + case SOLVE_FORWARD_SIMPLE: + g1 = (double *) mxMalloc(size*size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + begining = it_code; + if (steady_state) + { + cvg = false; + iter = 0; + while (!(cvg || (iter > maxit_))) + { + it_code = begining; + Per_y_ = it_*y_size; + compute_block_time(0, false, block_num); + double rr; + rr = r[0]; + cvg = (fabs(rr) < solve_tolf); + if (cvg) + continue; + y[Block_Contain[0].Variable] += -r[0]/g1[0]; + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, after %d iterations\n", Block_Count, iter); + mexPrintf("r[0]= %f\n", r[0]); + return false; + } + } + else + { + for (it_ = y_kmin; it_ < periods+y_kmin; it_++) + { + cvg = false; + iter = 0; + Per_y_ = it_*y_size; + while (!(cvg || (iter > maxit_))) + { + it_code = begining; + Per_y_ = it_*y_size; + compute_block_time(0, false, block_num); + double rr; + if (fabs(1+y[Per_y_+Block_Contain[0].Variable]) > eps) + rr = r[0]/(1+y[Per_y_+Block_Contain[0].Variable]); + else + rr = r[0]; + cvg = (fabs(rr) < solve_tolf); + if (cvg) + continue; + y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0]; + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter); + mexErrMsgTxt("End of bytecode"); + } + } + } + mxFree(g1); + mxFree(r); + break; + case SOLVE_BACKWARD_SIMPLE: + g1 = (double *) mxMalloc(size*size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + begining = it_code; + if (steady_state) + { + cvg = false; + iter = 0; + while (!(cvg || (iter > maxit_))) + { + it_code = begining; + Per_y_ = it_*y_size; + compute_block_time(0, false, block_num); + double rr; + rr = r[0]; + cvg = (fabs(rr) < solve_tolf); + if (cvg) + continue; + y[Block_Contain[0].Variable] += -r[0]/g1[0]; + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, after %d iterations\n", Block_Count, iter); + return false; + } + } + else + { + for (it_ = periods+y_kmin; it_ > y_kmin; it_--) + { + cvg = false; + iter = 0; + Per_y_ = it_*y_size; + while (!(cvg || (iter > maxit_))) + { + it_code = begining; + Per_y_ = it_*y_size; + compute_block_time(0, false, block_num); + double rr; + if (fabs(1+y[Per_y_+Block_Contain[0].Variable]) > eps) + rr = r[0]/(1+y[Per_y_+Block_Contain[0].Variable]); + else + rr = r[0]; + cvg = (fabs(rr) < solve_tolf); + if (cvg) + continue; + y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0]; + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter); + mexErrMsgTxt("End of bytecode"); + } + } + } + mxFree(g1); + mxFree(r); + break; + case SOLVE_FORWARD_COMPLETE: + fixe_u(&u, u_count_int, u_count_int); + Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false); + g1 = (double *) mxMalloc(size*size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + begining = it_code; + Per_u_ = 0; + if (steady_state) + { + if (!is_linear) + { + max_res_idx = 0; + cvg = false; + iter = 0; + while (!(cvg || (iter > maxit_))) + { + it_code = begining; + error_not_printed = true; + res2 = 0; + res1 = 0; + max_res = 0; + compute_block_time(0, false, block_num); + if (!(isnan(res1) || isinf(res1))) + { + for (i = 0; i < size; i++) + { + double rr; + rr = r[i]; + if (max_res < fabs(rr)) + { + max_res = fabs(rr); + max_res_idx = i; + } + res2 += rr*rr; + res1 += fabs(rr); + } + cvg = (max_res < solve_tolf); + } + else + cvg = false; + if (cvg) + continue; + result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true); + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, after %d iterations\n", Block_Count, iter); + return false; + } + } + else + { + it_code = begining; + Per_y_ = it_*y_size; + iter = 0; + res1 = res2 = max_res = 0; max_res_idx = 0; + error_not_printed = true; + compute_block_time(0, false, block_num); + cvg = false; + result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true); + } + } + else + { + if (!is_linear) + { + max_res_idx = 0; + for (it_ = y_kmin; it_ < periods+y_kmin; it_++) + { + cvg = false; + iter = 0; + Per_y_ = it_*y_size; + while (!(cvg || (iter > maxit_))) + { + it_code = begining; + error_not_printed = true; + res2 = 0; + res1 = 0; + max_res = 0; + compute_block_time(0, false, block_num); + if (!(isnan(res1) || isinf(res1))) + { + for (i = 0; i < size; i++) + { + double rr; + if (fabs(1+y[Per_y_+Block_Contain[i].Variable]) > eps) + rr = r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); + else + rr = r[i]; + if (max_res < fabs(rr)) + { + max_res = fabs(rr); + max_res_idx = i; + } + res2 += rr*rr; + res1 += fabs(rr); + } + cvg = (max_res < solve_tolf); + } + else + cvg = false; + if (cvg) + continue; + result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false); + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter); + mexErrMsgTxt("End of bytecode"); + } + } + } + else + { + for (it_ = y_kmin; it_ < periods+y_kmin; it_++) + { + it_code = begining; + Per_y_ = it_*y_size; + iter = 0; + res1 = res2 = max_res = 0; max_res_idx = 0; + error_not_printed = true; + compute_block_time(0, false, block_num); + cvg = false; + result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false); + } + } + } + mxFree(index_equa); + mxFree(index_vara); + memset(direction, 0, size_of_direction); + mxFree(g1); + mxFree(r); + mxFree(u); + break; + case SOLVE_BACKWARD_COMPLETE: + fixe_u(&u, u_count_int, u_count_int); + Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false); + g1 = (double *) mxMalloc(size*size*sizeof(double)); + r = (double *) mxMalloc(size*sizeof(double)); + begining = it_code; + if (steady_state) + { + if (!is_linear) + { + max_res_idx = 0; + cvg = false; + iter = 0; + while (!(cvg || (iter > maxit_))) + { + it_code = begining; + error_not_printed = true; + res2 = 0; + res1 = 0; + max_res = 0; + compute_block_time(0, false, block_num); + if (!(isnan(res1) || isinf(res1))) + { + for (i = 0; i < size; i++) + { + double rr; + rr = r[i]; + if (max_res < fabs(rr)) + { + max_res = fabs(rr); + max_res_idx = i; + } + res2 += rr*rr; + res1 += fabs(rr); + } + cvg = (max_res < solve_tolf); + } + else + cvg = false; + if (cvg) + continue; + result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true); + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, after %d iterations\n", Block_Count, iter); + return false; + } + } + else + { + it_code = begining; + Per_y_ = it_*y_size; + iter = 0; + res1 = res2 = max_res = 0; max_res_idx = 0; + error_not_printed = true; + compute_block_time(0, false, block_num); + cvg = false; + result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true); + } + } + else + { + if (!is_linear) + { + max_res_idx = 0; + for (it_ = periods+y_kmin; it_ > y_kmin; it_--) + { + cvg = false; + iter = 0; + Per_y_ = it_*y_size; + while (!(cvg || (iter > maxit_))) + { + it_code = begining; + error_not_printed = true; + res2 = 0; + res1 = 0; + max_res = 0; + compute_block_time(0, false, block_num); + if (!(isnan(res1) || isinf(res1))) + { + for (i = 0; i < size; i++) + { + double rr; + if (fabs(1+y[Per_y_+Block_Contain[i].Variable]) > eps) + rr = r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); + else + rr = r[i]; + if (max_res < fabs(rr)) + { + max_res = fabs(rr); + max_res_idx = i; + } + res2 += rr*rr; + res1 += fabs(rr); + } + cvg = (max_res < solve_tolf); + } + else + cvg = false; + if (cvg) + continue; + result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false); + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n", Block_Count, it_, iter); + mexErrMsgTxt("End of bytecode"); + } + } + } + else + { + for (it_ = periods+y_kmin; it_ > y_kmin; it_--) + { + it_code = begining; + Per_y_ = it_*y_size; + error_not_printed = true; + compute_block_time(0, false, block_num); + cvg = false; + result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false); + } + } + } + mxFree(index_equa); + mxFree(index_vara); + memset(direction, 0, size_of_direction); + mxFree(g1); + mxFree(r); + mxFree(u); + break; + case SOLVE_TWO_BOUNDARIES_SIMPLE: + case SOLVE_TWO_BOUNDARIES_COMPLETE: + if (steady_state) + { + mexPrintf("SOLVE_TWO_BOUNDARIES in a steady state model: impossible case\n"); + return false; + } + fixe_u(&u, u_count_int, u_count_int); + Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true); + u_count = u_count_int*(periods+y_kmax+y_kmin); + r = (double *) mxMalloc(size*sizeof(double)); + y_save = (double *) mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); + begining = it_code; + if (!Gaussian_Elimination) + { + } + giter = 0; + iter = 0; + if (!is_linear) + { + cvg = false; + int u_count_saved = u_count; + while (!(cvg || (iter > maxit_))) + { + res2 = 0; + res1 = 0; + max_res = 0; + max_res_idx = 0; + memcpy(y_save, y, y_size*sizeof(double)*(periods+y_kmax+y_kmin)); + for (it_ = y_kmin; it_ < periods+y_kmin; it_++) + { + Per_u_ = (it_-y_kmin)*u_count_int; + Per_y_ = it_*y_size; + it_code = begining; + compute_block_time(Per_u_, false, block_num); + if (isnan(res1) || isinf(res1)) + { + memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin)); + break; + } + for (i = 0; i < size; i++) + { + double rr; + if (fabs(1+y[Per_y_+Block_Contain[i].Variable]) > eps) + rr = r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); + else + rr = r[i]; + if (max_res < fabs(rr)) + { + max_res = fabs(rr); + max_res_idx = i; + } + res2 += rr*rr; + res1 += fabs(rr); + } + } + if (isnan(res1) || isinf(res1)) + cvg = false; + else + cvg = (max_res < solve_tolf); + u_count = u_count_saved; + simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods); + iter++; + } + if (!cvg) + { + mexPrintf("Convergence not achieved in block %d, after %d iterations\n", Block_Count, iter); + mexErrMsgTxt("End of bytecode"); + } + } + else + { + res1 = res2 = max_res = 0; max_res_idx = 0; + for (it_ = y_kmin; it_ < periods+y_kmin; it_++) + { + Per_u_ = (it_-y_kmin)*u_count_int; + Per_y_ = it_*y_size; + it_code = begining; + compute_block_time(Per_u_, false, block_num); + for (i = 0; i < size; i++) + { + double rr; + rr = r[i]; + if (max_res < fabs(rr)) + { + max_res = fabs(rr); + max_res_idx = i; + } + res2 += rr*rr; + res1 += fabs(rr); + } + } + cvg = false; + simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods); + } + mxFree(r); + mxFree(y_save); + mxFree(u); + mxFree(index_vara); + mxFree(index_equa); + memset(direction, 0, size_of_direction); + break; + default: + mexPrintf("Unknown type =%d\n", type); + mexEvalString("drawnow;"); + mexErrMsgTxt("End of bytecode"); } - return true; + return true; } bool @@ -1369,89 +1369,89 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s { bool result = true; - + int var; - if(steady_state) + if (steady_state) file_name += "_static"; - else - file_name += "_dynamic"; - CodeLoad code; - //First read and store in memory the code + else + file_name += "_dynamic"; + CodeLoad code; + //First read and store in memory the code code_liste = code.get_op_code(file_name); if (!code_liste.size()) { - mexPrintf("%s.cod Cannot be opened\n",file_name.c_str()); + mexPrintf("%s.cod Cannot be opened\n", file_name.c_str()); mexEvalString("drawnow;"); - filename+=" stopped"; + filename += " stopped"; mexEvalString("drawnow;"); mexErrMsgTxt(filename.c_str()); } //The big loop on intructions - Block_Count=-1; - bool go_on=true; + Block_Count = -1; + bool go_on = true; it_code = code_liste.begin(); - it_code_type Init_Code=it_code; + it_code_type Init_Code = it_code; while (go_on) { switch (it_code->first) { - case FBEGINBLOCK : - Block_Count++; + case FBEGINBLOCK: + Block_Count++; #ifdef DEBUG - mexPrintf("FBEGINBLOCK %d\n",Block_Count+1); + mexPrintf("FBEGINBLOCK %d\n", Block_Count+1); #endif - //it's a new block - { - FBEGINBLOCK_ *fb = (FBEGINBLOCK_*)it_code->second; - Block_Contain = fb->get_Block_Contain(); - it_code++; - if(evaluate) - evaluate_a_block(fb->get_size(), fb->get_type(), bin_basename, steady_state, Block_Count, - fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); - else - result = simulate_a_block(fb->get_size(), fb->get_type(), file_name, bin_basename,true, steady_state, Block_Count, + //it's a new block + { + FBEGINBLOCK_ *fb = (FBEGINBLOCK_ *) it_code->second; + Block_Contain = fb->get_Block_Contain(); + it_code++; + if (evaluate) + evaluate_a_block(fb->get_size(), fb->get_type(), bin_basename, steady_state, Block_Count, + fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); + else + result = simulate_a_block(fb->get_size(), fb->get_type(), file_name, bin_basename, true, steady_state, Block_Count, fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); - delete fb; - } - if(!result) - go_on = false; - break; - case FEND : + delete fb; + } + if (!result) + go_on = false; + break; + case FEND: #ifdef DEBUG - mexPrintf("FEND\n"); + mexPrintf("FEND\n"); #endif - go_on=false; - it_code++; - break; - case FDIMT : + go_on = false; + it_code++; + break; + case FDIMT: #ifdef DEBUG - mexPrintf("FDIMT size=%d\n",((FDIMT_*)it_code->second)->get_size()); + mexPrintf("FDIMT size=%d\n", ((FDIMT_ *) it_code->second)->get_size()); #endif - var = ((FDIMT_*)it_code->second)->get_size(); - if(T) - mxFree(T); - T=(double*)mxMalloc(var*(periods+y_kmin+y_kmax)*sizeof(double)); - it_code++; - break; - case FDIMST : + var = ((FDIMT_ *) it_code->second)->get_size(); + if (T) + mxFree(T); + T = (double *) mxMalloc(var*(periods+y_kmin+y_kmax)*sizeof(double)); + it_code++; + break; + case FDIMST: #ifdef DEBUG - mexPrintf("FDIMST\n"); + mexPrintf("FDIMST\n"); #endif - var = ((FDIMST_*)it_code->second)->get_size(); - if(T) - mxFree(T); - T=(double*)mxMalloc(var*sizeof(double)); - it_code++; - break; - default : - mexPrintf("Unknown command \n"); - mexEvalString("drawnow;"); - mexErrMsgTxt("End of bytecode"); - break; + var = ((FDIMST_ *) it_code->second)->get_size(); + if (T) + mxFree(T); + T = (double *) mxMalloc(var*sizeof(double)); + it_code++; + break; + default: + mexPrintf("Unknown command \n"); + mexEvalString("drawnow;"); + mexErrMsgTxt("End of bytecode"); + break; } } mxFree(Init_Code->second); - if(T) + if (T) mxFree(T); - return result; + return result; } diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh index d45fff7b7..818d1cf08 100644 --- a/mex/sources/bytecode/Interpreter.hh +++ b/mex/sources/bytecode/Interpreter.hh @@ -31,54 +31,52 @@ # include "linbcg.hh" #endif #ifndef DEBUG_EX - #include "mex.h" +# include "mex.h" #else - #include "mex_interface.hh" +# include "mex_interface.hh" #endif //#define DEBUGC using namespace std; - #define pow_ pow -typedef vector >::const_iterator it_code_type; +typedef vector >::const_iterator it_code_type; class Interpreter : SparseMatrix { - protected : - double pow1(double a, double b); - double log1(double a); - void compute_block_time(int Per_u_, bool evaluate, int block_num); - void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num, - const bool is_linear=false, const int symbol_table_endo_nbr=0, const int Block_List_Max_Lag=0, const int Block_List_Max_Lead=0, const int u_count_int=0); - bool simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num, - const bool is_linear=false, const int symbol_table_endo_nbr=0, const int Block_List_Max_Lag=0, const int Block_List_Max_Lead=0, const int u_count_int=0); - double *T; - vector Block_Contain; - vector > code_liste; - it_code_type it_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 solve_tolf; - bool GaussSeidel; - double *x, *params; - double *steady_y, *steady_x; - map ,int>, int> IM_i; - int equation, derivative_equation, derivative_variable; - string filename; - int minimal_solving_periods; - public : +protected: + double pow1(double a, double b); + double log1(double a); + void compute_block_time(int Per_u_, bool evaluate, int block_num); + void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num, + const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0); + bool simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num, + const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0); + double *T; + vector Block_Contain; + vector > code_liste; + it_code_type it_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 solve_tolf; + bool GaussSeidel; + double *x, *params; + double *steady_y, *steady_x; + map, int>, int> IM_i; + int equation, derivative_equation, derivative_variable; + string filename; + int minimal_solving_periods; +public: - Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_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, int minimal_solving_periods_arg); - bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate); + Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_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, int minimal_solving_periods_arg); + bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate); }; - #endif diff --git a/mex/sources/bytecode/Mem_Mngr.cc b/mex/sources/bytecode/Mem_Mngr.cc index 6afbd881d..37162a971 100644 --- a/mex/sources/bytecode/Mem_Mngr.cc +++ b/mex/sources/bytecode/Mem_Mngr.cc @@ -21,16 +21,16 @@ Mem_Mngr::Mem_Mngr() { - swp_f=false; - swp_f_b=0; + swp_f = false; + swp_f_b = 0; } void Mem_Mngr::Print_heap() { int i; mexPrintf("i :"); - for (i=0;i=0)) + gap = ((size_t)(pos)-(size_t)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem); + if ((gap < CHUNK_BLCK_SIZE) && (gap >= 0)) break; } - Chunk_Stack.push_back((NonZeroElem*)pos); + Chunk_Stack.push_back((NonZeroElem *) pos); } - void Mem_Mngr::Free_All() { - while(NZE_Mem_Allocated.size()) - { - mxFree(NZE_Mem_Allocated.back()); - NZE_Mem_Allocated.pop_back(); - } + while (NZE_Mem_Allocated.size()) + { + mxFree(NZE_Mem_Allocated.back()); + NZE_Mem_Allocated.pop_back(); + } mxFree(NZE_Mem_add); init_Mem(); } diff --git a/mex/sources/bytecode/Mem_Mngr.hh b/mex/sources/bytecode/Mem_Mngr.hh index 73bddc584..3cc7d6ffe 100644 --- a/mex/sources/bytecode/Mem_Mngr.hh +++ b/mex/sources/bytecode/Mem_Mngr.hh @@ -20,47 +20,46 @@ #ifndef MEM_MNGR_HH_INCLUDED #define MEM_MNGR_HH_INCLUDED - #include #include #ifndef DEBUG_EX - #include "mex.h" +# include "mex.h" #else - #include "mex_interface.hh" +# include "mex_interface.hh" #endif using namespace std; struct NonZeroElem - { - int u_index; - int r_index, c_index, lag_index; - NonZeroElem *NZE_R_N, *NZE_C_N; - }; +{ + int u_index; + int r_index, c_index, lag_index; + NonZeroElem *NZE_R_N, *NZE_C_N; +}; -typedef vector v_NonZeroElem; +typedef vector v_NonZeroElem; class Mem_Mngr { public: - 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); - bool 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); + bool swp_f; private: - v_NonZeroElem Chunk_Stack; - int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK; - int CHUNK_heap_pos; - NonZeroElem** NZE_Mem_add; - NonZeroElem* NZE_Mem; - vector NZE_Mem_Allocated; - int swp_f_b; - fstream SaveCode_swp; - string filename; + v_NonZeroElem Chunk_Stack; + int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK; + int CHUNK_heap_pos; + NonZeroElem **NZE_Mem_add; + NonZeroElem *NZE_Mem; + vector NZE_Mem_Allocated; + int swp_f_b; + fstream SaveCode_swp; + string filename; }; #endif diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 2e894bdff..7917ebb86 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -141,7 +141,7 @@ SparseMatrix::At_Col(int c, int lag, NonZeroElem **first) void SparseMatrix::Delete(const int r, const int c) { - NonZeroElem *first = FNZE_R[r], *firsta = NULL; + NonZeroElem *first = FNZE_R[r], *firsta = NULL; while (first->c_index != c) { @@ -165,7 +165,7 @@ SparseMatrix::Delete(const int r, const int c) if (firsta != NULL) firsta->NZE_C_N = first->NZE_C_N; if (first == FNZE_C[c]) - FNZE_C[c] = first->NZE_C_N; + FNZE_C[c] = first->NZE_C_N; u_liste.push_back(first->u_index); mem_mngr.mxFree_NZE(first); @@ -283,24 +283,24 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i mem_mngr.fixe_file_name(file_name); if (!SaveCode.is_open()) { - if(steady_state) + if (steady_state) SaveCode.open((file_name + "_static.bin").c_str(), ios::in | ios::binary); - else - SaveCode.open((file_name + "_dynamic.bin").c_str(), ios::in | ios::binary); + else + SaveCode.open((file_name + "_dynamic.bin").c_str(), ios::in | ios::binary); if (!SaveCode.is_open()) { - if(steady_state) + if (steady_state) mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_static.bin").c_str()); - else - mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_dynamic.bin").c_str()); + else + mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_dynamic.bin").c_str()); mexEvalString("st=fclose('all');clear all;"); mexErrMsgTxt("Exit from Dynare"); } } IM_i.clear(); - if(two_boundaries) - { - for (i = 0; i < u_count_init-Size; i++) + if (two_boundaries) + { + for (i = 0; i < u_count_init-Size; i++) { SaveCode.read(reinterpret_cast(&eq), sizeof(eq)); SaveCode.read(reinterpret_cast(&var), sizeof(var)); @@ -308,12 +308,12 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i SaveCode.read(reinterpret_cast(&j), sizeof(j)); IM_i[make_pair(make_pair(eq, var), lag)] = j; } - for (j=0;j(&eq), sizeof(eq)); SaveCode.read(reinterpret_cast(&var), sizeof(var)); @@ -321,7 +321,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i SaveCode.read(reinterpret_cast(&j), sizeof(j)); IM_i[make_pair(make_pair(eq, var), lag)] = j; } - } + } index_vara = (int *) mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int)); for (j = 0; j < Size; j++) SaveCode.read(reinterpret_cast(&index_vara[j]), sizeof(*index_vara)); @@ -415,7 +415,7 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, mapfirst+t*diff1[j]]; switch (save_op_s->operat) { @@ -751,7 +751,7 @@ SparseMatrix::compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, in i = j = 0; while (i < nop4) { - save_op_s = (t_save_op_s *) (&(save_op[i])); + save_op_s = (t_save_op_s *)(&(save_op[i])); if (save_op_s->lag < (periods_beg_t-t)) { up = &u[save_op_s->first+t*diff1[j]]; @@ -1018,14 +1018,14 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, double piv_abs; NonZeroElem *first, *firsta, *first_suba; double *piv_v; - int *pivj_v, *pivk_v, *NR; + int *pivj_v, *pivk_v, *NR; int l, N_max; bool one; - Clear_u(); - piv_v = (double*)mxMalloc(Size*sizeof(double)); - pivj_v = (int*)mxMalloc(Size*sizeof(int)); - pivk_v = (int*)mxMalloc(Size*sizeof(int)); - NR = (int*)mxMalloc(Size*sizeof(int)); + Clear_u(); + piv_v = (double *) mxMalloc(Size*sizeof(double)); + pivj_v = (int *) mxMalloc(Size*sizeof(int)); + pivk_v = (int *) mxMalloc(Size*sizeof(int)); + NR = (int *) mxMalloc(Size*sizeof(int)); error_not_printed = true; u_count_alloc_save = u_count_alloc; if (isnan(res1) || isinf(res1)) @@ -1034,63 +1034,63 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, { for (j = 0; j < y_size; j++) { - bool select=false; - for(int i = 0; i variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - else - mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); + else + mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); } - mexPrintf("res1=%5.25\n",res1); + mexPrintf("res1=%5.25\n", res1); mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); mexPrintf("Change them!\n"); - mexEvalString("drawnow;"); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); + mexEvalString("drawnow;"); + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); mxFree(NR); - if(steady_state) + if (steady_state) return false; - else - { + else + { mexEvalString("st=fclose('all');clear all;"); filename += " stopped"; mexErrMsgTxt(filename.c_str()); - } + } } if (slowc_save < 1e-8) { for (j = 0; j < y_size; j++) { - bool select=false; - for(int i = 0; i variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - else - mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); + else + mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); } mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx); mexEvalString("drawnow;"); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - if(steady_state) - return false; + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); + if (steady_state) + return false; else { mexEvalString("st=fclose('all');clear all;"); - filename += " stopped"; + filename += " stopped"; mexErrMsgTxt(filename.c_str()); } } @@ -1098,24 +1098,24 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save); for (i = 0; i < y_size; i++) y[i+it_*y_size] = ya[i+it_*y_size] + slowc_save*direction[i+it_*y_size]; - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); iter--; return true; } - if (cvg) - { - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - return (true); + if (cvg) + { + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); + return (true); } - Simple_Init(it_, y_kmin, y_kmax, Size, IM_i); - NonZeroElem **bc; - bc = (NonZeroElem**)mxMalloc(Size*sizeof(*bc)); + Simple_Init(it_, y_kmin, y_kmax, Size, IM_i); + NonZeroElem **bc; + bc = (NonZeroElem **) mxMalloc(Size*sizeof(*bc)); for (i = 0; i < Size; i++) { /*finding the max-pivot*/ @@ -1183,7 +1183,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, { for (j = 0; j < l; j++) { - markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double(N_max))); + markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max))); if (markovitz > markovitz_max && NR[j] == 1) { piv = piv_v[j]; @@ -1199,14 +1199,14 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, line_done[pivj] = true; if (piv_abs < eps) { - mexPrintf("Error: singular system in Simulate_NG in block %d\n",blck+1); - mexEvalString("drawnow;"); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); + mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1); + mexEvalString("drawnow;"); + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); mxFree(bc); - if(steady_state) + if (steady_state) return false; else { @@ -1236,7 +1236,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, } //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) for (j = 0; j < nb_eq_todo; j++) - { + { first = bc[j]; int row = first->r_index; double first_elem = u[first->u_index]; @@ -1310,18 +1310,18 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, } } u[b[row]] -= u[b[pivj]]*first_elem; - } + } } double slowc_lbx = slowc, res1bx; for (i = 0; i < y_size; i++) ya[i+it_*y_size] = y[i+it_*y_size]; slowc_save = slowc; res1bx = simple_bksub(it_, Size, slowc_lbx); - End(Size); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); + End(Size); + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); mxFree(bc); return true; } @@ -1382,8 +1382,8 @@ SparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, mexPrintf("G1a red done\n"); SaveResult >> row; mexPrintf("row(2)=%d\n", row); - double *B; - B = (double*)mxMalloc(row*sizeof(double)); + double *B; + B = (double *) mxMalloc(row*sizeof(double)); for (int i = 0; i < row; i++) SaveResult >> B[i]; SaveResult.close(); @@ -1393,7 +1393,7 @@ SparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, { if (abs(u[b[i]]+B[i]) > epsilon) mexPrintf("Problem at i=%d u[b[i]]=%f B[i]=%f\n", i, u[b[i]], B[i]); - } + } mxFree(B); } @@ -1470,9 +1470,9 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax { for (j = 0; j < y_size; j++) mexPrintf("variable %d at time %d = %f\n", j+1, it_, y[j+it_*y_size]); - for(j = 0; j < Size; j++) - mexPrintf("residual(%d)=%5.25f\n",j, u[j]); - mexPrintf("res1=%5.25f\n",res1); + for (j = 0; j < Size; j++) + mexPrintf("residual(%d)=%5.25f\n", j, u[j]); + mexPrintf("res1=%5.25f\n", res1); mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); mexPrintf("Change them!\n"); mexEvalString("drawnow;"); @@ -1526,7 +1526,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax } else { - start_compare = max( y_kmin, minimal_solving_periods); + start_compare = max(y_kmin, minimal_solving_periods); restart = 0; } res1a = res1; @@ -1549,10 +1549,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax Init(periods, y_kmin, y_kmax, Size, IM_i); double *piv_v; int *pivj_v, *pivk_v, *NR; - piv_v = (double*)mxMalloc(Size*sizeof(double)); - pivj_v = (int*)mxMalloc(Size*sizeof(int)); - pivk_v = (int*)mxMalloc(Size*sizeof(int)); - NR = (int*)mxMalloc(Size*sizeof(int)); + piv_v = (double *) mxMalloc(Size*sizeof(double)); + pivj_v = (int *) mxMalloc(Size*sizeof(int)); + pivk_v = (int *) mxMalloc(Size*sizeof(int)); + NR = (int *) mxMalloc(Size*sizeof(int)); for (int t = 0; t < periods; t++) { if (record && symbolic) @@ -1670,7 +1670,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax nopa = int (1.5*nopa); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); } - save_op_s = (t_save_op_s *) (&(save_op[nop])); + 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; @@ -1687,7 +1687,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax /*divide all the non zeros elements of the line pivj by the max_pivot*/ int nb_var = At_Row(pivj, &first); NonZeroElem **bb; - bb = (NonZeroElem**)mxMalloc(nb_var*sizeof(first)); + bb = (NonZeroElem **) mxMalloc(nb_var*sizeof(first)); for (j = 0; j < nb_var; j++) { bb[j] = first; @@ -1707,13 +1707,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax nopa = int (1.5*nopa); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); } - save_op_s = (t_save_op_s *) (&(save_op[nop+j*2])); + save_op_s = (t_save_op_s *)(&(save_op[nop+j*2])); save_op_s->operat = IFDIV; save_op_s->first = first->u_index; save_op_s->lag = first->lag_index; } } - } + } mxFree(bb); nop += nb_var*2; u[b[pivj]] /= piv; @@ -1726,7 +1726,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax nopa = int (1.5*nopa); save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); } - save_op_s = (t_save_op_s *) (&(save_op[nop])); + 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; @@ -1739,7 +1739,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax int nb_var_piva = At_Row(pivj, &first_piva); NonZeroElem **bc; - bc = (NonZeroElem**)mxMalloc(nb_eq*sizeof(first)); + bc = (NonZeroElem **) mxMalloc(nb_eq*sizeof(first)); int nb_eq_todo = 0; for (j = 0; j < nb_eq && first; j++) { @@ -1766,7 +1766,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); } } - save_op_s_l = (t_save_op_s *) (&(save_op[nop])); + save_op_s_l = (t_save_op_s *)(&(save_op[nop])); save_op_s_l->operat = IFLD; save_op_s_l->first = first->u_index; save_op_s_l->lag = abs(first->lag_index); @@ -1817,7 +1817,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); } } - save_op_s_l = (t_save_op_s *) (&(save_op[nop])); + save_op_s_l = (t_save_op_s *)(&(save_op[nop])); save_op_s_l->operat = IFLESS; save_op_s_l->first = tmp_u_count; save_op_s_l->second = first_piv->u_index; @@ -1875,7 +1875,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); } } - save_op_s_l = (t_save_op_s *) (&(save_op[nop])); + save_op_s_l = (t_save_op_s *)(&(save_op[nop])); save_op_s_l->operat = IFSUB; save_op_s_l->first = first_sub->u_index; save_op_s_l->second = first_piv->u_index; @@ -1912,7 +1912,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax save_op = (int *) mxRealloc(save_op, nopa*sizeof(int)); } } - save_op_s_l = (t_save_op_s *) (&(save_op[nop])); + save_op_s_l = (t_save_op_s *)(&(save_op[nop])); save_op_s_l->operat = IFSUB; save_op_s_l->first = b[row]; save_op_s_l->second = b[pivj]; @@ -1920,7 +1920,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax } nop += 3; } - } + } mxFree(bc); } if (symbolic) @@ -1977,10 +1977,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax nop1 = nop; } } - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); } nop_all += nop; if (symbolic) diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index e6b2ef77d..cc93a568f 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/mex/sources/bytecode/SparseMatrix.hh @@ -37,35 +37,40 @@ using namespace std; extern unsigned long _nan[2]; extern double NAN; -inline bool isnan(double value) +inline bool +isnan(double value) { return _isnan(value); } -inline bool isinf(double value) +inline bool +isinf(double value) { - return (std::numeric_limits::has_infinity && - value == std::numeric_limits::infinity()); + return (std::numeric_limits::has_infinity + && value == std::numeric_limits::infinity()); } template -inline T asinh(T x) +inline T +asinh(T x) { return log(x+sqrt(x*x+1)); } template -inline T acosh(T x) +inline T +acosh(T x) { - if (!(x>=1.0)) + if (!(x >= 1.0)) return sqrt(-1.0); return log(x+sqrt(x*x-1.0)); } template -inline T atanh(T x) +inline T +atanh(T x) { - if(!(x>-1.0 && x<1.0)) + if (!(x > -1.0 && x < 1.0)) return sqrt(-1.0); return log((1.0+x)/(1.0-x))/2.0; } @@ -78,108 +83,105 @@ struct t_save_op_s 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-10; -const double very_big=1e24; -const int alt_symbolic_count_max=1; - +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-10; +const double very_big = 1e24; +const int alt_symbolic_count_max = 1; class SparseMatrix - { - public: - SparseMatrix(); - int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods); - bool simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state); - 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 Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries); - 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); +{ +public: + SparseMatrix(); + int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods); + bool simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state); + 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 Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries); + 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 Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map ,int>, int> &IM); - void End(int Size); - bool compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size +private: + void Init(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM); + void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM); + void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map, int>, int> &IM); + void End(int Size); + bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size #ifdef PROFILER - , long int *ndiv, long int *nsub + , long int *ndiv, long int *nsub #endif - ); - void Insert(const int r, const int c, const int u_index, const int lag_index); - void Delete(const int r,const int c); - 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(); - void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter); - void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int* pivot, int* b); - int complete(int beg_t, int Size, int periods, int *b); - double bksub( int tbreak, int last_period, int Size, double slowc_l + ); + void Insert(const int r, const int c, const int u_index, const int lag_index); + void Delete(const int r, const int c); + 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(); + void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter); + void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b); + 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 + , long int *nmul #endif - ); - double simple_bksub(int it_, int Size, double slowc_l); - 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; - fstream SaveCode; - string filename; - int max_u, min_u; - clock_t time00; + ); + double simple_bksub(int it_, int Size, double slowc_l); + 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; + fstream SaveCode; + string filename; + int max_u, min_u; + clock_t time00; - Mem_Mngr mem_mngr; - vector u_liste; - map,NonZeroElem*> Mapped_Array; - int *NbNZRow, *NbNZCol; - NonZeroElem **FNZE_R, **FNZE_C; - int nb_endo, u_count_init; + Mem_Mngr mem_mngr; + vector u_liste; + map, NonZeroElem *> Mapped_Array; + int *NbNZRow, *NbNZCol; + NonZeroElem **FNZE_R, **FNZE_C; + int nb_endo, u_count_init; - int *pivot, *pivotk, *pivot_save; - 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, nop1, nop2; - map ,int>, int> IM_i; + int *pivot, *pivotk, *pivot_save; + 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, nop1, nop2; + map, int>, int> IM_i; protected: - double *u, *y, *ya; - double res1, res2, max_res, max_res_idx; - double slowc, slowc_save, markowitz_c; - int y_kmin, y_kmax, y_size, periods, y_decal; - int *index_vara, *index_equa; - int u_count, tbreak_g; - int iter; - double *direction; - int start_compare; - int restart; - bool error_not_printed; - }; - - + double *u, *y, *ya; + double res1, res2, max_res, max_res_idx; + double slowc, slowc_save, markowitz_c; + int y_kmin, y_kmax, y_size, periods, y_decal; + int *index_vara, *index_equa; + int u_count, tbreak_g; + int iter; + double *direction; + int start_compare; + int restart; + bool error_not_printed; +}; #endif diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 43d7b6a07..40d2b8e41 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -26,18 +26,17 @@ #include "Interpreter.hh" #ifndef DEBUG_EX - #include "mex.h" +# include "mex.h" #else - #include "mex_interface.hh" +# include "mex_interface.hh" #endif #include "Mem_Mngr.hh" - #ifdef DEBUG_EX using namespace std; -#include +# include string Get_Argument(const char *argv) { @@ -45,23 +44,22 @@ Get_Argument(const char *argv) return f; } - int -main( int argc, const char* argv[] ) +main(int argc, const char *argv[]) { FILE *fid; bool steady_state = false; bool evaluate = false; - printf("argc=%d\n",argc); - if(argc<2) + printf("argc=%d\n", argc); + if (argc < 2) { - mexPrintf("model filename expected\n"); - mexEvalString("st=fclose('all');clear all;"); + mexPrintf("model filename expected\n"); + mexEvalString("st=fclose('all');clear all;"); mexErrMsgTxt("Exit from Dynare"); } float f_tmp; ostringstream tmp_out(""); - tmp_out << argv[1] << "_options.txt"; + tmp_out << argv[1] << "_options.txt"; cout << tmp_out.str().c_str() << "\n"; int nb_params; int i, row_y, col_y, row_x, col_x; @@ -71,13 +69,13 @@ main( int argc, const char* argv[] ) string file_name(argv[1]); - for(i=2;i0) + clock_t t1 = clock(); + if (!steady_state && !evaluate) + mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC)); + if (nlhs > 0) { plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL); pind = mxGetPr(plhs[0]); - if(evaluate) - for (i=0;i1) - { - plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); - pind = mxGetPr(plhs[1]); - if(result) - pind[0] = 0; - else - pind[0] = 1; - } + for (i = 0; i < row_y*col_y; i++) + pind[i] = y[i]; + if (nlhs > 1) + { + plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); + pind = mxGetPr(plhs[1]); + if (result) + pind[0] = 0; + else + pind[0] = 1; + } } - if(x) + if (x) mxFree(x); - if(y) + if (y) mxFree(y); - if(ya) + if (ya) mxFree(ya); - if(direction) + if (direction) mxFree(direction); } #endif diff --git a/mex/sources/bytecode/testing/mex_interface.cc b/mex/sources/bytecode/testing/mex_interface.cc index 40b76f0f0..60085ef44 100644 --- a/mex/sources/bytecode/testing/mex_interface.cc +++ b/mex/sources/bytecode/testing/mex_interface.cc @@ -10,9 +10,9 @@ mexPrintf(const char *str, ...) va_list args; int retval; - va_start (args, str); - retval = vprintf (str, args); - va_end (args); + va_start(args, str); + retval = vprintf(str, args); + va_end(args); return retval; } @@ -25,24 +25,23 @@ mexErrMsgTxt(const string str) } void -mxFree(void* to_release) +mxFree(void *to_release) { free(to_release); } -void* +void * mxMalloc(int amount) { return malloc(amount); } -void* -mxRealloc(void* to_extend, int amount) +void * +mxRealloc(void *to_extend, int amount) { return realloc(to_extend, amount); } - void mexEvalString(const string str) { diff --git a/mex/sources/bytecode/testing/mex_interface.hh b/mex/sources/bytecode/testing/mex_interface.hh index 3999121db..cb26b9877 100644 --- a/mex/sources/bytecode/testing/mex_interface.hh +++ b/mex/sources/bytecode/testing/mex_interface.hh @@ -5,9 +5,9 @@ #include using namespace std; -int mexPrintf(/*const string*/const char* str, ...); +int mexPrintf(/*const string*/ const char *str, ...); void mexErrMsgTxt(const string str); -void* mxMalloc(int amount); -void* mxRealloc(void* to_extend, int amount); -void mxFree(void* to_release); +void *mxMalloc(int amount); +void *mxRealloc(void *to_extend, int amount); +void mxFree(void *to_release); void mexEvalString(const string str); diff --git a/mex/sources/bytecode/testing/simulate_debug.m b/mex/sources/bytecode/testing/simulate_debug.m index 28e9ab2c8..1550585d8 100644 --- a/mex/sources/bytecode/testing/simulate_debug.m +++ b/mex/sources/bytecode/testing/simulate_debug.m @@ -1,37 +1,36 @@ -function simulate_debug() -global M_ oo_ options_; - fid = fopen([M_.fname '_options.txt'],'wt'); - fprintf(fid,'%d\n',options_.periods); - fprintf(fid,'%d\n',options_.maxit_); - fprintf(fid,'%6.20f\n',options_.slowc); - fprintf(fid,'%6.20f\n',options_.markowitz); - fprintf(fid,'%6.20f\n',options_.dynatol); - fprintf(fid,'%d\n',options_.minimal_solving_periods); - fclose(fid); - - fid = fopen([M_.fname '_M.txt'],'wt'); - fprintf(fid,'%d\n',M_.maximum_lag); - fprintf(fid,'%d\n',M_.maximum_lead); - fprintf(fid,'%d\n',M_.maximum_endo_lag); - fprintf(fid,'%d\n',M_.param_nbr); - fprintf(fid,'%d\n',size(oo_.exo_simul, 1)); - fprintf(fid,'%d\n',size(oo_.exo_simul, 2)); - fprintf(fid,'%d\n',M_.endo_nbr); - fprintf(fid,'%d\n',size(oo_.endo_simul, 2)); - fprintf(fid,'%d\n',M_.exo_det_nbr); - - fprintf(fid,'%d\n',size(oo_.steady_state,1)); - fprintf(fid,'%d\n',size(oo_.steady_state,2)); - fprintf(fid,'%d\n',size(oo_.exo_steady_state,1)); - fprintf(fid,'%d\n',size(oo_.exo_steady_state,2)); - - fprintf(fid,'%6.20f\n',M_.params); - - fclose(fid); - fid = fopen([M_.fname '_oo.txt'],'wt'); - fprintf(fid,'%6.20f\n',oo_.endo_simul); - fprintf(fid,'%6.20f\n',oo_.exo_simul); - fprintf(fid,'%6.20f\n',oo_.steady_state); - fprintf(fid,'%6.20f\n',oo_.exo_steady_state); - fclose(fid); - \ No newline at end of file +function simulate_debug() +global M_ oo_ options_; +fid = fopen([M_.fname '_options.txt'],'wt'); +fprintf(fid,'%d\n',options_.periods); +fprintf(fid,'%d\n',options_.maxit_); +fprintf(fid,'%6.20f\n',options_.slowc); +fprintf(fid,'%6.20f\n',options_.markowitz); +fprintf(fid,'%6.20f\n',options_.dynatol); +fprintf(fid,'%d\n',options_.minimal_solving_periods); +fclose(fid); + +fid = fopen([M_.fname '_M.txt'],'wt'); +fprintf(fid,'%d\n',M_.maximum_lag); +fprintf(fid,'%d\n',M_.maximum_lead); +fprintf(fid,'%d\n',M_.maximum_endo_lag); +fprintf(fid,'%d\n',M_.param_nbr); +fprintf(fid,'%d\n',size(oo_.exo_simul, 1)); +fprintf(fid,'%d\n',size(oo_.exo_simul, 2)); +fprintf(fid,'%d\n',M_.endo_nbr); +fprintf(fid,'%d\n',size(oo_.endo_simul, 2)); +fprintf(fid,'%d\n',M_.exo_det_nbr); + +fprintf(fid,'%d\n',size(oo_.steady_state,1)); +fprintf(fid,'%d\n',size(oo_.steady_state,2)); +fprintf(fid,'%d\n',size(oo_.exo_steady_state,1)); +fprintf(fid,'%d\n',size(oo_.exo_steady_state,2)); + +fprintf(fid,'%6.20f\n',M_.params); + +fclose(fid); +fid = fopen([M_.fname '_oo.txt'],'wt'); +fprintf(fid,'%6.20f\n',oo_.endo_simul); +fprintf(fid,'%6.20f\n',oo_.exo_simul); +fprintf(fid,'%6.20f\n',oo_.steady_state); +fprintf(fid,'%6.20f\n',oo_.exo_steady_state); +fclose(fid); diff --git a/mex/sources/dynblas.h b/mex/sources/dynblas.h index bdf878dac..15820e421 100644 --- a/mex/sources/dynblas.h +++ b/mex/sources/dynblas.h @@ -30,77 +30,77 @@ #define _DYNBLAS_H /* Starting from version 7.8, MATLAB BLAS expects ptrdiff_t arguments for integers */ -# if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708 -# ifdef __cplusplus -# include -# else -# include -# endif -typedef ptrdiff_t blas_int; -# else -typedef int blas_int; -# endif - -# if defined(MATLAB_MEX_FILE) && defined(_WIN32) -# define FORTRAN_WRAPPER(x) x -# else -# define FORTRAN_WRAPPER(x) x ## _ -# endif - +#if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708 # ifdef __cplusplus -extern "C" { +# include +# else +# include # endif +typedef ptrdiff_t blas_int; +#else +typedef int blas_int; +#endif + +#if defined(MATLAB_MEX_FILE) && defined(_WIN32) +# define FORTRAN_WRAPPER(x) x +#else +# define FORTRAN_WRAPPER(x) x ## _ +#endif + +#ifdef __cplusplus +extern "C" { +#endif typedef const char *BLCHAR; typedef const blas_int *CONST_BLINT; typedef const double *CONST_BLDOU; typedef double *BLDOU; -# define dgemm FORTRAN_WRAPPER(dgemm) +#define dgemm FORTRAN_WRAPPER(dgemm) void dgemm(BLCHAR transa, BLCHAR transb, CONST_BLINT m, CONST_BLINT n, CONST_BLINT k, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta, BLDOU c, CONST_BLINT ldc); -# define dgemv FORTRAN_WRAPPER(dgemv) +#define dgemv FORTRAN_WRAPPER(dgemv) void dgemv(BLCHAR trans, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU beta, BLDOU y, CONST_BLINT incy); -# define dtrsv FORTRAN_WRAPPER(dtrsv) +#define dtrsv FORTRAN_WRAPPER(dtrsv) void dtrsv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n, CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx); -# define dtrmv FORTRAN_WRAPPER(dtrmv) +#define dtrmv FORTRAN_WRAPPER(dtrmv) void dtrmv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n, CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx); -# define daxpy FORTRAN_WRAPPER(daxpy) +#define daxpy FORTRAN_WRAPPER(daxpy) void daxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx, BLDOU y, CONST_BLINT incy); -# define dcopy FORTRAN_WRAPPER(dcopy) +#define dcopy FORTRAN_WRAPPER(dcopy) void dcopy(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, BLDOU y, CONST_BLINT incy); -# define zaxpy FORTRAN_WRAPPER(zaxpy) +#define zaxpy FORTRAN_WRAPPER(zaxpy) void zaxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx, BLDOU y, CONST_BLINT incy); -# define dscal FORTRAN_WRAPPER(dscal) +#define dscal FORTRAN_WRAPPER(dscal) void dscal(CONST_BLINT n, CONST_BLDOU a, BLDOU x, CONST_BLINT incx); -# define dtrsm FORTRAN_WRAPPER(dtrsm) +#define dtrsm FORTRAN_WRAPPER(dtrsm) void dtrsm(BLCHAR side, BLCHAR uplo, BLCHAR transa, BLCHAR diag, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda, BLDOU b, CONST_BLINT ldb); -# define ddot FORTRAN_WRAPPER(ddot) +#define ddot FORTRAN_WRAPPER(ddot) double ddot(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU y, CONST_BLINT incy); -# ifdef __cplusplus +#ifdef __cplusplus } /* extern "C" */ -# endif +#endif #endif /* _DYNBLAS_H */ diff --git a/mex/sources/dynlapack.h b/mex/sources/dynlapack.h index 50ea84ed0..7bed4b314 100644 --- a/mex/sources/dynlapack.h +++ b/mex/sources/dynlapack.h @@ -6,7 +6,7 @@ * * When used in the context of a MATLAB MEX file, you must define MATLAB_MEX_FILE * and MATLAB_VERSION (for version 7.4, define it to 0x0704). - * + * * * Copyright (C) 2009 Dynare Team * @@ -30,88 +30,88 @@ #define _DYNLAPACK_H /* Starting from version 7.8, MATLAB LAPACK expects ptrdiff_t arguments for integers */ -# if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708 -# ifdef __cplusplus -# include -# else -# include -# endif -typedef ptrdiff_t lapack_int; -# else -typedef int lapack_int; -# endif - -# if defined(MATLAB_MEX_FILE) && defined(_WIN32) -# define FORTRAN_WRAPPER(x) x -# else -# define FORTRAN_WRAPPER(x) x ## _ -# endif - +#if defined(MATLAB_MEX_FILE) && MATLAB_VERSION >= 0x0708 # ifdef __cplusplus -extern "C" { +# include +# else +# include # endif +typedef ptrdiff_t lapack_int; +#else +typedef int lapack_int; +#endif + +#if defined(MATLAB_MEX_FILE) && defined(_WIN32) +# define FORTRAN_WRAPPER(x) x +#else +# define FORTRAN_WRAPPER(x) x ## _ +#endif + +#ifdef __cplusplus +extern "C" { +#endif typedef const char *LACHAR; typedef const lapack_int *CONST_LAINT; typedef lapack_int *LAINT; typedef const double *CONST_LADOU; typedef double *LADOU; - typedef lapack_int (*DGGESCRIT)(const double*, const double*, const double*); + typedef lapack_int (*DGGESCRIT)(const double *, const double *, const double *); -# define dgetrs FORTRAN_WRAPPER(dgetrs) +#define dgetrs FORTRAN_WRAPPER(dgetrs) void dgetrs(LACHAR trans, CONST_LAINT n, CONST_LAINT nrhs, CONST_LADOU a, CONST_LAINT lda, CONST_LAINT ipiv, LADOU b, CONST_LAINT ldb, LAINT info); -# define dgetrf FORTRAN_WRAPPER(dgetrf) +#define dgetrf FORTRAN_WRAPPER(dgetrf) void dgetrf(CONST_LAINT m, CONST_LAINT n, LADOU a, CONST_LAINT lda, LAINT ipiv, LAINT info); -# define dgees FORTRAN_WRAPPER(dgees) - void dgees(LACHAR jobvs, LACHAR sort, const void* select, +#define dgees FORTRAN_WRAPPER(dgees) + void dgees(LACHAR jobvs, LACHAR sort, const void *select, CONST_LAINT n, LADOU a, CONST_LAINT lda, LAINT sdim, LADOU wr, LADOU wi, LADOU vs, CONST_LAINT ldvs, - LADOU work, CONST_LAINT lwork, const void* bwork, LAINT info); + LADOU work, CONST_LAINT lwork, const void *bwork, LAINT info); -# define dgecon FORTRAN_WRAPPER(dgecon) +#define dgecon FORTRAN_WRAPPER(dgecon) void dgecon(LACHAR norm, CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda, CONST_LADOU anorm, LADOU rnorm, LADOU work, LAINT iwork, LAINT info); -# define dtrexc FORTRAN_WRAPPER(dtrexc) +#define dtrexc FORTRAN_WRAPPER(dtrexc) void dtrexc(LACHAR compq, CONST_LAINT n, LADOU t, CONST_LAINT ldt, LADOU q, CONST_LAINT ldq, LAINT ifst, LAINT ilst, LADOU work, LAINT info); -# define dtrsyl FORTRAN_WRAPPER(dtrsyl) +#define dtrsyl FORTRAN_WRAPPER(dtrsyl) void dtrsyl(LACHAR trana, LACHAR tranb, CONST_LAINT isgn, CONST_LAINT m, CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda, CONST_LADOU b, CONST_LAINT ldb, LADOU c, CONST_LAINT ldc, LADOU scale, LAINT info); -# define dpotrf FORTRAN_WRAPPER(dpotrf) +#define dpotrf FORTRAN_WRAPPER(dpotrf) void dpotrf(LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda, LAINT info); -# define dgges FORTRAN_WRAPPER(dgges) +#define dgges FORTRAN_WRAPPER(dgges) void dgges(LACHAR jobvsl, LACHAR jobvsr, LACHAR sort, DGGESCRIT delztg, CONST_LAINT n, LADOU a, CONST_LAINT lda, LADOU b, CONST_LAINT ldb, LAINT sdim, LADOU alphar, LADOU alphai, LADOU beta, LADOU vsl, CONST_LAINT ldvsl, LADOU vsr, CONST_LAINT ldvsr, LADOU work, CONST_LAINT lwork, LAINT bwork, LAINT info); -# define dsyev FORTRAN_WRAPPER(dsyev) +#define dsyev FORTRAN_WRAPPER(dsyev) void dsyev(LACHAR jobz, LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda, - LADOU w, LADOU work, CONST_LAINT lwork, LAINT info); + LADOU w, LADOU work, CONST_LAINT lwork, LAINT info); -# define dsyevr FORTRAN_WRAPPER(dsyevr) +#define dsyevr FORTRAN_WRAPPER(dsyevr) void dsyevr(LACHAR jobz, LACHAR range, LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda, LADOU lv, LADOU vu, CONST_LAINT il, CONST_LAINT iu, CONST_LADOU abstol, LAINT m, LADOU w, LADOU z, CONST_LAINT ldz, LAINT isuppz, LADOU work, CONST_LAINT lwork, LAINT iwork, CONST_LAINT liwork, LAINT info); -# ifdef __cplusplus +#ifdef __cplusplus } /* extern "C" */ -# endif +#endif #endif /* _DYNLAPACK_H */ diff --git a/mex/sources/dynmex.h b/mex/sources/dynmex.h index fca6b9feb..c75a1a9f5 100644 --- a/mex/sources/dynmex.h +++ b/mex/sources/dynmex.h @@ -20,16 +20,16 @@ #ifndef _DYNMEX_H #define _DYNMEX_H -# if !defined(MATLAB_MEX_FILE) && !defined(OCTAVE_MEX_FILE) -# error You must define either MATLAB_MEX_FILE or OCTAVE_MEX_FILE -# endif +#if !defined(MATLAB_MEX_FILE) && !defined(OCTAVE_MEX_FILE) +# error You must define either MATLAB_MEX_FILE or OCTAVE_MEX_FILE +#endif -# include +#include /* mwSize, mwIndex and mwSignedIndex appeared in MATLAB 7.3 */ -# if defined(MATLAB_MEX_FILE) && MATLAB_VERSION < 0x0703 +#if defined(MATLAB_MEX_FILE) && MATLAB_VERSION < 0x0703 typedef int mwIndex; typedef int mwSize; -# endif +#endif #endif diff --git a/mex/sources/k_order_perturbation/dynamic_dll.cpp b/mex/sources/k_order_perturbation/dynamic_dll.cpp index 781225f2e..414b9d50a 100644 --- a/mex/sources/k_order_perturbation/dynamic_dll.cpp +++ b/mex/sources/k_order_perturbation/dynamic_dll.cpp @@ -27,8 +27,8 @@ * _dynamic () function **************************************/ DynamicModelDLL::DynamicModelDLL(const string &modName, const int y_length, const int j_cols, - const int n_max_lag, const int n_exog, const string &sExt) throw (DynareException) - : length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog) + const int n_max_lag, const int n_exog, const string &sExt) throw (DynareException) : + length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog) { string fName; #if !defined(__CYGWIN32__) && !defined(_WIN32) @@ -42,12 +42,12 @@ DynamicModelDLL::DynamicModelDLL(const string &modName, const int y_length, cons dynamicHinstance = LoadLibrary(fName.c_str()); if (dynamicHinstance == NULL) throw 1; - Dynamic = (DynamicFn)GetProcAddress(dynamicHinstance, "Dynamic"); + Dynamic = (DynamicFn) GetProcAddress(dynamicHinstance, "Dynamic"); if (Dynamic == NULL) - { - FreeLibrary(dynamicHinstance); // Free the library - throw 2; - } + { + FreeLibrary(dynamicHinstance); // Free the library + throw 2; + } #else // Linux or Mac dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW); if ((dynamicHinstance == NULL) || dlerror()) @@ -58,7 +58,7 @@ DynamicModelDLL::DynamicModelDLL(const string &modName, const int y_length, cons Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic"); if ((Dynamic == NULL) || dlerror()) { - dlclose(dynamicHinstance); // Free the library + dlclose(dynamicHinstance); // Free the library cerr << dlerror() << endl; throw 2; } @@ -120,16 +120,16 @@ DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modPa dg1 = const_cast(g1->base()); } if (g2 != NULL) - dg2 = const_cast(g2->base()); + dg2 = const_cast(g2->base()); dresidual = const_cast(residual.base()); if (g3 != NULL) - dg3 = const_cast(g3->base()); + dg3 = const_cast(g3->base()); dresidual = const_cast(residual.base()); double *dy = const_cast(y.base()); double *dx = const_cast(x.base()); double *dbParams = const_cast(modParams->base()); - Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2, dg3); + Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2, dg3); } void diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cpp b/mex/sources/k_order_perturbation/k_ord_dynare.cpp index e8bdb3c58..20585ae0e 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.cpp +++ b/mex/sources/k_order_perturbation/k_ord_dynare.cpp @@ -36,35 +36,35 @@ KordpDynare::KordpDynare(const char **endo, int num_endo, const char **exo, int nexog, int npar, //const char** par, Vector *ysteady, TwoDMatrix *vcov, Vector *inParams, int nstat, - int npred, int nforw, int nboth, const int jcols, const Vector *nnzd, + int npred, int nforw, int nboth, const int jcols, const Vector *nnzd, const int nsteps, int norder, //const char* modName, - Journal &jr, DynamicModelDLL &dynamicDLL, double sstol, - const vector *var_order, const TwoDMatrix *llincidence, double criterium) throw (TLException) - : nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar), - nYs(npred + nboth), nYss(nboth + nforw), nY(num_endo), nJcols(jcols), NNZD(nnzd), nSteps(nsteps), + Journal &jr, DynamicModelDLL &dynamicDLL, double sstol, + const vector *var_order, const TwoDMatrix *llincidence, double criterium) throw (TLException) : + nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar), + nYs(npred + nboth), nYss(nboth + nforw), nY(num_endo), nJcols(jcols), NNZD(nnzd), nSteps(nsteps), nOrder(norder), journal(jr), dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params(inParams), md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder(var_order), - ll_Incidence(llincidence), qz_criterium(criterium) + ll_Incidence(llincidence), qz_criterium(criterium) { - dnl = new DynareNameList(*this, endo); - denl = new DynareExogNameList(*this, exo); - dsnl = new DynareStateNameList(*this, *dnl, *denl); + dnl = new DynareNameList(*this, endo); + denl = new DynareExogNameList(*this, exo); + dsnl = new DynareStateNameList(*this, *dnl, *denl); - JacobianIndices = ReorderDynareJacobianIndices(varOrder); + JacobianIndices = ReorderDynareJacobianIndices(varOrder); - // Initialise ModelDerivativeContainer(*this, this->md, nOrder); - for (int iord = 1; iord <= nOrder; iord++) - { - FSSparseTensor *t = new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY); - md.insert(t); - } + // Initialise ModelDerivativeContainer(*this, this->md, nOrder); + for (int iord = 1; iord <= nOrder; iord++) + { + FSSparseTensor *t = new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY); + md.insert(t); + } } -KordpDynare::KordpDynare(const KordpDynare &dynare) - : nStat(dynare.nStat), nBoth(dynare.nBoth), nPred(dynare.nPred), +KordpDynare::KordpDynare(const KordpDynare &dynare) : + nStat(dynare.nStat), nBoth(dynare.nBoth), nPred(dynare.nPred), nForw(dynare.nForw), nExog(dynare.nExog), nPar(dynare.nPar), nYs(dynare.nYs), nYss(dynare.nYss), nY(dynare.nY), nJcols(dynare.nJcols), - NNZD(dynare.NNZD), nSteps(dynare.nSteps), nOrder(dynare.nOrder), + NNZD(dynare.NNZD), nSteps(dynare.nSteps), nOrder(dynare.nOrder), journal(dynare.journal), dynamicDLL(dynare.dynamicDLL), ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol), @@ -100,7 +100,7 @@ KordpDynare::~KordpDynare() if (ll_Incidence) delete ll_Incidence; if (NNZD) - delete NNZD; + delete NNZD; } /** This clears the container of model derivatives and initializes it @@ -143,8 +143,8 @@ KordpDynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, void KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareException) { - TwoDMatrix *g2 = 0;// NULL; - TwoDMatrix *g3 = 0;// NULL; + TwoDMatrix *g2 = 0; // NULL; + TwoDMatrix *g3 = 0; // NULL; TwoDMatrix *g1 = new TwoDMatrix(nY, nJcols); // generate g1 for jacobian g1->zeros(); @@ -154,21 +154,21 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareEx if (nOrder > 1) { // generate g2 space for sparse Hessian 3x NNZH = 3x NNZD[1] - g2 = new TwoDMatrix((int) (*NNZD)[1],3); + g2 = new TwoDMatrix((int)(*NNZD)[1], 3); g2->zeros(); } if (nOrder > 2) { - g3 = new TwoDMatrix((int) (*NNZD)[2],3); + g3 = new TwoDMatrix((int)(*NNZD)[2], 3); g3->zeros(); } Vector out(nY); out.zeros(); const Vector *llxYYp; // getting around the constantness if ((nJcols - nExog) > yy.length()) - llxYYp = (LLxSteady(yy)); + llxYYp = (LLxSteady(yy)); else - llxYYp = &yy; + llxYYp = &yy; const Vector &llxYY = *(llxYYp); dynamicDLL.eval(llxYY, xx, params, out, g1, g2, g3); @@ -180,13 +180,13 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareEx delete g1; if (nOrder > 1) { - populateDerivativesContainer(g2, 2, JacobianIndices); - delete g2; + populateDerivativesContainer(g2, 2, JacobianIndices); + delete g2; } if (nOrder > 2) { - populateDerivativesContainer(g3, 3, JacobianIndices); - delete g3; + populateDerivativesContainer(g3, 3, JacobianIndices); + delete g3; } delete llxYYp; } @@ -200,8 +200,8 @@ KordpDynare::calcDerivativesAtSteady() throw (DynareException) } /******************************************************************************* -* populateDerivatives to sparse Tensor and fit it in the Derivatives Container -*******************************************************************************/ + * populateDerivatives to sparse Tensor and fit it in the Derivatives Container + *******************************************************************************/ void KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector *vOrder) { @@ -213,80 +213,80 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vectorncols(); i++) - { - for (int j = 0; j < g->nrows(); j++) - { - double x; - if (s[0] < nJcols-nExog) - x = g->get(j, (*vOrder)[s[0]]); - else - x = g->get(j, s[0]); - if (x != 0.0) - mdTi->insert(s, j, x); - } - s[0]++; - } + { + for (int j = 0; j < g->nrows(); j++) + { + double x; + if (s[0] < nJcols-nExog) + x = g->get(j, (*vOrder)[s[0]]); + else + x = g->get(j, s[0]); + if (x != 0.0) + mdTi->insert(s, j, x); + } + s[0]++; + } } - else if (ord == 2) + else if (ord == 2) { - int nJcols1 = nJcols-nExog; - vector revOrder(nJcols1); - for (int i = 0; i < nJcols1; i++) - revOrder[(*vOrder)[i]] = i; - for (int i = 0; i < g->nrows(); i++) - { - int j = (int)g->get(i,0)-1; // hessian indices start with 1 - int i1 = (int)g->get(i,1) -1; - int s0 = (int)floor(((double) i1)/((double) nJcols)); - int s1 = i1- (nJcols*s0); - if (s0 < nJcols1) - s[0] = revOrder[s0]; - else - s[0] = s0; - if (s1 < nJcols1) - s[1] = revOrder[s1]; - else - s[1] = s1; - if (s[1] >= s[0]) - { - double x = g->get(i,2); - mdTi->insert(s, j, x); - } - } + int nJcols1 = nJcols-nExog; + vector revOrder(nJcols1); + for (int i = 0; i < nJcols1; i++) + revOrder[(*vOrder)[i]] = i; + for (int i = 0; i < g->nrows(); i++) + { + int j = (int) g->get(i, 0)-1; // hessian indices start with 1 + int i1 = (int) g->get(i, 1) -1; + int s0 = (int) floor(((double) i1)/((double) nJcols)); + int s1 = i1- (nJcols*s0); + if (s0 < nJcols1) + s[0] = revOrder[s0]; + else + s[0] = s0; + if (s1 < nJcols1) + s[1] = revOrder[s1]; + else + s[1] = s1; + if (s[1] >= s[0]) + { + double x = g->get(i, 2); + mdTi->insert(s, j, x); + } + } } - else if (ord == 3) + else if (ord == 3) { int nJcols1 = nJcols-nExog; int nJcols2 = nJcols*nJcols; vector revOrder(nJcols1); for (int i = 0; i < nJcols1; i++) - revOrder[(*vOrder)[i]] = i; + revOrder[(*vOrder)[i]] = i; for (int i = 0; i < g->nrows(); i++) - { - int j = (int)g->get(i,0)-1; - int i1 = (int)g->get(i,1) -1; - int s0 = (int)floor(((double) i1)/((double) nJcols2)); - int i2 = i1 - nJcols2*s0; - int s1 = (int)floor(((double) i2)/((double) nJcols)); - int s2 = i2 - nJcols*s1; - if (s0 < nJcols1) - s[0] = revOrder[s0]; - else - s[0] = s0; - if (s1 < nJcols1) - s[1] = revOrder[s1]; - else - s[1] = s1; - if (s2 < nJcols1) - s[2] = revOrder[s2]; - else - s[2] = s2; - if ((s[2] >= s[1]) && (s[1] >= s[0])) - { - double x = g->get(i,2); - mdTi->insert(s, j, x); - } - } + { + int j = (int) g->get(i, 0)-1; + int i1 = (int) g->get(i, 1) -1; + int s0 = (int) floor(((double) i1)/((double) nJcols2)); + int i2 = i1 - nJcols2*s0; + int s1 = (int) floor(((double) i2)/((double) nJcols)); + int s2 = i2 - nJcols*s1; + if (s0 < nJcols1) + s[0] = revOrder[s0]; + else + s[0] = s0; + if (s1 < nJcols1) + s[1] = revOrder[s1]; + else + s[1] = s1; + if (s2 < nJcols1) + s[2] = revOrder[s2]; + else + s[2] = s2; + if ((s[2] >= s[1]) && (s[1] >= s[0])) + { + double x = g->get(i, 2); + mdTi->insert(s, j, x); + } + } } // md container @@ -298,20 +298,20 @@ void KordpDynare::writeModelInfo(Journal &jr) const { // write info on variables - JournalRecordPair rp(journal); - rp << "Information on variables" << endrec; - JournalRecord rec1(journal); - rec1 << "Number of endogenous: " << ny() << endrec; - JournalRecord rec2(journal); - rec2 << "Number of exogenous: " << nexog() << endrec; - JournalRecord rec3(journal); - rec3 << "Number of static: " << nstat() << endrec; - JournalRecord rec4(journal); - rec4 << "Number of predetermined: " << npred()+nboth() << endrec; - JournalRecord rec5(journal); - rec5 << "Number of forward looking: " << nforw()+nboth() << endrec; - JournalRecord rec6(journal); - rec6 << "Number of both: " << nboth() << endrec; + JournalRecordPair rp(journal); + rp << "Information on variables" << endrec; + JournalRecord rec1(journal); + rec1 << "Number of endogenous: " << ny() << endrec; + JournalRecord rec2(journal); + rec2 << "Number of exogenous: " << nexog() << endrec; + JournalRecord rec3(journal); + rec3 << "Number of static: " << nstat() << endrec; + JournalRecord rec4(journal); + rec4 << "Number of predetermined: " << npred()+nboth() << endrec; + JournalRecord rec5(journal); + rec5 << "Number of forward looking: " << nforw()+nboth() << endrec; + JournalRecord rec6(journal); + rec6 << "Number of both: " << nboth() << endrec; } /********************************************************* @@ -328,35 +328,35 @@ KordpDynare::LLxSteady(const Vector &yS) throw (DynareException, TLException) // create temporary square 2D matrix size nEndo x nEndo (sparse) // for the lag, current and lead blocks of the jacobian Vector *llxSteady = new Vector(nJcols-nExog); - for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++) + for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++) + { + // populate (non-sparse) vector with ysteady values + for (int i = 0; i < nY; i++) { - // populate (non-sparse) vector with ysteady values - for (int i = 0; i < nY; i++) - { - if (ll_Incidence->get(ll_row, i)) - (*llxSteady)[((int) ll_Incidence->get(ll_row, i))-1] = yS[i]; - } + if (ll_Incidence->get(ll_row, i)) + (*llxSteady)[((int) ll_Incidence->get(ll_row, i))-1] = yS[i]; } + } return llxSteady; } /************************************ -* Reorder DynareJacobianIndices of variables in a vector according to -* given int * varOrder together with lead & lag incidence matrix and -* any the extra columns for exogenous vars, and then, -* reorders its blocks given by the varOrder and the Dynare++ expectations: + * Reorder DynareJacobianIndices of variables in a vector according to + * given int * varOrder together with lead & lag incidence matrix and + * any the extra columns for exogenous vars, and then, + * reorders its blocks given by the varOrder and the Dynare++ expectations: -* extra nboth+ npred (t-1) lags -* varOrder + * extra nboth+ npred (t-1) lags + * varOrder static: pred both forward -* extra both + nforw (t+1) leads, and -* extra exogen + * extra both + nforw (t+1) leads, and + * extra exogen -* so to match the jacobian organisation expected by the Appoximation class + * so to match the jacobian organisation expected by the Appoximation class both + nforw (t+1) leads static pred @@ -375,37 +375,37 @@ KordpDynare::ReorderDynareJacobianIndices(const vector *varOrder) throw (TL vector tmp(nY); int i, j, rjoff = nJcols-nExog-1; - for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++) - { - // reorder in orde-var order & populate temporary nEndo (sparse) vector with - // the lag, current and lead blocks of the jacobian respectively - for (i = 0; i < nY; i++) - tmp[i] = ((int) ll_Incidence->get(ll_row, (*varOrder)[i]-1)); - // write the reordered blocks back to the jacobian - // in reverse order - for (j = nY-1; j >= 0; j--) - if (tmp[j]) - { - (*JacobianIndices)[rjoff] = tmp[j] -1; - rjoff--; - if (rjoff < 0) - break; - } - } + for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++) + { + // reorder in orde-var order & populate temporary nEndo (sparse) vector with + // the lag, current and lead blocks of the jacobian respectively + for (i = 0; i < nY; i++) + tmp[i] = ((int) ll_Incidence->get(ll_row, (*varOrder)[i]-1)); + // write the reordered blocks back to the jacobian + // in reverse order + for (j = nY-1; j >= 0; j--) + if (tmp[j]) + { + (*JacobianIndices)[rjoff] = tmp[j] -1; + rjoff--; + if (rjoff < 0) + break; + } + } //add the indices for the nExog exogenous jacobians for (j = nJcols-nExog; j < nJcols; j++) - (*JacobianIndices)[j] = j; + (*JacobianIndices)[j] = j; return JacobianIndices; } /************************************ -* Reorder first set of columns of variables in a (jacobian) matrix -* according to order given in varsOrder together with the extras -* assuming tdx ncols() - nExog is eaqual or less than length of varOrder and -* of any of its elements too. -************************************/ + * Reorder first set of columns of variables in a (jacobian) matrix + * according to order given in varsOrder together with the extras + * assuming tdx ncols() - nExog is eaqual or less than length of varOrder and + * of any of its elements too. + ************************************/ void KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector *vOrder) throw (DynareException, TLException) @@ -419,8 +419,8 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector *vOrder) throw (Dyna tdx->zeros(); // empty original matrix // reorder the columns - for (int i = 0; i < tdx->ncols(); i++) - tdx->copyColumn(tmpR, (*vOrder)[i], i); + for (int i = 0; i < tdx->ncols(); i++) + tdx->copyColumn(tmpR, (*vOrder)[i], i); } void @@ -431,15 +431,15 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const int *vOrder) throw (TLException) TwoDMatrix &tmpR = tmp; tdx->zeros(); // empty original matrix // reorder the columns - for (int i = 0; i < tdx->ncols(); i++) - tdx->copyColumn(tmpR, vOrder[i], i); + for (int i = 0; i < tdx->ncols(); i++) + tdx->copyColumn(tmpR, vOrder[i], i); } /*********************************************************************** -* Recursive hierarchical block reordering of the higher order, input model -* derivatives inc. Hessian -* This is now obsolete but kept in in case it is needed -***********************************************************************/ + * Recursive hierarchical block reordering of the higher order, input model + * derivatives inc. Hessian + * This is now obsolete but kept in in case it is needed + ***********************************************************************/ void KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector *vOrder) throw (DynareException, TLException) @@ -476,8 +476,8 @@ KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector *vOrder) throw (Dy throw DynareException(__FILE__, __LINE__, "Size of order var is too small"); // reorder the columns - for (int i = 0; i < tdx->ncols(); i++) - tdx->copyColumn(tmpR, (*vOrder)[i], i); + for (int i = 0; i < tdx->ncols(); i++) + tdx->copyColumn(tmpR, (*vOrder)[i], i); } } @@ -514,7 +514,7 @@ DynareNameList::selectIndices(const vector &ns) const throw (Dynar DynareNameList::DynareNameList(const KordpDynare &dynare) { for (int i = 0; i < dynare.ny(); i++) - names.push_back(dynare.dnl->getName(i)); + names.push_back(dynare.dnl->getName(i)); } DynareNameList::DynareNameList(const KordpDynare &dynare, const char **namesp) { @@ -538,7 +538,7 @@ DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const Dynare const DynareExogNameList &denl) { for (int i = 0; i < dynare.nys(); i++) - names.push_back(dnl.getName(i+dynare.nstat())); + names.push_back(dnl.getName(i+dynare.nstat())); for (int i = 0; i < dynare.nexog(); i++) - names.push_back(denl.getName(i)); + names.push_back(denl.getName(i)); } diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.h b/mex/sources/k_order_perturbation/k_ord_dynare.h index 8b6180bce..7a5323cf7 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.h +++ b/mex/sources/k_order_perturbation/k_ord_dynare.h @@ -144,10 +144,10 @@ public: KordpDynare(const char **endo, int num_endo, const char **exo, int num_exo, int num_par, Vector *ySteady, TwoDMatrix *vCov, Vector *params, int nstat, int nPred, - int nforw, int nboth, const int nJcols, const Vector *NNZD, + int nforw, int nboth, const int nJcols, const Vector *NNZD, const int nSteps, const int ord, Journal &jr, DynamicModelDLL &dynamicDLL, double sstol, - const vector *varOrder, const TwoDMatrix *ll_Incidence, + const vector *varOrder, const TwoDMatrix *ll_Incidence, double qz_criterium) throw (TLException); /** Makes a deep copy of the object. */ @@ -291,8 +291,8 @@ class KordpVectorFunction : public ogu::VectorFunction protected: KordpDynare &d; public: - KordpVectorFunction(KordpDynare &dyn) - : d(dyn) + KordpVectorFunction(KordpDynare &dyn) : + d(dyn) { } virtual ~KordpVectorFunction() diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cpp b/mex/sources/k_order_perturbation/k_order_perturbation.cpp index 964187071..226a91995 100644 --- a/mex/sources/k_order_perturbation/k_order_perturbation.cpp +++ b/mex/sources/k_order_perturbation/k_order_perturbation.cpp @@ -18,24 +18,24 @@ */ /****************************************************** - // k_order_perturbation.cpp : Defines the entry point for the k-order perturbation application DLL. - // - // called from Dynare dr1_k_order.m, (itself called form resol.m instead of regular dr1.m) - // if options_.order < 2 % 1st order only - // [ysteady, ghx_u]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]); - // else % 2nd order - // [ysteady, ghx_u, g_2]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]); - // inputs: - // dr, - Dynare structure - // task, - check or not, not used - // M_ - Dynare structure - // options_ - Dynare structure - // oo_ - Dynare structure - // ['.' mexext] Matlab dll extension - // returns: - // ysteady steady state - // ghx_u - first order rules packed in one matrix - // g_2 - 2nd order rules packed in one matrix + // k_order_perturbation.cpp : Defines the entry point for the k-order perturbation application DLL. + // + // called from Dynare dr1_k_order.m, (itself called form resol.m instead of regular dr1.m) + // if options_.order < 2 % 1st order only + // [ysteady, ghx_u]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]); + // else % 2nd order + // [ysteady, ghx_u, g_2]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]); + // inputs: + // dr, - Dynare structure + // task, - check or not, not used + // M_ - Dynare structure + // options_ - Dynare structure + // oo_ - Dynare structure + // ['.' mexext] Matlab dll extension + // returns: + // ysteady steady state + // ghx_u - first order rules packed in one matrix + // g_2 - 2nd order rules packed in one matrix **********************************************************/ #include "k_ord_dynare.h" @@ -84,7 +84,7 @@ DynareMxArrayToString(const char *cNamesCharStr, const int len, const int width) { char **cNamesMX; cNamesMX = (char **) calloc(len, sizeof(char *)); - for(int i = 0; i < len; i++) + for (int i = 0; i < len; i++) cNamesMX[i] = (char *) calloc(width+1, sizeof(char)); for (int i = 0; i < width; i++) @@ -99,7 +99,7 @@ DynareMxArrayToString(const char *cNamesCharStr, const int len, const int width) else cNamesMX[j][i] = '\0'; } } - return (const char **)cNamesMX; + return (const char **) cNamesMX; } ////////////////////////////////////////////////////// @@ -112,7 +112,7 @@ DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width) { char *cNamesCharStr = mxArrayToString(mxFldp); const char **ret = DynareMxArrayToString(cNamesCharStr, len, width); - + return ret; } @@ -127,7 +127,7 @@ extern "C" { mexErrMsgTxt("Must have at least 5 input parameters."); if (nlhs == 0) mexErrMsgTxt("Must have at least 1 output parameter."); - + const mxArray *dr = prhs[0]; const int check_flag = (int) mxGetScalar(prhs[1]); const mxArray *M_ = prhs[2]; @@ -136,7 +136,7 @@ extern "C" { mxArray *mFname = mxGetField(M_, 0, "fname"); if (!mxIsChar(mFname)) - mexErrMsgTxt("Input must be of type char."); + mexErrMsgTxt("Input must be of type char."); string fName = mxArrayToString(mFname); const mxArray *mexExt = prhs[5]; string dfExt = mxArrayToString(mexExt); //Dynamic file extension, e.g.".dll" or .mexw32; @@ -147,12 +147,12 @@ extern "C" { kOrder = (int) mxGetScalar(mxFldp); else kOrder = 1; - + if (kOrder == 1 && nlhs != 1) mexErrMsgTxt("k_order_perturbation at order 1 requires exactly 1 argument in output"); else if (kOrder > 1 && nlhs != kOrder+1) mexErrMsgTxt("k_order_perturbation at order > 1 requires exactly order + 1 argument in output"); - + double qz_criterium = 1+1e-6; mxFldp = mxGetField(options_, 0, "qz_criterium"); if (mxIsNumeric(mxFldp)) @@ -201,11 +201,11 @@ extern "C" { mxFldp = mxGetField(dr, 0, "order_var"); dparams = (double *) mxGetData(mxFldp); npar = (int) mxGetM(mxFldp); - if (npar != nEndo) //(nPar != npar) - mexErrMsgTxt("Incorrect number of input var_order vars."); + if (npar != nEndo) //(nPar != npar) + mexErrMsgTxt("Incorrect number of input var_order vars."); vector *var_order_vp = (new vector(nEndo)); for (int v = 0; v < nEndo; v++) - (*var_order_vp)[v] = (int)(*(dparams++)); + (*var_order_vp)[v] = (int)(*(dparams++)); // the lag, current and lead blocks of the jacobian respectively mxFldp = mxGetField(M_, 0, "lead_lag_incidence"); @@ -213,15 +213,14 @@ extern "C" { npar = (int) mxGetN(mxFldp); int nrows = (int) mxGetM(mxFldp); - TwoDMatrix *llincidence = new TwoDMatrix(nrows, npar, dparams); if (npar != nEndo) mexErrMsgIdAndTxt("dynare:k_order_perturbation", "Incorrect length of lead lag incidences: ncol=%d != nEndo=%d.", npar, nEndo); - //get NNZH =NNZD(2) = the total number of non-zero Hessian elements + //get NNZH =NNZD(2) = the total number of non-zero Hessian elements mxFldp = mxGetField(M_, 0, "NNZDerivatives"); dparams = (double *) mxGetData(mxFldp); - Vector *NNZD = new Vector (dparams, (int) mxGetM(mxFldp)); + Vector *NNZD = new Vector(dparams, (int) mxGetM(mxFldp)); const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns @@ -235,7 +234,7 @@ extern "C" { const int widthExog = (int) mxGetN(mxFldp); const char **exoNamesMX = DynareMxArrayToString(mxFldp, nexo, widthExog); - if ((nEndo != nendo) || (nExog != nexo)) //(nPar != npar) + if ((nEndo != nendo) || (nExog != nexo)) //(nPar != npar) mexErrMsgTxt("Incorrect number of input parameters."); /* Fetch time index */ @@ -259,11 +258,11 @@ extern "C" { // make KordpDynare object KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, - jcols, NNZD, nSteps, kOrder, journal, dynamicDLL, + jcols, NNZD, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp, llincidence, qz_criterium); // construct main K-order approximation class - + Approximation app(dynare, journal, nSteps, false, qz_criterium); // run stochastic steady app.walkStochSteady(); @@ -301,11 +300,11 @@ extern "C" { if (kOrder == 1) { /* Set the output pointer to the output matrix ysteady. */ - map::const_iterator cit = mm.begin(); - ++cit; + map::const_iterator cit = mm.begin(); + ++cit; plhs[0] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL); - TwoDMatrix g((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[0])); - g = (const TwoDMatrix &)(*cit).second; + TwoDMatrix g((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[0])); + g = (const TwoDMatrix &)(*cit).second; } if (kOrder >= 2) { diff --git a/mex/sources/k_order_perturbation/tests/first_order.m b/mex/sources/k_order_perturbation/tests/first_order.m index 995ec3af7..f46631c22 100644 --- a/mex/sources/k_order_perturbation/tests/first_order.m +++ b/mex/sources/k_order_perturbation/tests/first_order.m @@ -1,129 +1,129 @@ -function [gy]=first_order(M_, dr, jacobia) -% Emulation of Dynare++ c++ first_order.cpp for testing pruposes - -% Copyright (C) 2009 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -% fd = jacobia_ -% reorder jacobia_ -[fd]=k_reOrderedJacobia(M_, jacobia) - -%ypart=dr; -ypart.ny=M_.endo_nbr; -ypart.nyss=dr.nfwrd+dr.nboth; -ypart.nys=dr.npred; -ypart.npred=dr.npred-dr.nboth; -ypart.nboth=dr.nboth; -ypart.nforw=dr.nfwrd; -ypart.nstat =dr.nstatic -nu=M_.exo_nbr - - off= 1; % = 0 in C - fyplus = fd(:,off:off+ypart.nyss-1); - off= off+ypart.nyss; - fyszero= fd(:,off:off+ypart.nstat-1); - off= off+ypart.nstat; - fypzero= fd(:,off:off+ypart.npred-1); - off= off+ypart.npred; - fybzero= fd(:,off:off+ypart.nboth-1); - off= off+ypart.nboth; - fyfzero= fd(:,off:off+ypart.nforw-1); - off= off+ypart.nforw; - fymins= fd(:,off:off+ypart.nys-1); - off= off+ypart.nys; - fuzero= fd(:,off:off+nu-1); - off=off+ nu; - - n= ypart.ny+ypart.nboth; - %TwoDMatrix - matD=zeros(n,n); -% matD.place(fypzero,0,0); - matD(1:n-ypart.nboth,1:ypart.npred)= fypzero; -% matD.place(fybzero,0,ypart.npred); - matD(1:n-ypart.nboth,ypart.npred+1:ypart.npred+ypart.nboth)=fybzero; -% matD.place(fyplus,0,ypart.nys()+ypart.nstat); - matD(1:n-ypart.nboth,ypart.nys+ypart.nstat+1:ypart.nys+ypart.nstat+ypart.nyss)=fyplus; - for i=1:ypart.nboth - matD(ypart.ny()+i,ypart.npred+i)= 1.0; - end - - matE=[fymins, fyszero, zeros(n-ypart.nboth,ypart.nboth), fyfzero; zeros(ypart.nboth,n)]; -% matE.place(fymins; -% matE.place(fyszero,0,ypart.nys()); -% matE.place(fyfzero,0,ypart.nys()+ypart.nstat+ypart.nboth); - - for i= 1:ypart.nboth - matE(ypart.ny()+i,ypart.nys()+ypart.nstat+i)= -1.0; - end - matE=-matE; %matE.mult(-1.0); - -% vsl=zeros(n,n); -% vsr=zeros(n,n); -% lwork= 100*n+16; -% work=zeros(1,lwork); -% bwork=zeros(1,n); - %int info; - -% LAPACK_dgges("N","V","S",order_eigs,&n,matE.getData().base(),&n, -% matD.getData().base(),&n,&sdim,alphar.base(),alphai.base(), -% beta.base(),vsl.getData().base(),&n,vsr.getData().base(),&n, -% work.base(),&lwork,&(bwork[0]),&info); - - [matE1,matD1,vsr,sdim,dr.eigval,info] = mjdgges(matE,matD,1); - - bk_cond= (sdim==ypart.nys); - -% ConstGeneralMatrix z11(vsr,0,0,ypart.nys(),ypart.nys()); - z11=vsr(1:ypart.nys,1:ypart.nys); -% ConstGeneralMatrix z12(vsr,0,ypart.nys(),ypart.nys(),n-ypart.nys()); - z12=vsr(1:ypart.nys(),ypart.nys+1:end);%, n-ypart.nys); -% ConstGeneralMatrix z21(vsr,ypart.nys(),0,n-ypart.nys(),ypart.nys()); - z21=vsr(ypart.nys+1:end,1:ypart.nys); -% ConstGeneralMatrix z22(vsr,ypart.nys(),ypart.nys(),n-ypart.nys(),n-ypart.nys()); - z22=vsr(ypart.nys+1:end,ypart.nys+1:end); - -% GeneralMatrix sfder(z12,"transpose"); - sfder=z12';%,"transpose"); -% z22.multInvLeftTrans(sfder); - sfder=z22'\sfder; - sfder=-sfder;% .mult(-1); - - %s11(matE,0,0,ypart.nys(),ypart.nys()); - s11=matE1(1:ypart.nys,1:ypart.nys); -% t11=(matD1,0,0,ypart.nys(),ypart.nys()); - t11=matD1(1:ypart.nys,1:ypart.nys); - dumm=(s11');%,"transpose"); - %z11.multInvLeftTrans(dumm); - dumm=z11'\dumm; - preder=(dumm');%,"transpose"); - %t11.multInvLeft(preder); - preder=t11\preder; - %preder.multLeft(z11); - preder= z11*preder; - -% gy.place(preder,ypart.nstat,0); -% gy=(zeros(ypart.nstat,size(preder,2)) ;preder); -% sder(sfder,0,0,ypart.nstat,ypart.nys()); - sder=sfder(1:ypart.nstat,1:ypart.nys); -% gy.place(sder,0,0); -% gy(1:ypart.nstat, 1:ypart.nys)=sder; -% gy=[sder;preder]; -% fder(sfder,ypart.nstat+ypart.nboth,0,ypart.nforw,ypart.nys()); - fder=sfder(ypart.nstat+ypart.nboth+1:ypart.nstat+ypart.nboth+ypart.nforw,1:ypart.nys); -% gy.place(fder,ypart.nstat+ypart.nys(),0); -% gy(ypart.nstat+ypart.nys,1)=fder; - gy=[sder;preder;fder]; +function [gy]=first_order(M_, dr, jacobia) +% Emulation of Dynare++ c++ first_order.cpp for testing pruposes + +% Copyright (C) 2009 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% fd = jacobia_ +% reorder jacobia_ +[fd]=k_reOrderedJacobia(M_, jacobia) + +%ypart=dr; +ypart.ny=M_.endo_nbr; +ypart.nyss=dr.nfwrd+dr.nboth; +ypart.nys=dr.npred; +ypart.npred=dr.npred-dr.nboth; +ypart.nboth=dr.nboth; +ypart.nforw=dr.nfwrd; +ypart.nstat =dr.nstatic +nu=M_.exo_nbr + +off= 1; % = 0 in C +fyplus = fd(:,off:off+ypart.nyss-1); +off= off+ypart.nyss; +fyszero= fd(:,off:off+ypart.nstat-1); +off= off+ypart.nstat; +fypzero= fd(:,off:off+ypart.npred-1); +off= off+ypart.npred; +fybzero= fd(:,off:off+ypart.nboth-1); +off= off+ypart.nboth; +fyfzero= fd(:,off:off+ypart.nforw-1); +off= off+ypart.nforw; +fymins= fd(:,off:off+ypart.nys-1); +off= off+ypart.nys; +fuzero= fd(:,off:off+nu-1); +off=off+ nu; + +n= ypart.ny+ypart.nboth; +%TwoDMatrix +matD=zeros(n,n); +% matD.place(fypzero,0,0); +matD(1:n-ypart.nboth,1:ypart.npred)= fypzero; +% matD.place(fybzero,0,ypart.npred); +matD(1:n-ypart.nboth,ypart.npred+1:ypart.npred+ypart.nboth)=fybzero; +% matD.place(fyplus,0,ypart.nys()+ypart.nstat); +matD(1:n-ypart.nboth,ypart.nys+ypart.nstat+1:ypart.nys+ypart.nstat+ypart.nyss)=fyplus; +for i=1:ypart.nboth + matD(ypart.ny()+i,ypart.npred+i)= 1.0; +end + +matE=[fymins, fyszero, zeros(n-ypart.nboth,ypart.nboth), fyfzero; zeros(ypart.nboth,n)]; +% matE.place(fymins; +% matE.place(fyszero,0,ypart.nys()); +% matE.place(fyfzero,0,ypart.nys()+ypart.nstat+ypart.nboth); + +for i= 1:ypart.nboth + matE(ypart.ny()+i,ypart.nys()+ypart.nstat+i)= -1.0; +end +matE=-matE; %matE.mult(-1.0); + +% vsl=zeros(n,n); +% vsr=zeros(n,n); +% lwork= 100*n+16; +% work=zeros(1,lwork); +% bwork=zeros(1,n); +%int info; + +% LAPACK_dgges("N","V","S",order_eigs,&n,matE.getData().base(),&n, +% matD.getData().base(),&n,&sdim,alphar.base(),alphai.base(), +% beta.base(),vsl.getData().base(),&n,vsr.getData().base(),&n, +% work.base(),&lwork,&(bwork[0]),&info); + +[matE1,matD1,vsr,sdim,dr.eigval,info] = mjdgges(matE,matD,1); + +bk_cond= (sdim==ypart.nys); + +% ConstGeneralMatrix z11(vsr,0,0,ypart.nys(),ypart.nys()); +z11=vsr(1:ypart.nys,1:ypart.nys); +% ConstGeneralMatrix z12(vsr,0,ypart.nys(),ypart.nys(),n-ypart.nys()); +z12=vsr(1:ypart.nys(),ypart.nys+1:end);%, n-ypart.nys); + % ConstGeneralMatrix z21(vsr,ypart.nys(),0,n-ypart.nys(),ypart.nys()); +z21=vsr(ypart.nys+1:end,1:ypart.nys); +% ConstGeneralMatrix z22(vsr,ypart.nys(),ypart.nys(),n-ypart.nys(),n-ypart.nys()); +z22=vsr(ypart.nys+1:end,ypart.nys+1:end); + +% GeneralMatrix sfder(z12,"transpose"); +sfder=z12';%,"transpose"); + % z22.multInvLeftTrans(sfder); +sfder=z22'\sfder; +sfder=-sfder;% .mult(-1); + +%s11(matE,0,0,ypart.nys(),ypart.nys()); +s11=matE1(1:ypart.nys,1:ypart.nys); +% t11=(matD1,0,0,ypart.nys(),ypart.nys()); +t11=matD1(1:ypart.nys,1:ypart.nys); +dumm=(s11');%,"transpose"); + %z11.multInvLeftTrans(dumm); +dumm=z11'\dumm; +preder=(dumm');%,"transpose"); + %t11.multInvLeft(preder); +preder=t11\preder; +%preder.multLeft(z11); +preder= z11*preder; + +% gy.place(preder,ypart.nstat,0); +% gy=(zeros(ypart.nstat,size(preder,2)) ;preder); +% sder(sfder,0,0,ypart.nstat,ypart.nys()); +sder=sfder(1:ypart.nstat,1:ypart.nys); +% gy.place(sder,0,0); +% gy(1:ypart.nstat, 1:ypart.nys)=sder; +% gy=[sder;preder]; +% fder(sfder,ypart.nstat+ypart.nboth,0,ypart.nforw,ypart.nys()); +fder=sfder(ypart.nstat+ypart.nboth+1:ypart.nstat+ypart.nboth+ypart.nforw,1:ypart.nys); +% gy.place(fder,ypart.nstat+ypart.nys(),0); +% gy(ypart.nstat+ypart.nys,1)=fder; +gy=[sder;preder;fder]; diff --git a/mex/sources/k_order_perturbation/tests/k_order_test_main.cpp b/mex/sources/k_order_perturbation/tests/k_order_test_main.cpp index 978a6981d..f1241c447 100644 --- a/mex/sources/k_order_perturbation/tests/k_order_test_main.cpp +++ b/mex/sources/k_order_perturbation/tests/k_order_test_main.cpp @@ -81,7 +81,7 @@ main(int argc, char *argv[]) const int nSteady = 16; //27 //31;//29, 16 (int)mxGetM(mxFldp); Vector *ySteady = new Vector(dYSparams, nSteady); - double nnzd[3]={ 77,217,0}; + double nnzd[3] = { 77, 217, 0}; const Vector *NNZD = new Vector(nnzd, 3); //mxFldp = mxGetField(dr, 0,"nstatic" ); @@ -111,7 +111,7 @@ main(int argc, char *argv[]) = { 5, 6, 8, 10, 11, 12, 14, 7, 13, 1, 2, 3, 4, 9, 15, 16 // 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17, 18 - }; + }; //Vector * varOrder = new Vector(var_order, nEndo); vector *var_order_vp = new vector(nEndo); //nEndo)); for (int v = 0; v < nEndo; v++) @@ -135,7 +135,7 @@ main(int argc, char *argv[]) 0, 18, 0, 0, 19, 26, 0, 20, 27 - }; + }; TwoDMatrix *llincidence = new TwoDMatrix(3, nEndo, ll_incidence); const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns diff --git a/mex/sources/kronecker/A_times_B_kronecker_C.cc b/mex/sources/kronecker/A_times_B_kronecker_C.cc index 13c0abd5f..6feb68f42 100644 --- a/mex/sources/kronecker/A_times_B_kronecker_C.cc +++ b/mex/sources/kronecker/A_times_B_kronecker_C.cc @@ -18,8 +18,8 @@ */ /* - * This mex file computes A*kron(B,C) or A*kron(B,B) without explicitely building kron(B,C) or kron(B,B), so that - * one can consider large matrices B and/or C. + * This mex file computes A*kron(B,C) or A*kron(B,B) without explicitely building kron(B,C) or kron(B,B), so that + * one can consider large matrices B and/or C. */ #include @@ -28,105 +28,107 @@ #include #ifdef USE_OMP - #include +# include #endif -void full_A_times_kronecker_B_C(double *A, double *B, double *C, double *D, - blas_int mA, blas_int nA, blas_int mB, blas_int nB, blas_int mC, blas_int nC) +void +full_A_times_kronecker_B_C(double *A, double *B, double *C, double *D, + blas_int mA, blas_int nA, blas_int mB, blas_int nB, blas_int mC, blas_int nC) { - #if USE_OMP - #pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) - for(long int colD = 0; colD 3) || (nrhs <2) ) + if ((nrhs > 3) || (nrhs < 2)) { mexErrMsgTxt("Two or Three input arguments required."); } - if (nlhs>1) + if (nlhs > 1) { mexErrMsgTxt("Too many output arguments."); } @@ -136,21 +138,21 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) nA = mxGetN(prhs[0]); mB = mxGetM(prhs[1]); nB = mxGetN(prhs[1]); - if (nrhs == 3)// A*kron(B,C) is to be computed. + if (nrhs == 3) // A*kron(B,C) is to be computed. { mC = mxGetM(prhs[2]); nC = mxGetN(prhs[2]); if (mB*mC != nA) - { - mexErrMsgTxt("Input dimension error!"); - } + { + mexErrMsgTxt("Input dimension error!"); + } } - else// A*kron(B,B) is to be computed. + else // A*kron(B,B) is to be computed. { if (mB*mB != nA) - { - mexErrMsgTxt("Input dimension error!"); - } + { + mexErrMsgTxt("Input dimension error!"); + } } // Get input matrices: double *B, *C, *A; @@ -164,11 +166,11 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) double *D; if (nrhs == 3) { - plhs[0] = mxCreateDoubleMatrix(mA,nB*nC,mxREAL); + plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL); } else { - plhs[0] = mxCreateDoubleMatrix(mA,nB*nB,mxREAL); + plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL); } D = mxGetPr(plhs[0]); // Computational part: diff --git a/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc b/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc index 7a72d36b8..0fe664d73 100644 --- a/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc +++ b/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc @@ -19,7 +19,7 @@ /* * This mex file computes A*kron(B,C) or A*kron(B,B) without explicitly building kron(B,C) or kron(B,B), so that - * one can consider large matrices A, B and/or C, and assuming that A is a the hessian of a dsge model + * one can consider large matrices A, B and/or C, and assuming that A is a the hessian of a dsge model * (dynare format). This mex file should not be used outside dr1.m. */ @@ -31,77 +31,79 @@ # include #endif -void sparse_hessian_times_B_kronecker_B(mwIndex *isparseA, mwIndex *jsparseA, double *vsparseA, - double *B, double *D, mwSize mA, mwSize nA, mwSize mB, mwSize nB) +void +sparse_hessian_times_B_kronecker_B(mwIndex *isparseA, mwIndex *jsparseA, double *vsparseA, + double *B, double *D, mwSize mA, mwSize nA, mwSize mB, mwSize nB) { - /* + /* ** Loop over the columns of kron(B,B) (or of the result matrix D). ** This loop is splitted into two nested loops because we use the - ** symmetric pattern of the hessian matrix. + ** symmetric pattern of the hessian matrix. */ - #if USE_OMP - #pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) - #endif - for(int j1B=0; j1B0) - { - memcpy(&D[(j2B*nB+j1B)*mA],&D[jj*mA],mA*sizeof(double)); - } - } + /* + ** Loop over the rows of kron(B,B) (column jj). + */ + for (unsigned long int ii = 0; ii < nA; ii++) + { + k1 = jsparseA[ii]; + k2 = jsparseA[ii+1]; + if (k1 < k2) // otherwise column ii of A does not have non zero elements (and there is nothing to compute). + { + ++nz_in_column_ii_of_A; + unsigned int i1B = (ii/mB); + unsigned int i2B = (ii%mB); + double bb = B[j1B*mB+i1B]*B[j2B*mB+i2B]; + /* + ** Loop over the non zero entries of A(:,ii). + */ + for (unsigned int k = k1; k < k2; k++) + { + unsigned int kk = isparseA[k]; + D[jj*mA+kk] = D[jj*mA+kk] + bb*vsparseA[iv]; + iv++; + } + } + } + if (nz_in_column_ii_of_A > 0) + { + memcpy(&D[(j2B*nB+j1B)*mA], &D[jj*mA], mA*sizeof(double)); + } + } } } -void sparse_hessian_times_B_kronecker_C(mwIndex *isparseA, mwIndex *jsparseA, double *vsparseA, - double *B, double *C, double *D, - mwSize mA, mwSize nA, mwSize mB, mwSize nB, mwSize mC, mwSize nC) +void +sparse_hessian_times_B_kronecker_C(mwIndex *isparseA, mwIndex *jsparseA, double *vsparseA, + double *B, double *C, double *D, + mwSize mA, mwSize nA, mwSize mB, mwSize nB, mwSize mC, mwSize nC) { - /* + /* ** Loop over the columns of kron(B,B) (or of the result matrix D). */ - #if USE_OMP - #pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) - #endif - for(long int jj=0; jj 3) || (nrhs <2) ) + if ((nrhs > 3) || (nrhs < 2)) { mexErrMsgTxt("Two or Three input arguments required."); } - if (nlhs>1) + if (nlhs > 1) { mexErrMsgTxt("Too many output arguments."); } @@ -156,21 +159,21 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) nA = mxGetN(prhs[0]); mB = mxGetM(prhs[1]); nB = mxGetN(prhs[1]); - if (nrhs == 3)// A*kron(B,C) is to be computed. + if (nrhs == 3) // A*kron(B,C) is to be computed. { mC = mxGetM(prhs[2]); nC = mxGetN(prhs[2]); if (mB*mC != nA) - { - mexErrMsgTxt("Input dimension error!"); - } + { + mexErrMsgTxt("Input dimension error!"); + } } - else// A*kron(B,B) is to be computed. + else // A*kron(B,B) is to be computed. { if (mB*mB != nA) - { - mexErrMsgTxt("Input dimension error!"); - } + { + mexErrMsgTxt("Input dimension error!"); + } } // Get input matrices: double *B, *C; @@ -180,23 +183,23 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) C = mxGetPr(prhs[2]); } // Sparse (dynare) hessian matrix. - mwIndex *isparseA = (mwIndex*)mxGetIr(prhs[0]); - mwIndex *jsparseA = (mwIndex*)mxGetJc(prhs[0]); + mwIndex *isparseA = (mwIndex *) mxGetIr(prhs[0]); + mwIndex *jsparseA = (mwIndex *) mxGetJc(prhs[0]); double *vsparseA = mxGetPr(prhs[0]); // Initialization of the ouput: double *D; if (nrhs == 3) { - plhs[0] = mxCreateDoubleMatrix(mA,nB*nC,mxREAL); + plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL); } else { - plhs[0] = mxCreateDoubleMatrix(mA,nB*nB,mxREAL); + plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL); } D = mxGetPr(plhs[0]); // Computational part: if (nrhs == 2) - { + { sparse_hessian_times_B_kronecker_B(isparseA, jsparseA, vsparseA, B, D, mA, nA, mB, nB); } else diff --git a/mex/sources/kronecker/tests/test_kron.m b/mex/sources/kronecker/tests/test_kron.m index 1e95a176e..9be8a775d 100644 --- a/mex/sources/kronecker/tests/test_kron.m +++ b/mex/sources/kronecker/tests/test_kron.m @@ -16,82 +16,82 @@ function test_kron(test) % % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . + +if ~nargin + test = 3; +end + + +if test == 1 - if ~nargin - test = 3; + percentage_of_non_zero_elements = 10e-4; + NumberOfVariables = 549;%100; + NumberOfEquations = 256;%50 + NumberOfColsInB = 50 ; + A = zeros(NumberOfEquations,NumberOfVariables^2); + for i = 1:NumberOfEquations + for j = 1:NumberOfVariables + for k = j:NumberOfVariables + if rand. - */ - -#include - -#include -#include - -double criterium; - -lapack_int -my_criteria(const double *alphar, const double *alphai, const double *beta) -{ - return( (*alphar * *alphar + *alphai * *alphai) < criterium * *beta * *beta); -} - -void mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r, double *eval_i, double *info) -{ - lapack_int i_n, i_info, i_sdim, one, lwork; - double *alphar, *alphai, *beta, *work, *par, *pai, *pb, *per, *pei; - double *junk; - lapack_int *bwork; - - one = 1; - i_n = (lapack_int)*n; - alphar = mxCalloc(i_n,sizeof(double)); - alphai = mxCalloc(i_n,sizeof(double)); - beta = mxCalloc(i_n,sizeof(double)); - lwork = 16*i_n+16; - work = mxCalloc(lwork,sizeof(double)); - bwork = mxCalloc(i_n,sizeof(lapack_int)); - /* made necessary by bug in Lapack */ - junk = mxCalloc(i_n*i_n,sizeof(double)); - - dgges("N", "V", "S", my_criteria, &i_n, a, &i_n, b, &i_n, &i_sdim, alphar, alphai, beta, junk, &i_n, z, &i_n, work, &lwork, bwork, &i_info); - - *sdim = i_sdim; - *info = i_info; - - par = alphar; - pai = alphai; - pb = beta; - pei = eval_i; - for(per = eval_r; per <= &eval_r[i_n-1]; ++per) - { - *per = *par / *pb; - *pei = *pai / *pb; - ++par; - ++pai; - ++pb; - ++pei; - } -} - -/* MATLAB interface */ -void mexFunction( int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[] ) - -{ - unsigned int m1,n1,m2,n2; - double *s, *t, *z, *sdim, *eval_r, *eval_i, *info, *a, *b; - double n; - - /* Check for proper number of arguments */ - - if (nrhs < 2 || nrhs > 3) { - mexErrMsgTxt("MJDGGES: two or three input arguments are required."); - } else if (nlhs > 6) { - mexErrMsgTxt("MJDGGES: too many output arguments."); - } - - /* Check that A and B are real matrices of the same dimension.*/ - - m1 = mxGetM(prhs[0]); - n1 = mxGetN(prhs[0]); - m2 = mxGetM(prhs[1]); - n2 = mxGetN(prhs[1]); - if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || - !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || - (m1 != n1) || (m2!= n1) || (m2 != n2)) { - mexErrMsgTxt("MJDGGES requires two square real matrices of the same dimension."); - } - - /* Create a matrix for the return argument */ - plhs[0] = mxCreateDoubleMatrix(n1, n1, mxREAL); - plhs[1] = mxCreateDoubleMatrix(n1, n1, mxREAL); - plhs[2] = mxCreateDoubleMatrix(n1, n1, mxREAL); - plhs[3] = mxCreateDoubleMatrix(1, 1, mxREAL); - plhs[4] = mxCreateDoubleMatrix(n1, 1, mxCOMPLEX); - plhs[5] = mxCreateDoubleMatrix(1, 1, mxREAL); - - /* Assign pointers to the various parameters */ - s = mxGetPr(plhs[0]); - t = mxGetPr(plhs[1]); - z = mxGetPr(plhs[2]); - sdim = mxGetPr(plhs[3]); - eval_r = mxGetPr(plhs[4]); - eval_i = mxGetPi(plhs[4]); - info = mxGetPr(plhs[5]); - - a = mxGetPr(prhs[0]); - b = mxGetPr(prhs[1]); - - /* set criterium for stable eigenvalues */ - if ( nrhs == 3) - { - criterium = *mxGetPr(prhs[2]); - } - else - { - criterium = 1+1e-6; - } - - /* keep a and b intact */ - memcpy(s,a,sizeof(double)*n1*n1); - memcpy(t,b,sizeof(double)*n1*n1); - - n = n1; - - /* Do the actual computations in a subroutine */ - mjdgges(s, t, z, &n, sdim, eval_r, eval_i, info); - - - return; - -} - -/* -07/30/03 MJ added user set criterium for stable eigenvalues - corrected error messages in mexfunction() -*/ +/* + * Copyright (C) 2006-2008 Dynare Team + * + * This file is part of Dynare. + * + * Dynare is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * Dynare is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with Dynare. If not, see . + */ + +#include + +#include +#include + +double criterium; + +lapack_int +my_criteria(const double *alphar, const double *alphai, const double *beta) +{ + return ((*alphar * *alphar + *alphai * *alphai) < criterium * *beta * *beta); +} + +void +mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r, double *eval_i, double *info) +{ + lapack_int i_n, i_info, i_sdim, one, lwork; + double *alphar, *alphai, *beta, *work, *par, *pai, *pb, *per, *pei; + double *junk; + lapack_int *bwork; + + one = 1; + i_n = (lapack_int)*n; + alphar = mxCalloc(i_n, sizeof(double)); + alphai = mxCalloc(i_n, sizeof(double)); + beta = mxCalloc(i_n, sizeof(double)); + lwork = 16*i_n+16; + work = mxCalloc(lwork, sizeof(double)); + bwork = mxCalloc(i_n, sizeof(lapack_int)); + /* made necessary by bug in Lapack */ + junk = mxCalloc(i_n*i_n, sizeof(double)); + + dgges("N", "V", "S", my_criteria, &i_n, a, &i_n, b, &i_n, &i_sdim, alphar, alphai, beta, junk, &i_n, z, &i_n, work, &lwork, bwork, &i_info); + + *sdim = i_sdim; + *info = i_info; + + par = alphar; + pai = alphai; + pb = beta; + pei = eval_i; + for (per = eval_r; per <= &eval_r[i_n-1]; ++per) + { + *per = *par / *pb; + *pei = *pai / *pb; + ++par; + ++pai; + ++pb; + ++pei; + } +} + +/* MATLAB interface */ +void +mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + unsigned int m1, n1, m2, n2; + double *s, *t, *z, *sdim, *eval_r, *eval_i, *info, *a, *b; + double n; + + /* Check for proper number of arguments */ + + if (nrhs < 2 || nrhs > 3) + { + mexErrMsgTxt("MJDGGES: two or three input arguments are required."); + } + else if (nlhs > 6) + { + mexErrMsgTxt("MJDGGES: too many output arguments."); + } + + /* Check that A and B are real matrices of the same dimension.*/ + + m1 = mxGetM(prhs[0]); + n1 = mxGetN(prhs[0]); + m2 = mxGetM(prhs[1]); + n2 = mxGetN(prhs[1]); + if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) + || !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) + || (m1 != n1) || (m2 != n1) || (m2 != n2)) + { + mexErrMsgTxt("MJDGGES requires two square real matrices of the same dimension."); + } + + /* Create a matrix for the return argument */ + plhs[0] = mxCreateDoubleMatrix(n1, n1, mxREAL); + plhs[1] = mxCreateDoubleMatrix(n1, n1, mxREAL); + plhs[2] = mxCreateDoubleMatrix(n1, n1, mxREAL); + plhs[3] = mxCreateDoubleMatrix(1, 1, mxREAL); + plhs[4] = mxCreateDoubleMatrix(n1, 1, mxCOMPLEX); + plhs[5] = mxCreateDoubleMatrix(1, 1, mxREAL); + + /* Assign pointers to the various parameters */ + s = mxGetPr(plhs[0]); + t = mxGetPr(plhs[1]); + z = mxGetPr(plhs[2]); + sdim = mxGetPr(plhs[3]); + eval_r = mxGetPr(plhs[4]); + eval_i = mxGetPi(plhs[4]); + info = mxGetPr(plhs[5]); + + a = mxGetPr(prhs[0]); + b = mxGetPr(prhs[1]); + + /* set criterium for stable eigenvalues */ + if (nrhs == 3) + { + criterium = *mxGetPr(prhs[2]); + } + else + { + criterium = 1+1e-6; + } + + /* keep a and b intact */ + memcpy(s, a, sizeof(double)*n1*n1); + memcpy(t, b, sizeof(double)*n1*n1); + + n = n1; + + /* Do the actual computations in a subroutine */ + mjdgges(s, t, z, &n, sdim, eval_r, eval_i, info); + + return; + +} + +/* + 07/30/03 MJ added user set criterium for stable eigenvalues + corrected error messages in mexfunction() +*/