diff --git a/matlab/model_info.m b/matlab/model_info.m
index c6b92f87b..211ea95d9 100644
--- a/matlab/model_info.m
+++ b/matlab/model_info.m
@@ -34,7 +34,7 @@ if(isfield(M_,'block_structure'))
end;
for j=1:size_block
if(j==1)
- fprintf('| %3d (%4d) | %10d | %30s | %14d | %-6d %24s |\n',i,M_.block_structure.block(i).num,size_block,Sym_type(M_.block_structure.block(i).Simulation_Type),M_.block_structure.block(i).equation(j),M_.block_structure.block(i).variable(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
+ fprintf('| %10d | %10d | %30s | %14d | %-6d %24s |\n',i,size_block,Sym_type(M_.block_structure.block(i).Simulation_Type),M_.block_structure.block(i).equation(j),M_.block_structure.block(i).variable(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
else
fprintf('| %10s | %10s | %30s | %14d | %-6d %24s |\n','','','',M_.block_structure.block(i).equation(j),M_.block_structure.block(i).variable(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
end;
diff --git a/matlab/read_data_.m b/matlab/read_data_.m
index 99ffe7faa..7cf9541e8 100644
--- a/matlab/read_data_.m
+++ b/matlab/read_data_.m
@@ -43,7 +43,7 @@ if size(oo_.endo_simul,2) < M_.maximum_lag+M_.maximum_lead+options_.periods
positions = [positions ; strmatch(chopped,M_.endo_names,'exact')];
end
Values=fscanf(fid,'%f',inf);
- Values=reshape(Values,M_.endo_nbr,size(Values,1)/M_.endo_nbr);
+ Values=reshape(Values,M_.orig_endo_nbr,size(Values,1)/M_.orig_endo_nbr);
oo_.endo_simul=Values(positions,:);
fclose(fid);
end
diff --git a/matlab/simul.m b/matlab/simul.m
index fb7d3c942..ba80313dc 100644
--- a/matlab/simul.m
+++ b/matlab/simul.m
@@ -63,8 +63,8 @@ end
if options_.block && ~options_.bytecode && (options_.stack_solve_algo == 0 || options_.stack_solve_algo == 5)
error('SIMUL: for the moment, you must use stack_solve_algo={1,2,3,4} when using block without bytecode option')
end
-if options_.bytecode && options_.stack_solve_algo ~= 5
- error('SIMUL: for the moment, you must use stack_solve_algo=5 with bytecode option')
+if options_.bytecode && (options_.stack_solve_algo ~= 1 && options_.stack_solve_algo ~= 2 && options_.stack_solve_algo ~= 3 && options_.stack_solve_algo ~= 4 && options_.stack_solve_algo ~= 5)
+ error('SIMUL: for the moment, you must use stack_solve_algo= 1, 2, 3, 4 or 5 with bytecode option')
end
if exist('OCTAVE_VERSION') && options_.stack_solve_algo == 2
diff --git a/matlab/steady_.m b/matlab/steady_.m
index 7e9422178..74195c931 100644
--- a/matlab/steady_.m
+++ b/matlab/steady_.m
@@ -29,7 +29,8 @@ function steady_()
% along with Dynare. If not, see .
global M_ oo_ it_ options_
-if options_.bytecode && options_.solve_algo ~= 5
+if options_.bytecode && ...
+ (options_.solve_algo ~= 1 && options_.solve_algo ~= 2 && options_.solve_algo ~= 3 && options_.solve_algo ~= 4 && options_.solve_algo ~= 5)
error('STEADY: for the moment, you must use solve_algo=5 with bytecode option')
end
if ~options_.bytecode && options_.solve_algo == 5
@@ -75,7 +76,9 @@ if options_.steadystate_flag
check1 = 1;
end
elseif options_.block && options_.bytecode
- [residuals, check1] = bytecode('evaluate','static');
+ [residuals, check1] = bytecode('evaluate','static',oo_.steady_state,...
+ [oo_.exo_steady_state; ...
+ oo_.exo_det_steady_state], M_.params, 1);
else
check1 = 0;
check1 = max(abs(feval([M_.fname '_static'],...
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index 666765033..5555cfdf1 100755
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -38,7 +38,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
double *direction_arg, int y_size_arg,
int nb_row_x_arg, int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
int maxit_arg_, double solve_tolf_arg, int size_of_direction_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg,
- string &filename_arg, int minimal_solving_periods_arg)
+ string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg)
{
params = params_arg;
y = y_arg;
@@ -64,6 +64,8 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
T = NULL;
error_not_printed = true;
minimal_solving_periods = minimal_solving_periods_arg;
+ stack_solve_algo = stack_solve_algo_arg;
+ solve_algo = solve_algo_arg;
}
string
@@ -105,7 +107,7 @@ Interpreter::add_underscore_to_fpe(string str)
string
-Interpreter::get_variable(SymbolType variable_type, int variable_num)
+Interpreter::get_variable(SymbolType variable_type, unsigned int variable_num)
{
ostringstream res;
#ifndef DEBUG_EX
@@ -364,6 +366,9 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate)
dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2();
dvar3 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2();
break;
+ default:
+ mexPrintf("Dérivatives %d not implemented yet\n", it_code->first);
+ mexErrMsgTxt("end of bytecode\n");
}
break;
case FLDV:
@@ -1098,15 +1103,31 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate)
}
void
-Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
+Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state)
{
- int var, lag = 0, op;
+ int var = 0, lag = 0, op;
+ unsigned int eq, pos_col;
ostringstream tmp_out;
double v1, v2, v3;
bool go_on = true;
double ll;
-
+ double rr;
+ double *jacob = NULL, *jacob_other_endo = NULL, *jacob_exo = NULL, *jacob_exo_det = NULL;
EQN_block = block_num;
+#ifdef DEBUG
+ mexPrintf("compute_block_time\n");
+#endif
+ #ifndef DEBUG_EX
+ if (evaluate && !steady_state)
+ {
+ jacob = mxGetPr(jacobian_block[block_num]);
+ mexPrintf("jacobian_block[%d]=%x\n",block_num, jacobian_block[block_num]);
+ jacob_other_endo = mxGetPr(jacobian_other_endo_block[block_num]);
+ jacob_exo = mxGetPr(jacobian_exo_block[block_num]);
+ jacob_exo_det = mxGetPr(jacobian_det_exo_block[block_num]);
+ }
+ #endif
+
//feclearexcept (FE_ALL_EXCEPT);
while (go_on)
{
@@ -1114,41 +1135,74 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
switch (it_code->first)
{
case FNUMEXPR:
+#ifdef DEBUG
+ mexPrintf("FNUMEXPR\n");
+#endif
it_code_expr = it_code;
switch (((FNUMEXPR_ *) it_code->second)->get_expression_type())
{
case TemporaryTerm:
+#ifdef DEBUG
+ mexPrintf("TemporaryTerm\n");
+#endif
EQN_type = TemporaryTerm;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
break;
case ModelEquation:
+#ifdef DEBUG
+ mexPrintf("ModelEquation\n");
+#endif
EQN_type = ModelEquation;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
break;
case FirstEndoDerivative:
+#ifdef DEBUG
+ mexPrintf("FirstEndoDerivative\n");
+#endif
EQN_type = FirstEndoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1();
break;
+ case FirstOtherEndoDerivative:
+#ifdef DEBUG
+ mexPrintf("FirstOtherEndoDerivative\n");
+#endif
+ EQN_type = FirstOtherEndoDerivative;
+ EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
+ EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
+ EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1();
+ break;
case FirstExoDerivative:
+#ifdef DEBUG
+ mexPrintf("FirstExoDerivative\n");
+#endif
EQN_type = FirstExoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1();
break;
case FirstExodetDerivative:
+#ifdef DEBUG
+ mexPrintf("FirstExodetDerivative\n");
+#endif
EQN_type = FirstExodetDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1();
break;
case FirstParamDerivative:
+#ifdef DEBUG
+ mexPrintf("FirstParamDerivative\n");
+#endif
EQN_type = FirstParamDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
break;
case SecondEndoDerivative:
+#ifdef DEBUG
+ mexPrintf("SecondEndoDerivative\n");
+#endif
EQN_type = SecondEndoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@@ -1157,6 +1211,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2();
break;
case SecondExoDerivative:
+#ifdef DEBUG
+ mexPrintf("SecondExoDerivative\n");
+#endif
EQN_type = SecondExoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@@ -1165,6 +1222,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2();
break;
case SecondExodetDerivative:
+#ifdef DEBUG
+ mexPrintf("SecondExodetDerivative\n");
+#endif
EQN_type = SecondExodetDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@@ -1173,12 +1233,18 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2();
break;
case SecondParamDerivative:
+#ifdef DEBUG
+ mexPrintf("SecondParamDerivative\n");
+#endif
EQN_type = SecondParamDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2();
break;
case ThirdEndoDerivative:
+#ifdef DEBUG
+ mexPrintf("ThirdEndoDerivative\n");
+#endif
EQN_type = ThirdEndoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@@ -1189,6 +1255,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3();
break;
case ThirdExoDerivative:
+#ifdef DEBUG
+ mexPrintf("ThirdExoDerivative\n");
+#endif
EQN_type = ThirdExoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@@ -1199,6 +1268,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3();
break;
case ThirdExodetDerivative:
+#ifdef DEBUG
+ mexPrintf("ThirdExodetDerivative\n");
+#endif
EQN_type = ThirdExodetDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@@ -1209,6 +1281,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3();
break;
case ThirdParamDerivative:
+#ifdef DEBUG
+ mexPrintf("ThirdParamDerivative\n");
+#endif
EQN_type = ThirdParamDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@@ -1223,14 +1298,18 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
{
case eParameter:
var = ((FLDV_ *) it_code->second)->get_pos();
- Stack.push(params[var]);
#ifdef DEBUG
+ mexPrintf("FLDV Param[var=%d]\n",var);
tmp_out << " params[" << var << "](" << params[var] << ")";
#endif
+ Stack.push(params[var]);
break;
case eEndogenous:
var = ((FLDV_ *) it_code->second)->get_pos();
lag = ((FLDV_ *) it_code->second)->get_lead_lag();
+#ifdef DEBUG
+ mexPrintf("FLDV y[var=%d, lag=%d, it_=%d], y_size=%d evaluate=%d\n",var, lag, it_, y_size, evaluate);
+#endif
if (evaluate)
Stack.push(ya[(it_+lag)*y_size+var]);
else
@@ -1242,10 +1321,11 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
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
+ mexPrintf("FLDV x[var=%d, lag=%d, it_=%d], nb_row_x=%d evaluate=%d\n",var, lag, it_, nb_row_x, evaluate);
tmp_out << " x[" << it_+lag << ", " << var << "](" << x[it_+lag+var*nb_row_x] << ")";
#endif
+ Stack.push(x[it_+lag+var*nb_row_x]);
break;
case eExogenousDet:
var = ((FLDV_ *) it_code->second)->get_pos();
@@ -1268,18 +1348,16 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
{
case eParameter:
var = ((FLDSV_ *) it_code->second)->get_pos();
- Stack.push(params[var]);
#ifdef DEBUG
+ mexPrintf("FLDSV Param[var=%d]\n",var);
tmp_out << " params[" << var << "](" << params[var] << ")";
#endif
+ Stack.push(params[var]);
break;
case eEndogenous:
var = ((FLDSV_ *) it_code->second)->get_pos();
- if (evaluate)
- Stack.push(ya[var]);
- else
- Stack.push(y[var]);
#ifdef DEBUG
+ mexPrintf("FLDSV y[var=%d]\n",var);
tmp_out << " y[" << var << "](" << y[var] << ")";
if(var<0 || var>= y_size)
{
@@ -1287,16 +1365,24 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
mexErrMsgTxt("End of bytecode");
}
#endif
+ if (evaluate)
+ Stack.push(ya[var]);
+ else
+ Stack.push(y[var]);
break;
case eExogenous:
var = ((FLDSV_ *) it_code->second)->get_pos();
- Stack.push(x[var]);
#ifdef DEBUG
+ mexPrintf("FLDSV x[var=%d]\n",var);
tmp_out << " x[" << var << "](" << x[var] << ")";
#endif
+ Stack.push(x[var]);
break;
case eExogenousDet:
var = ((FLDSV_ *) it_code->second)->get_pos();
+#ifdef DEBUG
+ mexPrintf("FLDSV xd[var=%d]\n",var);
+#endif
Stack.push(x[var]);
break;
case eModelLocalVariable:
@@ -1316,28 +1402,28 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
case eParameter:
var = ((FLDVS_ *) it_code->second)->get_pos();
#ifdef DEBUG
- mexPrintf("params[%d]=%f\n", var, params[var]);
+ mexPrintf("params[%d]\n", var);
#endif
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("FLDVS steady_y[%d]\n", var);
#endif
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("FLDVS x[%d] \n", var);
#endif
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("FLDVS xd[%d]\n", var);
#endif
Stack.push(x[var]);
break;
@@ -1364,6 +1450,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
//load a temporary variable in the processor
var = ((FLDST_ *) it_code->second)->get_pos();
#ifdef DEBUG
+ mexPrintf("FLDST T[%d]\n",var);
tmp_out << " T[" << var << "](" << T[var] << ")";
#endif
Stack.push(T[var]);
@@ -1373,6 +1460,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
var = ((FLDU_ *) it_code->second)->get_pos();
var += Per_u_;
#ifdef DEBUG
+ mexPrintf("FLDU u[%d]\n",var);
tmp_out << " u[" << var << "](" << u[var] << ")";
#endif
Stack.push(u[var]);
@@ -1381,6 +1469,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
//load u variable in the processor
var = ((FLDSU_ *) it_code->second)->get_pos();
#ifdef DEBUG
+ mexPrintf("FLDSU u[%d]\n",var);
tmp_out << " u[" << var << "](" << u[var] << ")";
#endif
Stack.push(u[var]);
@@ -1388,10 +1477,16 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
case FLDR:
//load u variable in the processor
var = ((FLDR_ *) it_code->second)->get_pos();
+#ifdef DEBUG
+ mexPrintf("FLDR r[%d]\n",var);
+#endif
Stack.push(r[var]);
break;
case FLDZ:
//load 0 in the processor
+#ifdef DEBUG
+ mexPrintf("FLDZ\n");
+#endif
Stack.push(0.0);
#ifdef DEBUG
tmp_out << " 0";
@@ -1401,6 +1496,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
//load a numerical constant in the processor
ll = ((FLDC_ *) it_code->second)->get_value();
#ifdef DEBUG
+ mexPrintf("FLDC = %f\n",ll);
tmp_out << " " << ll;
#endif
@@ -1412,6 +1508,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
{
case eParameter:
var = ((FSTPV_ *) it_code->second)->get_pos();
+#ifdef DEBUG
+ mexPrintf("FSTPV params[%d]\n",var);
+#endif
params[var] = Stack.top();
Stack.pop();
break;
@@ -1559,6 +1658,63 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
#endif
Stack.pop();
break;
+
+ case FSTPG3:
+ //store in derivative (g) variable from the processor
+#ifdef DEBUG
+ mexPrintf("FSTPG3\n");
+ mexEvalString("drawnow;");
+#endif
+ rr = Stack.top();
+ switch(EQN_type)
+ {
+ case FirstEndoDerivative:
+ eq = ((FSTPG3_ *) it_code->second)->get_row();
+ var = ((FSTPG3_ *) it_code->second)->get_col();
+ lag = ((FSTPG3_ *) it_code->second)->get_lag();
+ pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
+ mexPrintf("jacob[%d(size=%d*pos_col=%d + eq=%d )]=%f\n",eq + size*pos_col, size, pos_col, eq, rr);
+ jacob[eq + size*pos_col] = rr;
+ break;
+ case FirstOtherEndoDerivative:
+ //eq = ((FSTPG3_ *) it_code->second)->get_row();
+ eq = EQN_equation;
+ var = ((FSTPG3_ *) it_code->second)->get_col();
+ lag = ((FSTPG3_ *) it_code->second)->get_lag();
+ pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
+ mexPrintf("jacob_other_endo[%d(size=%d*pos_col=%d + eq=%d)]=%f\n",size*pos_col + eq, size, pos_col, eq, rr);
+ jacob_other_endo[eq + size*pos_col] = rr;
+ break;
+ case FirstExoDerivative:
+ //eq = ((FSTPG3_ *) it_code->second)->get_row();
+ eq = EQN_equation;
+ var = ((FSTPG3_ *) it_code->second)->get_col();
+ lag = ((FSTPG3_ *) it_code->second)->get_lag();
+ pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
+ mexPrintf("jacob_exo[%d(size=%d*pos_col=%dr + eq=%d)]=%f\n",size*pos_col+eq, size, pos_col, eq, rr);
+ jacob_exo[eq + size*pos_col] = rr;
+ break;
+ case FirstExodetDerivative:
+ //eq = ((FSTPG3_ *) it_code->second)->get_row();
+ eq = EQN_equation;
+ var = ((FSTPG3_ *) it_code->second)->get_col();
+ lag = ((FSTPG3_ *) it_code->second)->get_lag();
+ pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
+ mexPrintf("jacob_exo_det[%d(size=%d*pos_col=%dr + eq=%d)]=%f\n",size*pos_col+eq, size, pos_col, eq, rr);
+ jacob_exo_det[eq + size*pos_col] = rr;
+ break;
+ default:
+ mexPrintf("Variable %d not used yet\n", EQN_type);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+#ifdef DEBUG
+ 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();
@@ -1691,7 +1847,6 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
Stack.push(log1(v1, evaluate));
if (isnan(res1))
go_on = false;
-
#ifdef DEBUG
tmp_out << " |log(" << v1 << ")|";
#endif
@@ -1808,10 +1963,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
switch (op)
{
case oNormcdf:
-#ifndef _MSC_VER
- //mexPrintf("normcdf(v1=%f, v2=%f, v3=%f)=%f\n", v1, v2, v3, 0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
- Stack.push(0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
-# ifdef DEBUG
+#ifndef _MSC_VER Stack.push(0.5*(1+erf((v1-v2)/v3/M_SQRT2)));# ifdef DEBUG
tmp_out << " |normcdf(" << v1 << ", " << v2 << ", " << v3 << ")|";
# endif
#else
@@ -1837,10 +1989,30 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
break;
case FENDBLOCK:
//it's the block end
+#ifdef DEBUG
+ mexPrintf("FENDBLOCK\n");
+#endif
go_on = false;
break;
case FENDEQU:
break;
+ case FJMPIFEVAL:
+ if (evaluate)
+ {
+#ifdef DEBUG
+ mexPrintf("FJMPIFEVAL length=%d\n",((FJMPIFEVAL_ *) it_code->second)->get_pos());
+ mexEvalString("drawnow;");
+#endif
+ it_code += ((FJMPIFEVAL_ *) it_code->second)->get_pos()/* - 1*/;
+ }
+ break;
+ case FJMP:
+#ifdef DEBUG
+ mexPrintf("FJMP length=%d\n",((FJMP_ *) it_code->second)->get_pos());
+ mexEvalString("drawnow;");
+#endif
+ it_code += ((FJMP_ *) it_code->second)->get_pos() /*- 1 */;
+ break;
case FOK:
op = ((FOK_ *) it_code->second)->get_arg();
if (Stack.size() > 0)
@@ -1856,18 +2028,23 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
}
it_code++;
}
+#ifdef DEBUG
+ mexPrintf("==> end of compute_block_time Block = %d\n",block_num);
+ mexEvalString("drawnow;");
+#endif
}
void
Interpreter::evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
- const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag, const int Block_List_Max_Lead, const int u_count_int)
+ const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
+ const int Block_List_Max_Lead, const int u_count_int)
{
it_code_type begining;
switch (type)
{
case EVALUATE_FORWARD:
if (steady_state)
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
else
{
begining = it_code;
@@ -1875,7 +2052,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
}
}
break;
@@ -1884,7 +2061,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
r = (double *) mxMalloc(size*sizeof(double));
if (steady_state)
{
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[Block_Contain[j].Variable] += r[j];
}
@@ -1893,10 +2070,9 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
begining = it_code;
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
{
- it_code = begining;
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@@ -1906,11 +2082,11 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
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);
+ Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false, stack_solve_algo, solve_algo);
r = (double *) mxMalloc(size*sizeof(double));
if (steady_state)
{
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[Block_Contain[j].Variable] += r[j];
}
@@ -1921,7 +2097,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@@ -1930,7 +2106,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
break;
case EVALUATE_BACKWARD:
if (steady_state)
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
else
{
begining = it_code;
@@ -1938,7 +2114,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
}
}
break;
@@ -1947,7 +2123,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
r = (double *) mxMalloc(size*sizeof(double));
if (steady_state)
{
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[Block_Contain[j].Variable] += r[j];
}
@@ -1958,7 +2134,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@@ -1968,11 +2144,11 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
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);
+ Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false, stack_solve_algo, solve_algo);
r = (double *) mxMalloc(size*sizeof(double));
if (steady_state)
{
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[Block_Contain[j].Variable] += r[j];
}
@@ -1983,7 +2159,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, true, block_num);
+ compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@@ -1993,7 +2169,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
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);
+ Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true, stack_solve_algo, solve_algo);
u_count = u_count_int*(periods+y_kmax+y_kmin);
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@@ -2002,7 +2178,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
Per_u_ = (it_-y_kmin)*u_count_int;
Per_y_ = it_*y_size;
it_code = begining;
- compute_block_time(Per_u_, true, block_num);
+ compute_block_time(Per_u_, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@@ -2021,11 +2197,18 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
int giter;
bool result = true;
double *y_save;
+ res1 = 0;
+#ifdef DEBUG
+ mexPrintf("simulate_a_block\n");
+#endif
switch (type)
{
case EVALUATE_FORWARD:
+#ifdef DEBUG
+ mexPrintf("EVALUATE_FORWARD\n");
+#endif
if (steady_state)
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
else
{
begining = it_code;
@@ -2033,13 +2216,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
}
}
break;
case EVALUATE_BACKWARD:
+#ifdef DEBUG
+ mexPrintf("EVALUATE_BACKWARD\n");
+#endif
if (steady_state)
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
else
{
begining = it_code;
@@ -2047,11 +2233,14 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
}
}
break;
case SOLVE_FORWARD_SIMPLE:
+#ifdef DEBUG
+ mexPrintf("SOLVE_FORWARD_SIMPLE size=%d\n",size);
+#endif
g1 = (double *) mxMalloc(size*size*sizeof(double));
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@@ -2063,7 +2252,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
double rr;
rr = r[0];
cvg = (fabs(rr) < solve_tolf);
@@ -2099,7 +2288,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
double rr;
if (fabs(1+y[Per_y_+Block_Contain[0].Variable]) > eps)
rr = r[0]/(1+y[Per_y_+Block_Contain[0].Variable]);
@@ -2131,6 +2320,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
mxFree(r);
break;
case SOLVE_BACKWARD_SIMPLE:
+#ifdef DEBUG
+ mexPrintf("SOLVE_BACKWARD_SIMPLE\n");
+#endif
g1 = (double *) mxMalloc(size*size*sizeof(double));
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@@ -2142,7 +2334,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
double rr;
rr = r[0];
cvg = (fabs(rr) < solve_tolf);
@@ -2177,7 +2369,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
double rr;
if (fabs(1+y[Per_y_+Block_Contain[0].Variable]) > eps)
rr = r[0]/(1+y[Per_y_+Block_Contain[0].Variable]);
@@ -2209,8 +2401,11 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
mxFree(r);
break;
case SOLVE_FORWARD_COMPLETE:
+#ifdef DEBUG
+ mexPrintf("SOLVE_FORWARD_COMPLETE\n");
+#endif
fixe_u(&u, u_count_int, u_count_int);
- Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
+ Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false, stack_solve_algo, solve_algo);
g1 = (double *) mxMalloc(size*size*sizeof(double));
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@@ -2231,7 +2426,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
res2 = 0;
res1 = 0;
max_res = 0;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
if (!(isnan(res1) || isinf(res1)))
{
for (i = 0; i < size; i++)
@@ -2253,7 +2448,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo);
iter++;
if (iter > prev_iter)
{
@@ -2273,11 +2468,14 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
it_code = begining;
Per_y_ = it_*y_size;
iter = 0;
- res1 = res2 = max_res = 0; max_res_idx = 0;
+ res1 = 0;
+ res2 = 0;
+ max_res = 0;
+ max_res_idx = 0;
error_not_printed = true;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo);
if (!result)
{
mexPrintf("Convergence not achieved in block %d\n", Block_Count);
@@ -2304,7 +2502,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
res2 = 0;
res1 = 0;
max_res = 0;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
if (!(isnan(res1) || isinf(res1)))
{
for (i = 0; i < size; i++)
@@ -2329,7 +2527,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo);
iter++;
if (iter > prev_iter)
{
@@ -2352,11 +2550,13 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
it_code = begining;
Per_y_ = it_*y_size;
iter = 0;
- res1 = res2 = max_res = 0; max_res_idx = 0;
+ res1 = 0;
+ res2 = 0;
+ max_res = 0; max_res_idx = 0;
error_not_printed = true;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo);
}
}
}
@@ -2368,8 +2568,11 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
mxFree(u);
break;
case SOLVE_BACKWARD_COMPLETE:
+#ifdef DEBUG
+ mexPrintf("SOLVE_BACKWARD_COMPLETE\n");
+#endif
fixe_u(&u, u_count_int, u_count_int);
- Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
+ Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false, stack_solve_algo, solve_algo);
g1 = (double *) mxMalloc(size*size*sizeof(double));
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@@ -2389,7 +2592,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
res2 = 0;
res1 = 0;
max_res = 0;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
if (!(isnan(res1) || isinf(res1)))
{
for (i = 0; i < size; i++)
@@ -2411,7 +2614,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo);
iter++;
if (iter > prev_iter)
{
@@ -2431,11 +2634,13 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
it_code = begining;
Per_y_ = it_*y_size;
iter = 0;
- res1 = res2 = max_res = 0; max_res_idx = 0;
+ res1 = 0;
+ res2 = 0;
+ max_res = 0; max_res_idx = 0;
error_not_printed = true;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo);
if (!result)
{
mexPrintf("Convergence not achieved in block %d\n", Block_Count);
@@ -2462,7 +2667,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
res2 = 0;
res1 = 0;
max_res = 0;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
if (!(isnan(res1) || isinf(res1)))
{
for (i = 0; i < size; i++)
@@ -2487,7 +2692,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo);
iter++;
if (iter > prev_iter)
{
@@ -2510,9 +2715,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
it_code = begining;
Per_y_ = it_*y_size;
error_not_printed = true;
- compute_block_time(0, false, block_num);
+ compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
- result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
+ result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo);
}
}
}
@@ -2525,13 +2730,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
+#ifdef DEBUG
+ mexPrintf("SOLVE_TWO_BOUNDARIES\n");
+#endif
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);
+ Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true, stack_solve_algo, solve_algo);
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));
@@ -2560,7 +2768,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
Per_y_ = it_*y_size;
it_code = begining;
error_not_printed = true;
- compute_block_time(Per_u_, false, block_num);
+ compute_block_time(Per_u_, false, block_num, size, steady_state);
if (isnan(res1) || isinf(res1))
{
memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@@ -2588,7 +2796,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
cvg = (max_res < solve_tolf);
u_count = u_count_saved;
int prev_iter = iter;
- simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number);
+ simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number, stack_solve_algo);
iter++;
if (iter > prev_iter)
{
@@ -2605,13 +2813,15 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
else
{
- res1 = res2 = max_res = 0; max_res_idx = 0;
+ res1 = 0;
+ res2 = 0;
+ 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);
+ compute_block_time(Per_u_, false, block_num, size, steady_state);
for (i = 0; i < size; i++)
{
double rr;
@@ -2626,7 +2836,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
}
cvg = false;
- simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number);
+ simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number, stack_solve_algo);
}
mxFree(r);
mxFree(y_save);
@@ -2644,7 +2854,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
bool
-Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate)
+Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks)
{
bool result = true;
@@ -2686,8 +2896,23 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
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());
+ {
+#ifndef DEBUG_EX
+ if (!steady_state)
+ {
+ jacobian_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_jacob(),mxREAL));
+ mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_nb_col_jacob(), sizeof(jacobian_block[Block_Count]));
+ jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_exo_size(),mxREAL));
+ mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_exo_block.size()=%d fb->get_exo_size()=%d\n",Block_Count, fb->get_size(), fb->get_exo_size(), sizeof(jacobian_exo_block[Block_Count]), fb->get_exo_size());
+ jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_det_exo_size(),mxREAL));
+ mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_det_exo_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_det_exo_size(), sizeof(jacobian_det_exo_block[Block_Count]));
+ jacobian_other_endo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_other_endo_jacob(),mxREAL));
+ mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_other_endo_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_nb_col_other_endo_jacob(), sizeof(jacobian_other_endo_block[Block_Count]));
+ }
+#endif
+ 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());
@@ -2724,13 +2949,14 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
it_code++;
break;
default:
- mexPrintf("Unknown command \n");
+ mexPrintf("Unknown command %d\n",it_code->first);
mexEvalString("drawnow;");
mexErrMsgTxt("End of bytecode");
break;
}
}
mxFree(Init_Code->second);
+ nb_blocks = Block_Count+1;
if (T)
mxFree(T);
return result;
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index fe74981c8..31d7e2fd0 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -41,12 +41,15 @@
using namespace std;
#define pow_ pow
-
-typedef vector >::const_iterator it_code_type;
+typedef vector > code_liste_type;
+typedef code_liste_type::const_iterator it_code_type;
class Interpreter : SparseMatrix
{
private:
+ #ifndef DEBUG_EX
+ vector jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block;
+ #endif
ExpressionType EQN_type;
unsigned int EQN_equation, EQN_block, EQN_block_number;
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
@@ -58,9 +61,9 @@ private:
double log10_1(double a, bool evaluate);
string remove_white(string str);
string add_underscore_to_fpe(string str);
- string get_variable(SymbolType variable_type, int variable_num);
+ string get_variable(SymbolType variable_type, unsigned int variable_num);
string error_location(bool evaluate);
- void compute_block_time(int Per_u_, bool evaluate, int block_num);
+ void compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state);
string print_expression(it_code_type it_code, bool evaluate);
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);
@@ -68,7 +71,7 @@ private:
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;
+ code_liste_type code_liste;
it_code_type it_code;
stack Stack;
int Block_Count, Per_u_, Per_y_;
@@ -82,13 +85,20 @@ private:
int equation, derivative_equation, derivative_variable;
string filename;
int minimal_solving_periods;
+ int stack_solve_algo, solve_algo;
public:
~Interpreter();
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);
+ double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg);
+ bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks);
+ #ifndef DEBUG_EX
+ inline mxArray* get_jacob(int block_num) {return jacobian_block[block_num];};
+ inline mxArray* get_jacob_exo(int block_num) {return jacobian_exo_block[block_num];};
+ inline mxArray* get_jacob_exo_det(int block_num) {return jacobian_det_exo_block[block_num];};
+ inline mxArray* get_jacob_other_endo(int block_num) {return jacobian_other_endo_block[block_num];};
+ #endif
};
#endif
diff --git a/mex/sources/bytecode/Mem_Mngr.cc b/mex/sources/bytecode/Mem_Mngr.cc
index 37162a971..deb3d9bdf 100644
--- a/mex/sources/bytecode/Mem_Mngr.cc
+++ b/mex/sources/bytecode/Mem_Mngr.cc
@@ -27,7 +27,7 @@ Mem_Mngr::Mem_Mngr()
void
Mem_Mngr::Print_heap()
{
- int i;
+ unsigned int i;
mexPrintf("i :");
for (i = 0; i < CHUNK_SIZE; i++)
mexPrintf("%3d ", i);
@@ -61,7 +61,7 @@ Mem_Mngr::init_CHUNK_BLCK_SIZE(int u_count)
NonZeroElem *
Mem_Mngr::mxMalloc_NZE()
{
- long int i;
+ long unsigned int i;
if (!Chunk_Stack.empty()) /*An unused block of memory available inside the heap*/
{
NonZeroElem *p1 = Chunk_Stack.back();
@@ -102,7 +102,7 @@ Mem_Mngr::mxMalloc_NZE()
void
Mem_Mngr::mxFree_NZE(void *pos)
{
- int i;
+ unsigned int i;
size_t gap;
for (i = 0; i < Nb_CHUNK; i++)
{
diff --git a/mex/sources/bytecode/Mem_Mngr.hh b/mex/sources/bytecode/Mem_Mngr.hh
index 3cc7d6ffe..a70f8d431 100644
--- a/mex/sources/bytecode/Mem_Mngr.hh
+++ b/mex/sources/bytecode/Mem_Mngr.hh
@@ -52,8 +52,8 @@ public:
bool swp_f;
private:
v_NonZeroElem Chunk_Stack;
- int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK;
- int CHUNK_heap_pos;
+ unsigned int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK;
+ unsigned int CHUNK_heap_pos;
NonZeroElem **NZE_Mem_add;
NonZeroElem *NZE_Mem;
vector NZE_Mem_Allocated;
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 4af23d9c6..1d8a0511d 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -46,6 +46,8 @@ SparseMatrix::SparseMatrix()
start_compare = 0;
restart = 0;
IM_i.clear();
+ lu_inc_tol = 1e-10;
+ reduced = true;
}
int
@@ -278,7 +280,7 @@ SparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_
}
void
-SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries)
+SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries, int stack_solve_algo, int solve_algo)
{
unsigned int eq, var;
int i, j, lag;
@@ -303,26 +305,57 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
IM_i.clear();
if (two_boundaries)
{
- for (i = 0; i < u_count_init-Size; i++)
+ if (stack_solve_algo == 5)
{
- SaveCode.read(reinterpret_cast(&eq), sizeof(eq));
- SaveCode.read(reinterpret_cast(&var), sizeof(var));
- SaveCode.read(reinterpret_cast(&lag), sizeof(lag));
- SaveCode.read(reinterpret_cast(&j), sizeof(j));
- IM_i[make_pair(make_pair(eq, var), lag)] = j;
+ for (i = 0; i < u_count_init-Size; i++)
+ {
+ SaveCode.read(reinterpret_cast(&eq), sizeof(eq));
+ SaveCode.read(reinterpret_cast(&var), sizeof(var));
+ SaveCode.read(reinterpret_cast(&lag), sizeof(lag));
+ SaveCode.read(reinterpret_cast(&j), sizeof(j));
+ IM_i[make_pair(make_pair(eq, var), lag)] = j;
+ }
+ for (j = 0; j < Size; j++)
+ IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
}
- for (j = 0; j < Size; j++)
- IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
+ else if (stack_solve_algo >= 1 || stack_solve_algo <= 4)
+ {
+ for (i = 0; i < u_count_init-Size; i++)
+ {
+ SaveCode.read(reinterpret_cast(&eq), sizeof(eq));
+ SaveCode.read(reinterpret_cast(&var), sizeof(var));
+ SaveCode.read(reinterpret_cast(&lag), sizeof(lag));
+ SaveCode.read(reinterpret_cast(&j), sizeof(j));
+ IM_i[make_pair(make_pair(var - lag*Size, -lag), eq)] = j;
+ }
+ for (j = 0; j < Size; j++)
+ IM_i[make_pair(make_pair(Size*(periods+y_kmax), 0), j)] = j;
+ }
+
}
else
{
- for (i = 0; i < u_count_init; i++)
+ if ((stack_solve_algo == 5 && !steady_state) || (solve_algo == 5 && steady_state))
{
- SaveCode.read(reinterpret_cast(&eq), sizeof(eq));
- SaveCode.read(reinterpret_cast(&var), sizeof(var));
- SaveCode.read(reinterpret_cast(&lag), sizeof(lag));
- SaveCode.read(reinterpret_cast(&j), sizeof(j));
- IM_i[make_pair(make_pair(eq, var), lag)] = j;
+ for (i = 0; i < u_count_init; i++)
+ {
+ SaveCode.read(reinterpret_cast(&eq), sizeof(eq));
+ SaveCode.read(reinterpret_cast(&var), sizeof(var));
+ SaveCode.read(reinterpret_cast(&lag), sizeof(lag));
+ SaveCode.read(reinterpret_cast(&j), sizeof(j));
+ IM_i[make_pair(make_pair(eq, var), lag)] = j;
+ }
+ }
+ else if ( ((stack_solve_algo >= 1 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 1 || solve_algo <= 4) && steady_state) )
+ {
+ for (i = 0; i < u_count_init; i++)
+ {
+ SaveCode.read(reinterpret_cast(&eq), sizeof(eq));
+ SaveCode.read(reinterpret_cast(&var), sizeof(var));
+ SaveCode.read(reinterpret_cast(&lag), sizeof(lag));
+ SaveCode.read(reinterpret_cast(&j), sizeof(j));
+ IM_i[make_pair(make_pair(var - lag*Size, -lag), eq)] = j;
+ }
}
}
index_vara = (int *) mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int));
@@ -415,7 +448,230 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map, int>, int> &IM)
+SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map, int>, int> &IM, mxArray *A_m, mxArray *b_m)
+{
+ int i, eq, var;
+ double *b = mxGetPr(b_m);
+ if (!b)
+ {
+ mexPrintf("Can't retrieve b vector in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ mwIndex *Ai = mxGetIr(A_m);
+ if (!Ai)
+ {
+ mexPrintf("Can't allocate Ai index vector in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ mwIndex *Aj = mxGetJc(A_m);
+ if (!Aj)
+ {
+ mexPrintf("Can't allocate Aj index vector in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ double *A = mxGetPr(A_m);
+ if (!A)
+ {
+ mexPrintf("Can't retrieve A matrix in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ map, int>, int>::iterator it4;
+ for (i = 0; i < y_size*(periods+y_kmin); i++)
+ ya[i] = y[i];
+#if DEBUG
+ unsigned int max_nze = mxGetNzmax(A_m);
+#endif
+ unsigned int NZE = 0;
+ int last_var = 0;
+ for (i = 0; i < Size; i++)
+ b[i] = u[i];
+ Aj[0] = 0;
+ last_var = 0;
+ it4 = IM.begin();
+ while (it4 != IM.end())
+ {
+ var = it4->first.first.first;
+ if (var != last_var)
+ {
+ Aj[1+last_var ] = NZE;
+ last_var = var;
+ }
+ eq = it4->first.second;
+ int index = it4->second;
+#if DEBUG
+ if (index<0 || index >= u_count_alloc)
+ {
+ mexPrintf("index (%d) out of range for u vector (0)\n",index);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (NZE >= max_nze)
+ {
+ mexPrintf("Exceeds the capacity of A_m sparse matrix in LU (0)\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+#endif
+ A[NZE] = u[index];
+ Ai[NZE] = eq;
+ NZE++;
+#if DEBUG
+ if ((index+lag*u_count_init) < 0 || (index+lag*u_count_init) >= u_count_alloc)
+ {
+ mexPrintf("index (%d) out of range for u vector (1)\n",index+lag*u_count_init);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (eq < 0 || eq >= (Size*periods))
+ {
+ mexPrintf("index (%d) out of range for b vector (0)\n",eq);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (var+Size*(y_kmin+t+lag) < 0 || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
+ {
+ mexPrintf("index (%d) out of range for index_vara vector (0)\n",var+Size*(y_kmin+t+lag));
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (index_vara[var+Size*(y_kmin+t+lag)] < 0 || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
+ {
+ mexPrintf("index (%d) out of range for y vector max=%d (0)\n",index_vara[var+Size*(y_kmin+t+lag)], y_size*(periods+y_kmin+y_kmax));
+ mexErrMsgTxt("end of bytecode\n");
+ }
+#endif
+ it4++;
+ }
+ Aj[Size] = NZE;
+}
+
+
+
+void
+SparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM, mxArray *A_m, mxArray *b_m)
+{
+ int t, i, eq, var, lag, ti_y_kmin, ti_y_kmax;
+ double *b = mxGetPr(b_m);
+ if (!b)
+ {
+ mexPrintf("Can't retrieve b vector in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ mwIndex *Ai = mxGetIr(A_m);
+ if (!Ai)
+ {
+ mexPrintf("Can't allocate Ai index vector in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ mwIndex *Aj = mxGetJc(A_m);
+ if (!Aj)
+ {
+ mexPrintf("Can't allocate Aj index vector in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ double *A = mxGetPr(A_m);
+ if (!A)
+ {
+ mexPrintf("Can't retrieve A matrix in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ map, int>, int>::iterator it4;
+ for (i = 0; i < y_size*(periods+y_kmin); i++)
+ ya[i] = y[i];
+#if DEBUG
+ unsigned int max_nze = mxGetNzmax(A_m);
+#endif
+ unsigned int NZE = 0;
+ int last_var = 0;
+ for (i = 0; i < periods*Size; i++)
+ b[i] = 0;
+ Aj[0] = 0;
+ for (t = 0; t < periods; t++)
+ {
+ last_var = 0;
+ it4 = IM.begin();
+ while (it4 != IM.end())
+ {
+ var = it4->first.first.first;
+ if (var != last_var)
+ {
+ Aj[1+last_var + t * Size] = NZE;
+ last_var = var;
+ }
+ eq = it4->first.second+Size*t;
+ lag = -it4->first.first.second;
+ int index = it4->second+ (t-lag) * u_count_init;
+ if (var < (periods+y_kmax)*Size)
+ {
+ ti_y_kmin = -min(t, y_kmin);
+ ti_y_kmax = min(periods-(t +1 ), y_kmax);
+ int ti_new_y_kmax = min(t, y_kmax);
+ int ti_new_y_kmin = -min(periods-(t+1), y_kmin);
+ if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/
+ {
+#if DEBUG
+ if (index<0 || index >= u_count_alloc)
+ {
+ mexPrintf("index (%d) out of range for u vector (0)\n",index);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (NZE >= max_nze)
+ {
+ mexPrintf("Exceeds the capacity of A_m sparse matrix in LU (0)\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+#endif
+ A[NZE] = u[index];
+ Ai[NZE] = eq - lag * Size;
+ NZE++;
+ }
+ if (lag > ti_y_kmax || lag < ti_y_kmin)
+ {
+#if DEBUG
+ if ((index+lag*u_count_init) < 0 || (index+lag*u_count_init) >= u_count_alloc)
+ {
+ mexPrintf("index (%d) out of range for u vector (1)\n",index+lag*u_count_init);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (eq < 0 || eq >= (Size*periods))
+ {
+ mexPrintf("index (%d) out of range for b vector (0)\n",eq);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (var+Size*(y_kmin+t+lag) < 0 || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
+ {
+ mexPrintf("index (%d) out of range for index_vara vector (0)\n",var+Size*(y_kmin+t+lag));
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (index_vara[var+Size*(y_kmin+t+lag)] < 0 || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
+ {
+ mexPrintf("index (%d) out of range for y vector max=%d (0)\n",index_vara[var+Size*(y_kmin+t+lag)], y_size*(periods+y_kmin+y_kmax));
+ mexErrMsgTxt("end of bytecode\n");
+ }
+#endif
+ b[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]];
+ }
+ }
+ else /* ...and store it in the u vector*/
+ {
+#if DEBUG
+ if (index < 0 || index >= u_count_alloc)
+ {
+ mexPrintf("index (%d) out of range for u vector (2)\n",index);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ if (eq < 0 || eq >= (Size*periods))
+ {
+ mexPrintf("index (%d) out of range for b vector (1)\n",eq);
+ mexErrMsgTxt("end of bytecode\n");
+ }
+#endif
+ b[eq] += u[index];
+ }
+ it4++;
+ }
+ }
+ Aj[Size*periods] = NZE;
+}
+
+
+void
+SparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM)
{
int t, i, eq, var, lag, ti_y_kmin, ti_y_kmax;
double tmp_b = 0.0;
@@ -524,47 +780,6 @@ SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM)
-{
- int t, eq, var, lag, ti_y_kmin, ti_y_kmax;
- double tmp_b = 0.0;
- map, int>, int>::iterator it4;
- //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(it4, ti_y_kmin, ti_y_kmax, eq, var, lag, tmp_b) schedule(dynamic)
- for (t = 0; t < periods; t++)
- {
- ti_y_kmin = -min(t, y_kmin);
- ti_y_kmax = min(periods-(t+1), y_kmax);
- it4 = IM.begin();
- eq = -1;
- while (it4 != IM.end())
- {
- var = it4->first.first.second;
- if (eq != it4->first.first.first+Size*t)
- tmp_b = 0;
- eq = it4->first.first.first+Size*t;
- if (var < (periods+y_kmax)*Size)
- {
- lag = it4->first.second;
- if (lag <= ti_y_kmax && lag >= ti_y_kmin)
- {
- var += Size*t;
- }
- else
- {
- tmp_b += u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
- }
- }
- else
- {
- b[eq] = it4->second+u_count_init*t;
- u[b[eq]] += tmp_b;
- }
- it4++;
- }
- }
-}
-
int
SparseMatrix::Get_u()
{
@@ -619,7 +834,7 @@ SparseMatrix::Print_u()
}
void
-SparseMatrix::End(int Size)
+SparseMatrix::End_GE(int Size)
{
mem_mngr.Free_All();
mxFree(FNZE_R);
@@ -995,7 +1210,7 @@ SparseMatrix::simple_bksub(int it_, int Size, double slowc_l)
}
bool
-SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int Block_number)
+SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int Block_number, int solve_algo)
{
int i, j, k;
int pivj = 0, pivk = 0;
@@ -1005,6 +1220,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
int *pivj_v, *pivk_v, *NR;
int l, N_max;
bool one;
+ mxArray *b_m, *A_m;
Clear_u();
piv_v = (double *) mxMalloc(Size*sizeof(double));
pivj_v = (int *) mxMalloc(Size*sizeof(int));
@@ -1131,230 +1347,209 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
mexPrintf(" abs. error=%.10e \n", double (res1));
mexPrintf("-----------------------------------\n");
}
- 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*/
- double piv = piv_abs = 0;
- int nb_eq = At_Col(i, &first);
- l = 0; N_max = 0;
- one = false;
- piv_abs = 0;
- for (j = 0; j < nb_eq; j++)
- {
- if (!line_done[first->r_index])
- {
- k = first->u_index;
- int jj = first->r_index;
- int NRow_jj = NRow(jj);
- piv_v[l] = u[k];
- double piv_fabs = fabs(u[k]);
- pivj_v[l] = jj;
- pivk_v[l] = k;
- NR[l] = NRow_jj;
- if (NRow_jj == 1 && !one)
+ if (solve_algo == 5)
+ Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
+ else
+ {
+ b_m = mxCreateDoubleMatrix(Size,1,mxREAL);
+ if (!b_m)
+ {
+ mexPrintf("Can't allocate b_m matrix in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ A_m = mxCreateSparse(Size, Size, IM_i.size()*2, mxREAL);
+ if (!A_m)
+ {
+ mexPrintf("Can't allocate A_m matrix in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ Init_Matlab_Sparse_Simple(Size, IM_i, A_m, b_m);
+ }
+ if (solve_algo == 5)
+ {
+ NonZeroElem **bc;
+ bc = (NonZeroElem **) mxMalloc(Size*sizeof(*bc));
+ for (i = 0; i < Size; i++)
+ {
+ /*finding the max-pivot*/
+ double piv = piv_abs = 0;
+ int nb_eq = At_Col(i, &first);
+ l = 0; N_max = 0;
+ one = false;
+ piv_abs = 0;
+ for (j = 0; j < nb_eq; j++)
+ {
+ if (!line_done[first->r_index])
{
- one = true;
- piv_abs = piv_fabs;
- N_max = NRow_jj;
- }
- if (!one)
- {
- if (piv_fabs > piv_abs)
- piv_abs = piv_fabs;
- if (NRow_jj > N_max)
- N_max = NRow_jj;
- }
- else
- {
- if (NRow_jj == 1)
+ k = first->u_index;
+ int jj = first->r_index;
+ int NRow_jj = NRow(jj);
+
+ piv_v[l] = u[k];
+ double piv_fabs = fabs(u[k]);
+ pivj_v[l] = jj;
+ pivk_v[l] = k;
+ NR[l] = NRow_jj;
+ if (NRow_jj == 1 && !one)
+ {
+ one = true;
+ piv_abs = piv_fabs;
+ N_max = NRow_jj;
+ }
+ if (!one)
{
if (piv_fabs > piv_abs)
piv_abs = piv_fabs;
if (NRow_jj > N_max)
N_max = NRow_jj;
}
+ else
+ {
+ if (NRow_jj == 1)
+ {
+ if (piv_fabs > piv_abs)
+ piv_abs = piv_fabs;
+ if (NRow_jj > N_max)
+ N_max = NRow_jj;
+ }
+ }
+ l++;
+ }
+ first = first->NZE_C_N;
+ }
+ if (piv_abs < eps)
+ {
+ if (Block_number > 1)
+ mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1);
+ else
+ mexPrintf("Error: singular system in Simulate_NG\n");
+ mexEvalString("drawnow;");
+ mxFree(piv_v);
+ mxFree(pivj_v);
+ mxFree(pivk_v);
+ mxFree(NR);
+ mxFree(bc);
+ if (steady_state)
+ return false;
+ else
+ {
+ mexEvalString("st=fclose('all');clear all;");
+ filename += " stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
+ }
+ double markovitz = 0, markovitz_max = -9e70;
+ int NR_max = 0;
+ if (!one)
+ {
+ for (j = 0; j < l; j++)
+ {
+ if (N_max > 0 && NR[j] > 0)
+ {
+ if (fabs(piv_v[j]) > 0)
+ {
+ if (markowitz_c > 0)
+ markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ else
+ markovitz = fabs(piv_v[j])/piv_abs;
+ }
+ else
+ markovitz = 0;
+ }
+ else
+ markovitz = fabs(piv_v[j])/piv_abs;
+ if (markovitz > markovitz_max)
+ {
+ piv = piv_v[j];
+ pivj = pivj_v[j]; //Line number
+ pivk = pivk_v[j]; //positi
+ markovitz_max = markovitz;
+ NR_max = NR[j];
+ }
}
- l++;
}
- first = first->NZE_C_N;
- }
- if (piv_abs < eps)
- {
- if (Block_number > 1)
- mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1);
- else
- mexPrintf("Error: singular system in Simulate_NG\n");
- mexEvalString("drawnow;");
- mxFree(piv_v);
- mxFree(pivj_v);
- mxFree(pivk_v);
- mxFree(NR);
- mxFree(bc);
- if (steady_state)
- return false;
else
{
- mexEvalString("st=fclose('all');clear all;");
- filename += " stopped";
- mexErrMsgTxt(filename.c_str());
- }
- }
- double markovitz = 0, markovitz_max = -9e70;
- int NR_max = 0;
- if (!one)
- {
- for (j = 0; j < l; j++)
- {
- if (N_max > 0 && NR[j] > 0)
+ for (j = 0; j < l; j++)
{
- if (fabs(piv_v[j]) > 0)
+ if (N_max > 0 && NR[j] > 0)
{
- if (markowitz_c > 0)
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ if (fabs(piv_v[j]) > 0)
+ {
+ if (markowitz_c > 0)
+ markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ else
+ markovitz = fabs(piv_v[j])/piv_abs;
+ }
else
- markovitz = fabs(piv_v[j])/piv_abs;
+ markovitz = 0;
}
else
- markovitz = 0;
- }
- else
- markovitz = fabs(piv_v[j])/piv_abs;
- if (markovitz > markovitz_max)
- {
- piv = piv_v[j];
- pivj = pivj_v[j]; //Line number
- pivk = pivk_v[j]; //positi
- markovitz_max = markovitz;
- NR_max = NR[j];
- }
- }
- }
- else
- {
- for (j = 0; j < l; j++)
- {
- if (N_max > 0 && NR[j] > 0)
- {
- if (fabs(piv_v[j]) > 0)
+ markovitz = fabs(piv_v[j])/piv_abs;
+ if (NR[j] == 1)
{
- if (markowitz_c > 0)
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
- else
- markovitz = fabs(piv_v[j])/piv_abs;
+ piv = piv_v[j];
+ pivj = pivj_v[j]; //Line number
+ pivk = pivk_v[j]; //positi
+ markovitz_max = markovitz;
+ NR_max = NR[j];
}
- else
- markovitz = 0;
- }
- else
- markovitz = fabs(piv_v[j])/piv_abs;
- if (NR[j] == 1)
- {
- piv = piv_v[j];
- pivj = pivj_v[j]; //Line number
- pivk = pivk_v[j]; //positi
- markovitz_max = markovitz;
- NR_max = NR[j];
}
}
- }
- /*if (fabs(piv) < eps)
- mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n",NR_max, N_max, piv, piv_abs, markovitz_max);
- if (NR_max == 0)
- mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n",piv, markovitz_max);*/
- pivot[i] = pivj;
- pivotk[i] = pivk;
- pivotv[i] = piv;
- line_done[pivj] = true;
+ pivot[i] = pivj;
+ pivotk[i] = pivk;
+ pivotv[i] = piv;
+ line_done[pivj] = true;
- /*divide all the non zeros elements of the line pivj by the max_pivot*/
- int nb_var = At_Row(pivj, &first);
- for (j = 0; j < nb_var; j++)
- {
- u[first->u_index] /= piv;
- first = first->NZE_R_N;
- }
- u[b[pivj]] /= piv;
- /*substract the elements on the non treated lines*/
- nb_eq = At_Col(i, &first);
- NonZeroElem *first_piva;
- int nb_var_piva = At_Row(pivj, &first_piva);
- int nb_eq_todo = 0;
- for (j = 0; j < nb_eq && first; j++)
- {
- if (!line_done[first->r_index])
- bc[nb_eq_todo++] = first;
- first = first->NZE_C_N;
- }
- //#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];
-
- int nb_var_piv = nb_var_piva;
- NonZeroElem *first_piv = first_piva;
- NonZeroElem *first_sub;
- int nb_var_sub = At_Row(row, &first_sub);
- int l_sub = 0, l_piv = 0;
- int sub_c_index = first_sub->c_index, piv_c_index = first_piv->c_index;
- while (l_sub < nb_var_sub || l_piv < nb_var_piv)
+ /*divide all the non zeros elements of the line pivj by the max_pivot*/
+ int nb_var = At_Row(pivj, &first);
+ for (j = 0; j < nb_var; j++)
{
- if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv))
+ u[first->u_index] /= piv;
+ first = first->NZE_R_N;
+ }
+ u[b[pivj]] /= piv;
+ /*substract the elements on the non treated lines*/
+ nb_eq = At_Col(i, &first);
+ NonZeroElem *first_piva;
+ int nb_var_piva = At_Row(pivj, &first_piva);
+ int nb_eq_todo = 0;
+ for (j = 0; j < nb_eq && first; j++)
+ {
+ if (!line_done[first->r_index])
+ bc[nb_eq_todo++] = first;
+ first = first->NZE_C_N;
+ }
+ //#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];
+
+ int nb_var_piv = nb_var_piva;
+ NonZeroElem *first_piv = first_piva;
+ NonZeroElem *first_sub;
+ int nb_var_sub = At_Row(row, &first_sub);
+ int l_sub = 0, l_piv = 0;
+ int sub_c_index = first_sub->c_index, piv_c_index = first_piv->c_index;
+ while (l_sub < nb_var_sub || l_piv < nb_var_piv)
{
- first_sub = first_sub->NZE_R_N;
- if (first_sub)
- sub_c_index = first_sub->c_index;
- else
- sub_c_index = Size;
- l_sub++;
- }
- else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub)
- {
- int tmp_u_count = Get_u();
- Insert(row, first_piv->c_index, tmp_u_count, 0);
- u[tmp_u_count] = -u[first_piv->u_index]*first_elem;
- first_piv = first_piv->NZE_R_N;
- if (first_piv)
- piv_c_index = first_piv->c_index;
- else
- piv_c_index = Size;
- l_piv++;
- }
- else
- {
- if (i == sub_c_index)
+ if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv))
{
- firsta = first;
- first_suba = first_sub->NZE_R_N;
- Delete(first_sub->r_index, first_sub->c_index);
- first = firsta->NZE_C_N;
- first_sub = first_suba;
- if (first_sub)
- sub_c_index = first_sub->c_index;
- else
- sub_c_index = Size;
- l_sub++;
- first_piv = first_piv->NZE_R_N;
- if (first_piv)
- piv_c_index = first_piv->c_index;
- else
- piv_c_index = Size;
- l_piv++;
- }
- else
- {
- u[first_sub->u_index] -= u[first_piv->u_index]*first_elem;
first_sub = first_sub->NZE_R_N;
if (first_sub)
sub_c_index = first_sub->c_index;
else
sub_c_index = Size;
l_sub++;
+ }
+ else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub)
+ {
+ int tmp_u_count = Get_u();
+ Insert(row, first_piv->c_index, tmp_u_count, 0);
+ u[tmp_u_count] = -u[first_piv->u_index]*first_elem;
first_piv = first_piv->NZE_R_N;
if (first_piv)
piv_c_index = first_piv->c_index;
@@ -1362,22 +1557,66 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
piv_c_index = Size;
l_piv++;
}
+ else
+ {
+ if (i == sub_c_index)
+ {
+ firsta = first;
+ first_suba = first_sub->NZE_R_N;
+ Delete(first_sub->r_index, first_sub->c_index);
+ first = firsta->NZE_C_N;
+ first_sub = first_suba;
+ if (first_sub)
+ sub_c_index = first_sub->c_index;
+ else
+ sub_c_index = Size;
+ l_sub++;
+ first_piv = first_piv->NZE_R_N;
+ if (first_piv)
+ piv_c_index = first_piv->c_index;
+ else
+ piv_c_index = Size;
+ l_piv++;
+ }
+ else
+ {
+ u[first_sub->u_index] -= u[first_piv->u_index]*first_elem;
+ first_sub = first_sub->NZE_R_N;
+ if (first_sub)
+ sub_c_index = first_sub->c_index;
+ else
+ sub_c_index = Size;
+ l_sub++;
+ first_piv = first_piv->NZE_R_N;
+ if (first_piv)
+ piv_c_index = first_piv->c_index;
+ else
+ piv_c_index = Size;
+ l_piv++;
+ }
+ }
}
+ u[b[row]] -= u[b[pivj]]*first_elem;
}
- 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_GE(Size);
+ mxFree(piv_v);
+ mxFree(pivj_v);
+ mxFree(pivk_v);
+ mxFree(NR);
+ mxFree(bc);
}
- 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);
- mxFree(bc);
+ else if (solve_algo == 1 || solve_algo == 4)
+ Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc);
+ else if (solve_algo == 2)
+ Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck);
+ else if (solve_algo == 3)
+ Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck);
return true;
}
@@ -1456,7 +1695,7 @@ void
SparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b)
{
const double epsilon = 1e-10;
- Init(periods, y_kmin, y_kmax, Size, IM_i);
+ Init_GE(periods, y_kmin, y_kmax, Size, IM_i);
NonZeroElem *first;
int cal_y = y_kmin*Size;
mexPrintf(" ");
@@ -1492,8 +1731,170 @@ SparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size,
mexErrMsgTxt(filename.c_str());
}
+void
+SparseMatrix::Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l)
+{
+ int n = mxGetM(A_m);
+ mxArray *z;
+ mxArray *rhs[2];
+ rhs[0] = A_m;
+ rhs[1] = b_m;
+ mexCallMATLAB(1,&z,2, rhs, "mldivide");
+ double *res = mxGetPr(z);
+ for (int i = 0; i < n; i++)
+ {
+ int eq = index_vara[i+Size*y_kmin];
+ double yy = - (res[i] + y[eq]);
+ direction[eq] = yy;
+ y[eq] += slowc_l * yy;
+ }
+ mxDestroyArray(A_m);
+ mxDestroyArray(b_m);
+ mxDestroyArray(z);
+}
+
+void
+SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block)
+{
+ int n = mxGetM(A_m);
+ /*[L1, U1]=luinc(g1a,luinc_tol);*/
+ mxArray *lhs0[2];
+ mxArray *rhs0[2];
+ rhs0[0] = A_m;
+ rhs0[1] = mxCreateDoubleScalar(lu_inc_tol);
+ mexCallMATLAB(2, lhs0, 2, rhs0, "luinc");
+ mxArray *L1 = lhs0[0];
+ mxArray *U1 = lhs0[1];
+ /*[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);*/
+ mxArray *rhs[7];
+ rhs[0] = A_m;
+ rhs[1] = b_m;
+ rhs[2] = mxCreateDoubleScalar(Size);
+ rhs[3] = mxCreateDoubleScalar(1e-6);
+ rhs[4] = mxCreateDoubleScalar(n);
+ rhs[5] = L1;
+ rhs[6] = U1;
+ mxArray *lhs[2];
+ mexCallMATLAB(2,lhs, 7, rhs, "gmres");
+ mxArray *z = lhs[0];
+ mxArray *flag = lhs[1];
+ double *flag1 = mxGetPr(flag);
+ mxDestroyArray(rhs0[1]);
+ mxDestroyArray(rhs[2]);
+ mxDestroyArray(rhs[3]);
+ mxDestroyArray(rhs[4]);
+ mxDestroyArray(rhs[5]);
+ mxDestroyArray(rhs[6]);
+ if (*flag1 > 0 || reduced)
+ {
+ ostringstream tmp;
+ if (*flag1 == 1)
+ {
+ tmp << "Error in simul: No convergence inside GMRES, in block " << block;
+ mexWarnMsgTxt(tmp.str().c_str());
+ }
+ else if (*flag1 == 2)
+ {
+ tmp << "Error in simul: Preconditioner is ill-conditioned, in block " << block;
+ mexWarnMsgTxt(tmp.str().c_str());
+ }
+ else if (*flag1 == 3)
+ {
+ tmp << "Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block " << block;
+ mexWarnMsgTxt(tmp.str().c_str());
+ }
+ lu_inc_tol /= 10;
+ reduced = false;
+ }
+ else
+ {
+ double *res = mxGetPr(z);
+ for (int i = 0; i < n; i++)
+ {
+ int eq = index_vara[i+Size*y_kmin];
+ double yy = - (res[i] + y[eq]);
+ direction[eq] = yy;
+ y[eq] += slowc * yy;
+ }
+ }
+ mxDestroyArray(A_m);
+ mxDestroyArray(b_m);
+ mxDestroyArray(z);
+ mxDestroyArray(flag);
+}
+
+
+void
+SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block)
+{
+ int n = mxGetM(A_m);
+ /*[L1, U1]=luinc(g1a,luinc_tol);*/
+ mxArray *lhs0[2];
+ mxArray *rhs0[2];
+ rhs0[0] = A_m;
+ rhs0[1] = mxCreateDoubleScalar(lu_inc_tol);
+ mexCallMATLAB(2, lhs0, 2, rhs0, "luinc");
+ mxArray *L1 = lhs0[0];
+ mxArray *U1 = lhs0[1];
+ /*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
+ mxArray *rhs[6];
+ rhs[0] = A_m;
+ rhs[1] = b_m;
+ rhs[2] = mxCreateDoubleScalar(1e-6);
+ rhs[3] = mxCreateDoubleScalar(n);
+ rhs[4] = L1;
+ rhs[5] = U1;
+ mxArray *lhs[2];
+ mexCallMATLAB(2,lhs, 6, rhs, "bicgstab");
+ mxArray *z = lhs[0];
+ mxArray *flag = lhs[1];
+ double *flag1 = mxGetPr(flag);
+ mxDestroyArray(rhs0[1]);
+ mxDestroyArray(rhs[2]);
+ mxDestroyArray(rhs[3]);
+ mxDestroyArray(rhs[4]);
+ mxDestroyArray(rhs[5]);
+ if (*flag1 > 0 || reduced)
+ {
+ ostringstream tmp;
+ if (*flag1 == 1)
+ {
+ tmp << "Error in simul: No convergence inside BiCGStab, in block " << block;
+ mexWarnMsgTxt(tmp.str().c_str());
+ }
+ else if (*flag1 == 2)
+ {
+ tmp << "Error in simul: Preconditioner is ill-conditioned, in block " << block;
+ mexWarnMsgTxt(tmp.str().c_str());
+ }
+ else if (*flag1 == 3)
+ {
+ tmp << "Error in simul: BiCGStab stagnated (Two consecutive iterates were the same.), in block " << block;
+ mexWarnMsgTxt(tmp.str().c_str());
+ }
+ lu_inc_tol /= 10;
+ reduced = false;
+ }
+ else
+ {
+ double *res = mxGetPr(z);
+ for (int i = 0; i < n; i++)
+ {
+ int eq = index_vara[i+Size*y_kmin];
+ double yy = - (res[i] + y[eq]);
+ direction[eq] = yy;
+ y[eq] += slowc * yy;
+ }
+ }
+ mxDestroyArray(A_m);
+ mxDestroyArray(b_m);
+ mxDestroyArray(z);
+ mxDestroyArray(flag);
+}
+
+
int
-SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number)
+SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number, int stack_solve_algo)
{
/*Triangularisation at each period of a block using a simple gaussian Elimination*/
t_save_op_s *save_op_s;
@@ -1513,6 +1914,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
clock_t t1 = clock();
nop1 = 0;
error_not_printed = true;
+ mxArray *b_m, *A_m;
if (iter > 0)
{
mexPrintf("Sim : %f ms\n", (1000.0*(double (clock())-double (time00)))/double (CLOCKS_PER_SEC));
@@ -1565,7 +1967,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
filename += " stopped";
mexErrMsgTxt(filename.c_str());
}
- if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)))
+ if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)) && (stack_solve_algo == 1 || stack_solve_algo == 5))
{
if (try_at_iteration == 0)
{
@@ -1583,6 +1985,16 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
glambda2 = res2;
try_at_iteration ++;
+ if (slowc_save<=0.1)
+ {
+ for (i = 0; i < y_size*(periods+y_kmin); i++)
+ y[i] = ya[i]+direction[i];
+ g0 = res2;
+ gp0 = -res2;
+ try_at_iteration = 0;
+ iter--;
+ return 0;
+ }
}
else
{
@@ -1596,7 +2008,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
iter--;
return 0;
}
- /*if (isnan(res1) || isinf(res1))
+ if (isnan(res1) || isinf(res1))
{
if (iter == 0)
{
@@ -1629,48 +2041,75 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
y[i] = ya[i]+slowc_save*direction[i];
iter--;
return (0);
- }*/
- u_count += u_count_init;
- if (alt_symbolic && alt_symbolic_count < alt_symbolic_count_max)
- {
- mexPrintf("Pivoting method will be applied only to the first periods.\n");
- alt_symbolic = false;
- symbolic = true;
- markowitz_c = markowitz_c_s;
- alt_symbolic_count++;
}
- if (((res1/res1a-1) > -0.3) && symbolic && iter > 0)
+ u_count += u_count_init;
+ if (stack_solve_algo == 5)
{
- if (restart > 2)
+ if (alt_symbolic && alt_symbolic_count < alt_symbolic_count_max)
{
- mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied to all periods.\n");
- symbolic = false;
- alt_symbolic = true;
- markowitz_c_s = markowitz_c;
- markowitz_c = 0;
+ mexPrintf("Pivoting method will be applied only to the first periods.\n");
+ alt_symbolic = false;
+ symbolic = true;
+ markowitz_c = markowitz_c_s;
+ alt_symbolic_count++;
+ }
+ if (((res1/res1a-1) > -0.3) && symbolic && iter > 0)
+ {
+ if (restart > 2)
+ {
+ mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied to all periods.\n");
+ symbolic = false;
+ alt_symbolic = true;
+ markowitz_c_s = markowitz_c;
+ markowitz_c = 0;
+ }
+ else
+ {
+ mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied for a longer period.\n");
+ start_compare = min(tbreak_g, periods);
+ restart++;
+ }
}
else
{
- mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied for a longer period.\n");
- start_compare = min(tbreak_g, periods);
- restart++;
+ start_compare = max(y_kmin, minimal_solving_periods);
+ restart = 0;
}
}
- else
- {
- start_compare = max(y_kmin, minimal_solving_periods);
- restart = 0;
- }
res1a = res1;
if (print_it)
{
+ if (iter == 0)
+ {
+ switch(stack_solve_algo)
+ {
+ case 1:
+ mexPrintf("MODEL SIMULATION: (method=Sparse LU)\n");
+ break;
+ case 2:
+ mexPrintf("MODEL SIMULATION: (method=GMRES)\n");
+ break;
+ case 3:
+ mexPrintf("MODEL SIMULATION: (method=BiCGStab)\n");
+ break;
+ case 4:
+ mexPrintf("MODEL SIMULATION: (method=Sparse LU & optimal path length)\n");
+ break;
+ case 5:
+ mexPrintf("MODEL SIMULATION: (method=ByteCode own solver)\n");
+ break;
+ default:
+ mexPrintf("MODEL SIMULATION: (method=Unknown - %d - )\n", stack_solve_algo);
+ }
+ }
mexPrintf("-----------------------------------\n");
mexPrintf(" Simulate iteration no %d \n", iter+1);
mexPrintf(" max. error=%.10e \n", double (max_res));
mexPrintf(" sqr. error=%.10e \n", double (res2));
mexPrintf(" abs. error=%.10e \n", double (res1));
mexPrintf("-----------------------------------\n");
+ mexEvalString("drawnow;");
}
if (cvg)
{
@@ -1678,250 +2117,173 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
else
{
- 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));
- for (int t = 0; t < periods; t++)
+ if (stack_solve_algo == 5)
+ Init_GE(periods, y_kmin, y_kmax, Size, IM_i);
+ else
{
- if (record && symbolic)
+ b_m = mxCreateDoubleMatrix(periods*Size,1,mxREAL);
+ if (!b_m)
{
- if (save_op)
- {
- mxFree(save_op);
- save_op = NULL;
- }
- save_op = (int *) mxMalloc(nop*sizeof(int));
- nopa = nop;
+ mexPrintf("Can't allocate b_m matrix in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
}
- nop = 0;
- Clear_u();
- int ti = t*Size;
- for (i = ti; i < Size+ti; i++)
+ A_m = mxCreateSparse(periods*Size, periods*Size, IM_i.size()* periods*2, mxREAL);
+ if (!A_m)
{
- /*finding the max-pivot*/
- double piv = piv_abs = 0;
- int nb_eq = At_Col(i, 0, &first);
- if ((symbolic && t <= start_compare) || !symbolic)
+ mexPrintf("Can't allocate A_m matrix in LU solver\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ Init_Matlab_Sparse(periods, y_kmin, y_kmax, Size, IM_i, A_m, b_m);
+ }
+ if (stack_solve_algo == 5)
+ {
+ 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));
+ for (int t = 0; t < periods; t++)
+ {
+ if (record && symbolic)
{
- int l = 0, N_max = 0;
- bool one = false;
- piv_abs = 0;
- for (j = 0; j < nb_eq; j++)
+ if (save_op)
{
- if (!line_done[first->r_index])
+ mxFree(save_op);
+ save_op = NULL;
+ }
+ save_op = (int *) mxMalloc(nop*sizeof(int));
+ nopa = nop;
+ }
+ nop = 0;
+ Clear_u();
+ int ti = t*Size;
+ for (i = ti; i < Size+ti; i++)
+ {
+ /*finding the max-pivot*/
+ double piv = piv_abs = 0;
+ int nb_eq = At_Col(i, 0, &first);
+ if ((symbolic && t <= start_compare) || !symbolic)
+ {
+ int l = 0, N_max = 0;
+ bool one = false;
+ piv_abs = 0;
+ for (j = 0; j < nb_eq; j++)
{
- k = first->u_index;
- int jj = first->r_index;
- int NRow_jj = NRow(jj);
- piv_v[l] = u[k];
- double piv_fabs = fabs(u[k]);
- pivj_v[l] = jj;
- pivk_v[l] = k;
- NR[l] = NRow_jj;
- if (NRow_jj == 1 && !one)
+ if (!line_done[first->r_index])
{
- one = true;
- piv_abs = piv_fabs;
- N_max = NRow_jj;
- }
- if (!one)
- {
- if (piv_fabs > piv_abs)
- piv_abs = piv_fabs;
- if (NRow_jj > N_max)
- N_max = NRow_jj;
- }
- else
- {
- if (NRow_jj == 1)
+ k = first->u_index;
+ int jj = first->r_index;
+ int NRow_jj = NRow(jj);
+ piv_v[l] = u[k];
+ double piv_fabs = fabs(u[k]);
+ pivj_v[l] = jj;
+ pivk_v[l] = k;
+ NR[l] = NRow_jj;
+ if (NRow_jj == 1 && !one)
+ {
+ one = true;
+ piv_abs = piv_fabs;
+ N_max = NRow_jj;
+ }
+ if (!one)
{
if (piv_fabs > piv_abs)
piv_abs = piv_fabs;
if (NRow_jj > N_max)
N_max = NRow_jj;
}
- }
- l++;
- }
- first = first->NZE_C_N;
- }
- double markovitz = 0, markovitz_max = -9e70;
- int NR_max = 0;
- if (!one)
- {
- for (j = 0; j < l; j++)
- {
- if (N_max > 0 && NR[j] > 0)
- {
- if (fabs(piv_v[j]) > 0)
+ else
{
- if (markowitz_c > 0)
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ if (NRow_jj == 1)
+ {
+ if (piv_fabs > piv_abs)
+ piv_abs = piv_fabs;
+ if (NRow_jj > N_max)
+ N_max = NRow_jj;
+ }
+ }
+ l++;
+ }
+ first = first->NZE_C_N;
+ }
+ double markovitz = 0, markovitz_max = -9e70;
+ int NR_max = 0;
+ if (!one)
+ {
+ for (j = 0; j < l; j++)
+ {
+ if (N_max > 0 && NR[j] > 0)
+ {
+ if (fabs(piv_v[j]) > 0)
+ {
+ if (markowitz_c > 0)
+ markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ else
+ markovitz = fabs(piv_v[j])/piv_abs;
+ }
else
- markovitz = fabs(piv_v[j])/piv_abs;
+ markovitz = 0;
}
else
- markovitz = 0;
- }
- else
- markovitz = fabs(piv_v[j])/piv_abs;
- if (markovitz > markovitz_max)
- {
- piv = piv_v[j];
- pivj = pivj_v[j]; //Line number
- pivk = pivk_v[j]; //positi
- markovitz_max = markovitz;
- NR_max = NR[j];
+ markovitz = fabs(piv_v[j])/piv_abs;
+ if (markovitz > markovitz_max)
+ {
+ piv = piv_v[j];
+ pivj = pivj_v[j]; //Line number
+ pivk = pivk_v[j]; //positi
+ markovitz_max = markovitz;
+ NR_max = NR[j];
+ }
}
}
+ else
+ {
+ for (j = 0; j < l; j++)
+ {
+ if (N_max > 0 && NR[j] > 0)
+ {
+ if (fabs(piv_v[j]) > 0)
+ {
+ if (markowitz_c > 0)
+ markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
+ else
+ markovitz = fabs(piv_v[j])/piv_abs;
+ }
+ else
+ markovitz = 0;
+ }
+ else
+ markovitz = fabs(piv_v[j])/piv_abs;
+ if (NR[j] == 1)
+ {
+ piv = piv_v[j];
+ pivj = pivj_v[j]; //Line number
+ pivk = pivk_v[j]; //positi
+ markovitz_max = markovitz;
+ NR_max = NR[j];
+ }
+ }
+ }
+ if (fabs(piv) < eps)
+ mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n",NR_max, N_max, piv, piv_abs, markovitz_max);
+ if (NR_max == 0)
+ mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n",piv, markovitz_max);
+ pivot[i] = pivj;
+ pivot_save[i] = pivj;
+ pivotk[i] = pivk;
+ pivotv[i] = piv;
}
else
{
- for (j = 0; j < l; j++)
- {
- if (N_max > 0 && NR[j] > 0)
- {
- if (fabs(piv_v[j]) > 0)
- {
- if (markowitz_c > 0)
- markovitz = exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double (NR[j])/double (N_max)));
- else
- markovitz = fabs(piv_v[j])/piv_abs;
- }
- else
- markovitz = 0;
- }
- else
- markovitz = fabs(piv_v[j])/piv_abs;
- if (NR[j] == 1)
- {
- piv = piv_v[j];
- pivj = pivj_v[j]; //Line number
- pivk = pivk_v[j]; //positi
- markovitz_max = markovitz;
- NR_max = NR[j];
- }
- }
+ pivj = pivot[i-Size]+Size;
+ pivot[i] = pivj;
+ At_Pos(pivj, i, &first);
+ pivk = first->u_index;
+ piv = u[pivk];
+ piv_abs = fabs(piv);
}
- if (fabs(piv) < eps)
- mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n",NR_max, N_max, piv, piv_abs, markovitz_max);
- if (NR_max == 0)
- mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n",piv, markovitz_max);
- pivot[i] = pivj;
- pivot_save[i] = pivj;
- pivotk[i] = pivk;
- pivotv[i] = piv;
- }
- else
- {
- pivj = pivot[i-Size]+Size;
- pivot[i] = pivj;
- At_Pos(pivj, i, &first);
- pivk = first->u_index;
- piv = u[pivk];
- piv_abs = fabs(piv);
- }
- line_done[pivj] = true;
- if (symbolic)
- {
- if (record)
- {
- if (nop+1 >= nopa)
- {
- nopa = long (mem_increasing_factor*(double)nopa);
- save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
- }
- save_op_s = (t_save_op_s *)(&(save_op[nop]));
- save_op_s->operat = IFLD;
- save_op_s->first = pivk;
- save_op_s->lag = 0;
- }
- nop += 2;
- }
- if (piv_abs < eps)
- {
- if (Block_number>1)
- mexPrintf("Error: singular system in Simulate_NG1 in block %d\n", blck);
- else
- mexPrintf("Error: singular system in Simulate_NG1\n");
- mexEvalString("drawnow;");
- filename += " stopped";
- mexErrMsgTxt(filename.c_str());
- }
- /*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));
- for (j = 0; j < nb_var; j++)
- {
- bb[j] = first;
- first = first->NZE_R_N;
- }
-
- for (j = 0; j < nb_var; j++)
- {
- first = bb[j];
- u[first->u_index] /= piv;
- if (symbolic)
- {
- if (record)
- {
- if (nop+j*2+1 >= nopa)
- {
- nopa = long (mem_increasing_factor*(double)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->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;
- if (symbolic)
- {
- if (record)
- {
- if (nop+1 >= nopa)
- {
- nopa = long (mem_increasing_factor*(double)nopa);
- save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
- }
- save_op_s = (t_save_op_s *)(&(save_op[nop]));
- save_op_s->operat = IFDIV;
- save_op_s->first = b[pivj];
- save_op_s->lag = 0;
- }
- nop += 2;
- }
- /*substract the elements on the non treated lines*/
- nb_eq = At_Col(i, &first);
- NonZeroElem *first_piva;
- int nb_var_piva = At_Row(pivj, &first_piva);
-
- NonZeroElem **bc;
- bc = (NonZeroElem **) mxMalloc(nb_eq*sizeof(first));
- int nb_eq_todo = 0;
- for (j = 0; j < nb_eq && first; j++)
- {
- if (!line_done[first->r_index])
- bc[nb_eq_todo++] = first;
- first = first->NZE_C_N;
- }
- //#pragma omp parallel for num_threads(2) shared(nb_var_piva, first_piva, nopa, nop, save_op, record)
- for (j = 0; j < nb_eq_todo; j++)
- {
- t_save_op_s *save_op_s_l;
- first = bc[j];
- int row = first->r_index;
- double first_elem = u[first->u_index];
+ line_done[pivj] = true;
if (symbolic)
{
if (record)
@@ -1931,123 +2293,158 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
nopa = long (mem_increasing_factor*(double)nopa);
save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
}
- 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);
+ save_op_s = (t_save_op_s *)(&(save_op[nop]));
+ save_op_s->operat = IFLD;
+ save_op_s->first = pivk;
+ save_op_s->lag = 0;
}
nop += 2;
}
-
- int nb_var_piv = nb_var_piva;
- NonZeroElem *first_piv = first_piva;
- NonZeroElem *first_sub;
- int nb_var_sub = At_Row(row, &first_sub);
- int l_sub = 0;
- int l_piv = 0;
- int sub_c_index = first_sub->c_index;
- int piv_c_index = first_piv->c_index;
- int tmp_lag = first_sub->lag_index;
- while (l_sub < nb_var_sub || l_piv < nb_var_piv)
+ if (piv_abs < eps)
{
- if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv))
- {
- //There is no nonzero element at row pivot for this column=> Nothing to do for the current element got to next column
- first_sub = first_sub->NZE_R_N;
- if (first_sub)
- sub_c_index = first_sub->c_index;
- else
- sub_c_index = Size*periods;
- l_sub++;
- }
- else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub)
- {
- // There is an nonzero element at row pivot but not at the current row=> insert a negative element in the current row
- tmp_u_count = Get_u();
- lag = first_piv->c_index/Size-row/Size;
- //#pragma omp critical
- {
- Insert(row, first_piv->c_index, tmp_u_count, lag);
- }
- u[tmp_u_count] = -u[first_piv->u_index]*first_elem;
- if (symbolic)
- {
- if (record)
- {
- if (nop+2 >= nopa)
- {
- nopa = long (mem_increasing_factor*(double)nopa);
- save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
- }
- 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;
- save_op_s_l->lag = max(first_piv->lag_index, abs(tmp_lag));
- }
- nop += 3;
- }
- first_piv = first_piv->NZE_R_N;
- if (first_piv)
- piv_c_index = first_piv->c_index;
- else
- piv_c_index = Size*periods;
- l_piv++;
- }
- else /*first_sub->c_index==first_piv->c_index*/
- {
- if (i == sub_c_index)
- {
+ if (Block_number>1)
+ mexPrintf("Error: singular system in Simulate_NG1 in block %d\n", blck);
+ else
+ mexPrintf("Error: singular system in Simulate_NG1\n");
+ mexEvalString("drawnow;");
+ filename += " stopped";
+ mexErrMsgTxt(filename.c_str());
+ }
+ /*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));
+ for (j = 0; j < nb_var; j++)
+ {
+ bb[j] = first;
+ first = first->NZE_R_N;
+ }
- //#pragma omp barrier
- //#pragma omp single
- //#pragma omp critical
- {
- NonZeroElem *firsta = first;
- NonZeroElem *first_suba = first_sub->NZE_R_N;
- Delete(first_sub->r_index, first_sub->c_index);
- first = firsta->NZE_C_N;
- first_sub = first_suba;
- }
-
- if (first_sub)
- sub_c_index = first_sub->c_index;
- else
- sub_c_index = Size*periods;
- l_sub++;
- first_piv = first_piv->NZE_R_N;
- if (first_piv)
- piv_c_index = first_piv->c_index;
- else
- piv_c_index = Size*periods;
- l_piv++;
- }
- else
+ for (j = 0; j < nb_var; j++)
+ {
+ first = bb[j];
+ u[first->u_index] /= piv;
+ if (symbolic)
+ {
+ if (record)
{
- u[first_sub->u_index] -= u[first_piv->u_index]*first_elem;
- if (symbolic)
+ if (nop+j*2+1 >= nopa)
{
- if (record)
- {
- if (nop+3 >= nopa)
- {
- nopa = long (mem_increasing_factor*(double)nopa);
- save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
- }
- 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;
- save_op_s_l->lag = max(abs(tmp_lag), first_piv->lag_index);
- }
- nop += 3;
+ nopa = long (mem_increasing_factor*(double)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->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;
+ if (symbolic)
+ {
+ if (record)
+ {
+ if (nop+1 >= nopa)
+ {
+ nopa = long (mem_increasing_factor*(double)nopa);
+ save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
+ }
+ save_op_s = (t_save_op_s *)(&(save_op[nop]));
+ save_op_s->operat = IFDIV;
+ save_op_s->first = b[pivj];
+ save_op_s->lag = 0;
+ }
+ nop += 2;
+ }
+ /*substract the elements on the non treated lines*/
+ nb_eq = At_Col(i, &first);
+ NonZeroElem *first_piva;
+ int nb_var_piva = At_Row(pivj, &first_piva);
+
+ NonZeroElem **bc;
+ bc = (NonZeroElem **) mxMalloc(nb_eq*sizeof(first));
+ int nb_eq_todo = 0;
+ for (j = 0; j < nb_eq && first; j++)
+ {
+ if (!line_done[first->r_index])
+ bc[nb_eq_todo++] = first;
+ first = first->NZE_C_N;
+ }
+ //#pragma omp parallel for num_threads(2) shared(nb_var_piva, first_piva, nopa, nop, save_op, record)
+ for (j = 0; j < nb_eq_todo; j++)
+ {
+ t_save_op_s *save_op_s_l;
+ first = bc[j];
+ int row = first->r_index;
+ double first_elem = u[first->u_index];
+ if (symbolic)
+ {
+ if (record)
+ {
+ if (nop+1 >= nopa)
+ {
+ nopa = long (mem_increasing_factor*(double)nopa);
+ save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
+ }
+ 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);
+ }
+ nop += 2;
+ }
+
+ int nb_var_piv = nb_var_piva;
+ NonZeroElem *first_piv = first_piva;
+ NonZeroElem *first_sub;
+ int nb_var_sub = At_Row(row, &first_sub);
+ int l_sub = 0;
+ int l_piv = 0;
+ int sub_c_index = first_sub->c_index;
+ int piv_c_index = first_piv->c_index;
+ int tmp_lag = first_sub->lag_index;
+ while (l_sub < nb_var_sub || l_piv < nb_var_piv)
+ {
+ if (l_sub < nb_var_sub && (sub_c_index < piv_c_index || l_piv >= nb_var_piv))
+ {
+ //There is no nonzero element at row pivot for this column=> Nothing to do for the current element got to next column
first_sub = first_sub->NZE_R_N;
if (first_sub)
sub_c_index = first_sub->c_index;
else
sub_c_index = Size*periods;
l_sub++;
+ }
+ else if (sub_c_index > piv_c_index || l_sub >= nb_var_sub)
+ {
+ // There is an nonzero element at row pivot but not at the current row=> insert a negative element in the current row
+ tmp_u_count = Get_u();
+ lag = first_piv->c_index/Size-row/Size;
+ //#pragma omp critical
+ {
+ Insert(row, first_piv->c_index, tmp_u_count, lag);
+ }
+ u[tmp_u_count] = -u[first_piv->u_index]*first_elem;
+ if (symbolic)
+ {
+ if (record)
+ {
+ if (nop+2 >= nopa)
+ {
+ nopa = long (mem_increasing_factor*(double)nopa);
+ save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
+ }
+ 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;
+ save_op_s_l->lag = max(first_piv->lag_index, abs(tmp_lag));
+ }
+ nop += 3;
+ }
first_piv = first_piv->NZE_R_N;
if (first_piv)
piv_c_index = first_piv->c_index;
@@ -2055,108 +2452,169 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
piv_c_index = Size*periods;
l_piv++;
}
- }
- }
- u[b[row]] -= u[b[pivj]]*first_elem;
-
- if (symbolic)
- {
- if (record)
- {
- if (nop+3 >= nopa)
+ else /*first_sub->c_index==first_piv->c_index*/
{
- nopa = long (mem_increasing_factor*(double)nopa);
- save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
+ if (i == sub_c_index)
+ {
+ NonZeroElem *firsta = first;
+ NonZeroElem *first_suba = first_sub->NZE_R_N;
+ Delete(first_sub->r_index, first_sub->c_index);
+ first = firsta->NZE_C_N;
+ first_sub = first_suba;
+ if (first_sub)
+ sub_c_index = first_sub->c_index;
+ else
+ sub_c_index = Size*periods;
+ l_sub++;
+ first_piv = first_piv->NZE_R_N;
+ if (first_piv)
+ piv_c_index = first_piv->c_index;
+ else
+ piv_c_index = Size*periods;
+ l_piv++;
+ }
+ else
+ {
+ u[first_sub->u_index] -= u[first_piv->u_index]*first_elem;
+ if (symbolic)
+ {
+ if (record)
+ {
+ if (nop+3 >= nopa)
+ {
+ nopa = long (mem_increasing_factor*(double)nopa);
+ save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
+ }
+ 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;
+ save_op_s_l->lag = max(abs(tmp_lag), first_piv->lag_index);
+ }
+ nop += 3;
+ }
+ first_sub = first_sub->NZE_R_N;
+ if (first_sub)
+ sub_c_index = first_sub->c_index;
+ else
+ sub_c_index = Size*periods;
+ l_sub++;
+ first_piv = first_piv->NZE_R_N;
+ if (first_piv)
+ piv_c_index = first_piv->c_index;
+ else
+ piv_c_index = Size*periods;
+ l_piv++;
+ }
}
- 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];
- save_op_s_l->lag = abs(tmp_lag);
}
- nop += 3;
- }
- }
- mxFree(bc);
- }
- if (symbolic)
- {
- if (record && (nop == nop1))
- {
- if (save_opa && save_opaa)
- {
- if (compare(save_op, save_opa, save_opaa, t, periods, nop, Size))
+ u[b[row]] -= u[b[pivj]]*first_elem;
+
+ if (symbolic)
{
- tbreak = t;
- tbreak_g = tbreak;
- break;
+ if (record)
+ {
+ if (nop+3 >= nopa)
+ {
+ nopa = long (mem_increasing_factor*(double)nopa);
+ save_op = (int *) mxRealloc(save_op, nopa*sizeof(int));
+ }
+ 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];
+ save_op_s_l->lag = abs(tmp_lag);
+ }
+ nop += 3;
}
}
- if (save_opa)
- {
- if (save_opaa)
- {
- mxFree(save_opaa);
- save_opaa = NULL;
- }
- save_opaa = (int *) mxMalloc(nop1*sizeof(int));
- memcpy(save_opaa, save_opa, nop1*sizeof(int));
- }
- if (save_opa)
- {
- mxFree(save_opa);
- save_opa = NULL;
- }
- save_opa = (int *) mxMalloc(nop*sizeof(int));
- memcpy(save_opa, save_op, nop*sizeof(int));
+ mxFree(bc);
}
- else
+ if (symbolic)
{
- if (nop == nop1)
- record = true;
- else
+ if (record && (nop == nop1))
{
- record = false;
+ if (save_opa && save_opaa)
+ {
+ if (compare(save_op, save_opa, save_opaa, t, periods, nop, Size))
+ {
+ tbreak = t;
+ tbreak_g = tbreak;
+ break;
+ }
+ }
+ if (save_opa)
+ {
+ if (save_opaa)
+ {
+ mxFree(save_opaa);
+ save_opaa = NULL;
+ }
+ save_opaa = (int *) mxMalloc(nop1*sizeof(int));
+ memcpy(save_opaa, save_opa, nop1*sizeof(int));
+ }
if (save_opa)
{
mxFree(save_opa);
save_opa = NULL;
}
- if (save_opaa)
+ save_opa = (int *) mxMalloc(nop*sizeof(int));
+ memcpy(save_opa, save_op, nop*sizeof(int));
+ }
+ else
+ {
+ if (nop == nop1)
+ record = true;
+ else
{
- mxFree(save_opaa);
- save_opaa = NULL;
+ record = false;
+ if (save_opa)
+ {
+ mxFree(save_opa);
+ save_opa = NULL;
+ }
+ if (save_opaa)
+ {
+ mxFree(save_opaa);
+ save_opaa = NULL;
+ }
}
}
+ nop2 = nop1;
+ nop1 = nop;
}
- nop2 = nop1;
- nop1 = nop;
}
- }
- mxFree(piv_v);
- mxFree(pivj_v);
- mxFree(pivk_v);
- mxFree(NR);
- }
- nop_all += nop;
- if (symbolic)
- {
- if (save_op)
- mxFree(save_op);
- if (save_opa)
- mxFree(save_opa);
- if (save_opaa)
- mxFree(save_opaa);
- }
+ mxFree(piv_v);
+ mxFree(pivj_v);
+ mxFree(pivk_v);
+ mxFree(NR);
+ nop_all += nop;
+ if (symbolic)
+ {
+ if (save_op)
+ mxFree(save_op);
+ if (save_opa)
+ mxFree(save_opa);
+ if (save_opaa)
+ mxFree(save_opaa);
+ }
- /*The backward substitution*/
- double slowc_lbx = slowc, res1bx;
- for (i = 0; i < y_size*(periods+y_kmin); i++)
- ya[i] = y[i];
- slowc_save = slowc;
- res1bx = bksub(tbreak, last_period, Size, slowc_lbx);
- t01 = clock();
- End(Size);
+ /*The backward substitution*/
+ double slowc_lbx = slowc, res1bx;
+ for (i = 0; i < y_size*(periods+y_kmin); i++)
+ ya[i] = y[i];
+ slowc_save = slowc;
+ res1bx = bksub(tbreak, last_period, Size, slowc_lbx);
+ t01 = clock();
+ End_GE(Size);
+ }
+ else if (stack_solve_algo == 1 || stack_solve_algo == 4)
+ Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc);
+ else if (stack_solve_algo == 2)
+ Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck);
+ else if (stack_solve_algo == 3)
+ Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck);
+ }
if (print_it)
{
clock_t t2 = clock();
diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh
index c7dd24a83..57b3e5ddd 100644
--- a/mex/sources/bytecode/SparseMatrix.hh
+++ b/mex/sources/bytecode/SparseMatrix.hh
@@ -33,7 +33,10 @@ using namespace std;
#ifdef _MSC_VER
# include
-
+#include
+#define M_SQRT2 1.4142135623730950488016887242097
+#define M_PI 3.1415926535897932384626433832795
+#define erf(x) boost::math::erf(x)
extern unsigned long _nan[2];
extern double NAN;
@@ -100,18 +103,23 @@ 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, int Block_number);
- 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, int Block_number);
+ 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, int Block_number, int stack_solve_algo);
+ 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, int Block_number, int solve_algo);
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_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries, int stack_solve_algo, int solve_algo);
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);
double g0, gp0, glambda2, try_at_iteration;
+
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 Init_GE(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM);
+ void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map, int>, int> &IM, mxArray *A_m, mxArray *b_m);
+ void Init_Matlab_Sparse_Simple(int Size, map, int>, int> &IM, mxArray *A_m, mxArray *b_m);
void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map, int>, int> &IM);
- void End(int Size);
+ void End_GE(int Size);
+ void Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l);
+ void Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block);
+ void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block);
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
@@ -173,6 +181,8 @@ private:
protected:
int u_count_alloc, u_count_alloc_save;
double *u, *y, *ya;
+ vector jac;
+ double *jcb;
double res1, res2, max_res, max_res_idx;
double slowc, slowc_save, prev_slowc_save, markowitz_c;
int y_kmin, y_kmax, y_size, periods, y_decal;
@@ -184,6 +194,8 @@ protected:
int restart;
bool error_not_printed;
double g_lambda1, g_lambda2, gp_0;
+ double lu_inc_tol;
+ bool reduced;
};
#endif
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index 0b0b52810..f37dc2d0f 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -54,15 +54,8 @@ main(int argc, const char *argv[])
FILE *fid;
bool steady_state = false;
bool evaluate = false;
- printf("argc=%d\n", argc);
- /*fexcept_t *flagp;
- flagp = (fexcept_t*) mxMalloc(sizeof(fexcept_t));
- if (fegetexceptflag(flagp, FE_ALL_EXCEPT))
- mexPrintf("fegetexceptflag failed\n");
- if (fesetexceptflag(flagp,FE_INVALID | FE_DIVBYZERO))
- mexPrintf("fesetexceptflag failed\n");
- mxFree(flagp);
- feclearexcept (FE_ALL_EXCEPT);*/
+ bool block = false;
+
if (argc < 2)
{
mexPrintf("model filename expected\n");
@@ -72,7 +65,7 @@ main(int argc, const char *argv[])
float f_tmp;
ostringstream tmp_out("");
tmp_out << argv[1] << "_options.txt";
- cout << tmp_out.str().c_str() << "\n";
+
int nb_params;
int i, row_y, col_y, row_x, col_x;
double *yd, *xd;
@@ -81,6 +74,7 @@ main(int argc, const char *argv[])
string file_name(argv[1]);
+ int count_array_argument = 0;
for (i = 2; i < argc; i++)
{
if (Get_Argument(argv[i]) == "static")
@@ -89,6 +83,8 @@ main(int argc, const char *argv[])
steady_state = false;
else if (Get_Argument(argv[i]) == "evaluate")
evaluate = true;
+ else if (Get_Argument(argv[i]) == "block")
+ block = true;
else
{
mexPrintf("Unknown argument : ");
@@ -99,6 +95,7 @@ main(int argc, const char *argv[])
mexErrMsgTxt(f.c_str());
}
}
+
fid = fopen(tmp_out.str().c_str(), "r");
int periods = 1;
if (!steady_state)
@@ -143,10 +140,12 @@ main(int argc, const char *argv[])
params[i] = f_tmp;
}
fclose(fid);
+
yd = (double *) malloc(row_y*col_y*sizeof(yd[0]));
xd = (double *) malloc(row_x*col_x*sizeof(xd[0]));
tmp_out.str("");
tmp_out << argv[1] << "_oo.txt";
+
fid = fopen(tmp_out.str().c_str(), "r");
for (i = 0; i < col_y*row_y; i++)
{
@@ -193,7 +192,8 @@ main(int argc, const char *argv[])
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods);
string f(file_name);
- interprete.compute_blocks(f, f, steady_state, evaluate);
+ int nb_blocks = 0;
+ interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks);
clock_t t1 = clock();
if (!evaluate)
mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
@@ -235,29 +235,56 @@ void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mxArray *M_, *oo_, *options_;
- int i, row_y, col_y, row_x, col_x, nb_row_xd;
+ mxArray *block_structur = NULL;
+ int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd;
int y_kmin = 0, y_kmax = 0, y_decal = 0, periods = 1;
double *pind;
double *direction;
bool steady_state = false;
bool evaluate = false;
- /*fexcept_t *flagp;
- flagp = (fexcept_t*) mxMalloc(sizeof(fexcept_t));
- if(fegetexceptflag(flagp, FE_ALL_EXCEPT))
- mexPrintf("fegetexceptflag failed\n");
- if(fesetexceptflag(flagp,FE_INVALID | FE_DIVBYZERO))
- mexPrintf("fesetexceptflag failed\n");
- mxFree(flagp);
- feclearexcept (FE_ALL_EXCEPT);*/
+ bool block = false;
+ double *params = NULL;
+ double *yd = NULL, *xd = NULL;
+ int count_array_argument = 0;
for (i = 0; i < nrhs; i++)
{
- if (Get_Argument(prhs[i]) == "static")
+ if (!mxIsChar(prhs[i]))
+ {
+ switch (count_array_argument)
+ {
+ case 0:
+ yd = mxGetPr(prhs[i]);
+ row_y = mxGetM(prhs[i]);
+ col_y = mxGetN(prhs[i]);
+ break;
+ case 1:
+ xd = mxGetPr(prhs[i]);
+ row_x = mxGetM(prhs[i]);
+ col_x = mxGetN(prhs[i]);
+ break;
+ case 2:
+ params = mxGetPr(prhs[i]);
+ break;
+ case 3:
+ periods = mxGetScalar(prhs[i]);
+ break;
+ case 4:
+ block_structur = mxDuplicateArray(prhs[i]);
+ break;
+ default:
+ mexPrintf("Unknown argument\n");
+ }
+ count_array_argument++;
+ }
+ else if (Get_Argument(prhs[i]) == "static")
steady_state = true;
else if (Get_Argument(prhs[i]) == "dynamic")
steady_state = false;
else if (Get_Argument(prhs[i]) == "evaluate")
evaluate = true;
+ else if (Get_Argument(prhs[i]) == "block")
+ block = true;
else
{
mexPrintf("Unknown argument : ");
@@ -267,6 +294,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexErrMsgTxt(f.c_str());
}
}
+ if (count_array_argument > 0 && count_array_argument < 4)
+ {
+ mexPrintf("Missing arguments. All the following arguments have to be indicated y, x, params, it_\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+
M_ = mexGetVariable("global", "M_");
if (M_ == NULL)
{
@@ -286,22 +319,30 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexPrintf("Global variable not found : ");
mexErrMsgTxt("options_ \n");
}
- double *params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
- double *yd, *xd, *steady_yd = NULL, *steady_xd = NULL;
+ if (!count_array_argument)
+ params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
+ double *steady_yd = NULL, *steady_xd = NULL;
if (!steady_state)
{
- yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
- row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
- col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));;
- xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
- row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
- col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
+ if (!count_array_argument)
+ {
+ yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
+ row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
+ col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));;
+ }
+ if (!count_array_argument)
+ {
+ xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
+ row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
+ col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
+ }
nb_row_xd = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_det_nbr"))))));
y_kmin = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_lag"))))));
y_kmax = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_lead"))))));
y_decal = max(0, y_kmin-int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_endo_lag")))))));
- periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "periods"))))));
+ if (!count_array_argument)
+ periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "periods"))))));
steady_yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
steady_row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
@@ -313,18 +354,26 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
}
else
{
- yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
- row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
- col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
- xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
- row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
- col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
+ if (!count_array_argument)
+ {
+ yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
+ row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
+ col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
+ }
+ if (!count_array_argument)
+ {
+ xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
+ row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
+ col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
+ }
nb_row_xd = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_det_nbr"))))));
}
int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "maxit_"))))));
double slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "slowc")))));
double markowitz_c = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "markowitz")))));
int minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "minimal_solving_periods")))));
+ int stack_solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "stack_solve_algo")))));
+ int solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_algo")))));
double solve_tolf;
if (steady_state)
solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_tolf"))));
@@ -356,9 +405,10 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int nb_row_x = row_x;
clock_t t0 = clock();
- Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods);
+ Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo);
string f(fname);
- bool result = interprete.compute_blocks(f, f, steady_state, evaluate);
+ int nb_blocks = 0;
+ bool result = interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks);
clock_t t1 = clock();
if (!steady_state && !evaluate)
mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
@@ -380,6 +430,62 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
pind[0] = 0;
else
pind[0] = 1;
+ if (nlhs > 2)
+ {
+ int jacob_field_number = 0, jacob_exo_field_number = 0, jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0;
+ if (!block_structur)
+ {
+ const char *field_names[] = {"jacob","jacob_exo","jacob_exo_det","jacob_other_endo"};
+ jacob_field_number=0;
+ jacob_exo_field_number=1;
+ jacob_exo_det_field_number=2;
+ jacob_other_endo_field_number=2;
+ mwSize dims[1] = {nb_blocks };
+ plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
+ mexPrintf("the structure has been created\n");
+ }
+ else if (!mxIsStruct(block_structur))
+ {
+ mexPrintf("The third output argument must be a structure\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ else
+ {
+ mexPrintf("Adding Fields\n");
+ jacob_field_number = mxAddField(block_structur, "jacob");
+ if (jacob_field_number == -1)
+ {
+ mexPrintf("Cannot add extra field to the structArray\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ jacob_exo_field_number = mxAddField(block_structur, "jacob_exo");
+ if (jacob_exo_field_number == -1)
+ {
+ mexPrintf("Cannot add extra field to the structArray\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ jacob_exo_det_field_number = mxAddField(block_structur, "jacob_exo_det");
+ if (jacob_exo_det_field_number == -1)
+ {
+ mexPrintf("Cannot add extra field to the structArray\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ jacob_other_endo_field_number = mxAddField(block_structur, "jacob_other_endo");
+ if (jacob_other_endo_field_number == -1)
+ {
+ mexPrintf("Cannot add extra field to the structArray\n");
+ mexErrMsgTxt("end of bytecode\n");
+ }
+ }
+ for (int i = 0; i < nb_blocks; i++)
+ {
+ mxSetFieldByNumber(block_structur,i,jacob_field_number,interprete.get_jacob(i));
+ mxSetFieldByNumber(block_structur,i,jacob_exo_field_number,interprete.get_jacob_exo(i));
+ mxSetFieldByNumber(block_structur,i,jacob_exo_det_field_number,interprete.get_jacob_exo_det(i));
+ mxSetFieldByNumber(block_structur,i,jacob_other_endo_field_number,interprete.get_jacob_other_endo(i));
+ }
+ plhs[2] = block_structur;
+ }
}
}
if (x)
diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh
index 6939ade9a..65d347cfa 100644
--- a/preprocessor/CodeInterpreter.hh
+++ b/preprocessor/CodeInterpreter.hh
@@ -57,47 +57,53 @@ using namespace std;
*/
enum Tags
{
- FLDZ, //!< Stores zero in the stack - 0
- FLDC, //!< Stores a constant term in the stack - 1
+ FLDZ, //!< Stores zero in the stack - 0 (0)
+ FLDC, //!< Stores a constant term in the stack - 1 (1)
- FDIMT, //!< Defines the number of temporary terms - dynamic context (the period has to be indicated) - 2
- FDIMST, //!< Defines the number of temporary terms - static context (the period hasn't to be indicated) - 3
- FLDT, //!< Stores a temporary term in the stack - dynamic context (the period has to be indicated) - 4
- FLDST, //!< Stores a temporary term in the stack - static context (the period hasn't to be indicated) - 5
- FSTPT, //!< Loads a temporary term from the stack - dynamic context (the period has to be indicated) - 6
- FSTPST, //!< Loads a temporary term from the stack - static context (the period hasn't to be indicated) - 7
+ FDIMT, //!< Defines the number of temporary terms - dynamic context (the period has to be indicated) - 2 (2)
+ FDIMST, //!< Defines the number of temporary terms - static context (the period hasn't to be indicated) - 3 (3)
+ FLDT, //!< Stores a temporary term in the stack - dynamic context (the period has to be indicated) - 4 (4)
+ FLDST, //!< Stores a temporary term in the stack - static context (the period hasn't to be indicated) - 5 (5)
+ FSTPT, //!< Loads a temporary term from the stack - dynamic context (the period has to be indicated) - 6 (6)
+ FSTPST, //!< Loads a temporary term from the stack - static context (the period hasn't to be indicated) - 7 (7)
- FLDU, //!< Stores an element of the vector U in the stack - dynamic context (the period has to be indicated) - 8
- FLDSU, //!< Stores an element of the vector U in the stack - static context (the period hasn't to be indicated) - 9
- FSTPU, //!< Loads an element of the vector U from the stack - dynamic context (the period has to be indicated) - A
- FSTPSU, //!< Loads an element of the vector U from the stack - static context (the period hasn't to be indicated) - B
+ FLDU, //!< Stores an element of the vector U in the stack - dynamic context (the period has to be indicated) - 8 (8)
+ FLDSU, //!< Stores an element of the vector U in the stack - static context (the period hasn't to be indicated) - 9 (9)
+ FSTPU, //!< Loads an element of the vector U from the stack - dynamic context (the period has to be indicated) - A (10)
+ FSTPSU, //!< Loads an element of the vector U from the stack - static context (the period hasn't to be indicated) - B (11)
- FLDV, //!< Stores a variable (described in SymbolType) in the stack - dynamic context (the period has to be indicated) - C
- FLDSV, //!< Stores a variable (described in SymbolType) in the stack - static context (the period hasn't to be indicated) - D
- FLDVS, //!< Stores a variable (described in SymbolType) in the stack - dynamic context but inside the STEADYSTATE function (the period hasn't to be indicated) - E
- FSTPV, //!< Loads a variable (described in SymbolType) from the stack - dynamic context (the period has to be indicated) - F
- FSTPSV, //!< Loads a variable (described in SymbolType) from the stack - static context (the period hasn't to be indicated) - 10
+ FLDV, //!< Stores a variable (described in SymbolType) in the stack - dynamic context (the period has to be indicated) - C (12)
+ FLDSV, //!< Stores a variable (described in SymbolType) in the stack - static context (the period hasn't to be indicated) - D (13)
+ FLDVS, //!< Stores a variable (described in SymbolType) in the stack - dynamic context but inside the STEADYSTATE function (the period hasn't to be indicated) - E (14)
+ FSTPV, //!< Loads a variable (described in SymbolType) from the stack - dynamic context (the period has to be indicated) - F (15)
+ FSTPSV, //!< Loads a variable (described in SymbolType) from the stack - static context (the period hasn't to be indicated) - 10 (16)
- FLDR, //!< Stores a residual in the stack - 11
- FSTPR, //!< Loads a residual from the stack - 12
+ FLDR, //!< Stores a residual in the stack - 11 (17)
+ FSTPR, //!< Loads a residual from the stack - 12 (18)
- FSTPG, //!< Loads a derivative from the stack - 13
- FSTPG2, //!< Loads a derivative matrix from the stack - 13
+ FSTPG, //!< Loads a derivative from the stack - 13 (19)
+ FSTPG2, //!< Loads a derivative matrix for static model from the stack - 14 (20)
+ FSTPG3, //!< Loads a derivative matrix for a dynamic model from the stack - 15 (21)
+ FSTPG4, //!< Loads a second order derivative matrix for a dynamic model from the stack - 16 (22)
- FUNARY, //!< A Unary operator - 14
- FBINARY, //!< A binary operator - 15
- FTRINARY, //!< A trinary operator - 15'
- FCUML, //!< Cumulates the result - 16
+ FUNARY, //!< A Unary operator - 17 (23)
+ FBINARY, //!< A binary operator - 18 (24)
+ FTRINARY, //!< A trinary operator - 19 (25)
- FBEGINBLOCK, //!< Defines the begining of a model block - 17
- FENDBLOCK, //!< Defines the end of a model block - 18
- FENDEQU, //!< Defines the last equation of the block. For block that has to be solved, the derivatives appear just after this flag - 19
- FEND, //!< Defines the end of the model code - 1A
+ FCUML, //!< Cumulates the result - 1A (26)
- FOK, //!< Used for debugging purpose - 1B
+ FJMPIFEVAL, //!< Jump if evaluate = true - 1B (27)
+ FJMP, //!< Jump - 1C (28)
- FNUMEXPR //!< Store the expression type and references
+ FBEGINBLOCK, //!< Defines the begining of a model block - 1D (29)
+ FENDBLOCK, //!< Defines the end of a model block - 1E (30)
+ FENDEQU, //!< Defines the last equation of the block. For block that has to be solved, the derivatives appear just after this flag - 1F (31)
+ FEND, //!< Defines the end of the model code - 20 (32)
+
+ FOK, //!< Used for debugging purpose - 21 (33)
+
+ FNUMEXPR //!< Store the expression type and references - 22 (34)
};
@@ -148,6 +154,7 @@ enum ExpressionType
TemporaryTerm,
ModelEquation,
FirstEndoDerivative,
+ FirstOtherEndoDerivative,
FirstExoDerivative,
FirstExodetDerivative,
FirstParamDerivative,
@@ -224,9 +231,10 @@ public:
{
};
inline void
- write(ostream &CompileCode)
+ write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast(this), sizeof(*this));
+ instruction_number++;
};
};
@@ -244,9 +252,10 @@ public:
{
};
inline void
- write(ostream &CompileCode)
+ write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast(this), sizeof(TagWithOneArgument));
+ instruction_number++;
};
};
@@ -265,9 +274,10 @@ public:
{
};
inline void
- write(ostream &CompileCode)
+ write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast(this), sizeof(*this));
+ instruction_number++;
};
};
@@ -287,12 +297,40 @@ public:
{
};
inline void
- write(ostream &CompileCode)
+ write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast(this), sizeof(*this));
+ instruction_number++;
};
};
+
+template < class T1, class T2, class T3, class T4 >
+class TagWithFourArguments
+{
+protected:
+ uint8_t op_code;
+ T1 arg1;
+ T2 arg2;
+ T3 arg3;
+ T4 arg4;
+public:
+ inline TagWithFourArguments(uint8_t op_code_arg) : op_code(op_code_arg)
+ {
+ };
+ inline TagWithFourArguments(uint8_t op_code_arg, T1 arg_arg1, T2 arg_arg2, T3 arg_arg3, T4 arg_arg4) : op_code(op_code_arg), arg1(arg_arg1), arg2(arg_arg2), arg3(arg_arg3), arg4(arg_arg4)
+ {
+ };
+ inline void
+ write(ostream &CompileCode, unsigned int &instruction_number)
+ {
+ CompileCode.write(reinterpret_cast(this), sizeof(*this));
+ instruction_number++;
+ };
+};
+
+
+
class FLDZ_ : public TagWithoutArgument
{
public:
@@ -578,6 +616,36 @@ public:
};
};
+class FSTPG3_ : public TagWithFourArguments
+{
+public:
+ inline FSTPG3_() : TagWithFourArguments::TagWithFourArguments(FSTPG3, 0, 0, 0, 0)
+ {
+ };
+ inline FSTPG3_(const unsigned int pos_arg1, const unsigned int pos_arg2, const int pos_arg3, const unsigned int pos_arg4) : TagWithFourArguments::TagWithFourArguments(FSTPG3, pos_arg1, pos_arg2, pos_arg3, pos_arg4)
+ {
+ };
+ inline unsigned int
+ get_row()
+ {
+ return arg1;
+ };
+ inline unsigned int
+ get_col()
+ {
+ return arg2;
+ };
+ inline int
+ get_lag()
+ {
+ return arg2;
+ };
+ inline unsigned int
+ get_col_pos()
+ {
+ return arg4;
+ };
+};
class FUNARY_ : public TagWithOneArgument
{
@@ -643,6 +711,38 @@ public:
};
};
+class FJMPIFEVAL_ : public TagWithOneArgument
+{
+public:
+ inline FJMPIFEVAL_() : TagWithOneArgument::TagWithOneArgument(FJMPIFEVAL)
+ {
+ };
+ inline FJMPIFEVAL_(unsigned int arg_pos) : TagWithOneArgument::TagWithOneArgument(FJMPIFEVAL, arg_pos)
+ {
+ };
+ inline unsigned int
+ get_pos()
+ {
+ return arg1;
+ }
+};
+
+class FJMP_ : public TagWithOneArgument
+{
+public:
+ inline FJMP_() : TagWithOneArgument::TagWithOneArgument(FJMP)
+ {
+ };
+ inline FJMP_(unsigned int arg_pos) : TagWithOneArgument::TagWithOneArgument(FJMP, arg_pos)
+ {
+ };
+ inline unsigned int
+ get_pos()
+ {
+ return arg1;
+ }
+};
+
class FLDVS_ : public TagWithTwoArguments
{
public:
@@ -873,9 +973,10 @@ public:
return lag3;
};
inline void
- write(ostream &CompileCode)
+ write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast(this), sizeof(FNUMEXPR_));
+ instruction_number++;
};
};
@@ -887,27 +988,52 @@ private:
uint8_t type;
vector variable;
vector equation;
+ vector other_endogenous;
+ vector exogenous;
+ vector det_exogenous;
bool is_linear;
vector Block_Contain_;
int endo_nbr;
int Max_Lag;
int Max_Lead;
int u_count_int;
+ int nb_col_jacob;
+ unsigned int det_exo_size, exo_size, other_endo_size;
+ unsigned int nb_col_other_endo_jacob;
public:
inline FBEGINBLOCK_()
{
op_code = FBEGINBLOCK; size = 0; type = UNKNOWN; /*variable = NULL; equation = NULL;*/
- is_linear = false; endo_nbr = 0; Max_Lag = 0; Max_Lead = 0; u_count_int = 0;
+ is_linear = false; endo_nbr = 0; Max_Lag = 0; Max_Lead = 0; u_count_int = 0; nb_col_jacob = 0;
};
inline FBEGINBLOCK_(unsigned int size_arg, BlockSimulationType type_arg, int unsigned first_element, int unsigned block_size,
const vector &variable_arg, const vector &equation_arg,
- bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg)
+ bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg, int nb_col_jacob_arg,
+ unsigned int det_exo_size_arg, unsigned int exo_size_arg, unsigned int other_endo_size_arg, unsigned int nb_col_other_endo_jacob_arg,
+ const vector &det_exogenous_arg, const vector &exogenous_arg, const vector &other_endogenous_arg)
{
op_code = FBEGINBLOCK; size = size_arg; type = type_arg;
variable = vector(variable_arg.begin()+first_element, variable_arg.begin()+(first_element+block_size));
equation = vector(equation_arg.begin()+first_element, equation_arg.begin()+(first_element+block_size));
- is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg; /*Block_Contain.clear();*/
+ det_exogenous = vector(det_exogenous_arg);
+ exogenous = vector(exogenous_arg);
+ other_endogenous = vector(other_endogenous_arg);
+ is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg;
+ nb_col_jacob = nb_col_jacob_arg; det_exo_size = det_exo_size_arg; exo_size = exo_size_arg; other_endo_size = other_endo_size_arg;
+ nb_col_other_endo_jacob = nb_col_other_endo_jacob_arg;
};
+ inline FBEGINBLOCK_(unsigned int size_arg, BlockSimulationType type_arg, int unsigned first_element, int unsigned block_size,
+ const vector &variable_arg, const vector &equation_arg,
+ bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg, int nb_col_jacob_arg)
+ {
+ op_code = FBEGINBLOCK; size = size_arg; type = type_arg;
+ variable = vector(variable_arg.begin()+first_element, variable_arg.begin()+(first_element+block_size));
+ equation = vector(equation_arg.begin()+first_element, equation_arg.begin()+(first_element+block_size));
+ is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg;
+ nb_col_jacob = nb_col_jacob_arg;
+ det_exo_size = 0; exo_size = 0; other_endo_size = 0;
+ nb_col_other_endo_jacob = 0;
+ }
inline unsigned int
get_size()
{
@@ -948,8 +1074,33 @@ public:
{
return Block_Contain_;
};
+ inline int
+ get_nb_col_jacob()
+ {
+ return nb_col_jacob;
+ };
+ inline unsigned int
+ get_exo_size()
+ {
+ return exo_size;
+ };
+ inline unsigned int
+ get_det_exo_size()
+ {
+ return det_exo_size;
+ };
+ inline unsigned int
+ get_other_endo_size()
+ {
+ return other_endo_size;
+ };
+ inline unsigned int
+ get_nb_col_other_endo_jacob()
+ {
+ return nb_col_other_endo_jacob;
+ };
inline void
- write(ostream &CompileCode)
+ write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast(&op_code), sizeof(op_code));
CompileCode.write(reinterpret_cast(&size), sizeof(size));
@@ -968,6 +1119,19 @@ public:
CompileCode.write(reinterpret_cast(&Max_Lead), sizeof(Max_Lead));
CompileCode.write(reinterpret_cast(&u_count_int), sizeof(u_count_int));
}
+ CompileCode.write(reinterpret_cast(&nb_col_jacob), sizeof(nb_col_jacob));
+ CompileCode.write(reinterpret_cast(&det_exo_size), sizeof(det_exo_size));
+ CompileCode.write(reinterpret_cast(&exo_size), sizeof(exo_size));
+ CompileCode.write(reinterpret_cast(&other_endo_size), sizeof(other_endo_size));
+ CompileCode.write(reinterpret_cast(&nb_col_other_endo_jacob), sizeof(nb_col_other_endo_jacob));
+
+ for (unsigned int i = 0; i < det_exo_size; i++)
+ CompileCode.write(reinterpret_cast(&det_exogenous[i]), sizeof(det_exogenous[0]));
+ for (unsigned int i = 0; i < exo_size; i++)
+ CompileCode.write(reinterpret_cast(&exogenous[i]), sizeof(exogenous[0]));
+ for (unsigned int i = 0; i < other_endo_size; i++)
+ CompileCode.write(reinterpret_cast(&other_endogenous[i]), sizeof(other_endogenous[0]));
+ instruction_number++;
};
#ifdef BYTE_CODE
@@ -993,6 +1157,30 @@ public:
memcpy(&Max_Lead, code, sizeof(Max_Lead)); code += sizeof(Max_Lead);
memcpy(&u_count_int, code, sizeof(u_count_int)); code += sizeof(u_count_int);
}
+ memcpy(&nb_col_jacob, code, sizeof(nb_col_jacob)); code += sizeof(nb_col_jacob);
+ memcpy(&det_exo_size, code, sizeof(det_exo_size)); code += sizeof(det_exo_size);
+ memcpy(&exo_size, code, sizeof(exo_size)); code += sizeof(exo_size);
+ memcpy(&other_endo_size, code, sizeof(other_endo_size)); code += sizeof(other_endo_size);
+ memcpy(&nb_col_other_endo_jacob, code, sizeof(nb_col_other_endo_jacob)); code += sizeof(nb_col_other_endo_jacob);
+
+ for (unsigned int i = 0; i < det_exo_size; i++)
+ {
+ unsigned int tmp_i;
+ memcpy(&tmp_i, code, sizeof(tmp_i)); code += sizeof(tmp_i);
+ det_exogenous.push_back(tmp_i);
+ }
+ for (unsigned int i = 0; i < exo_size; i++)
+ {
+ unsigned int tmp_i;
+ memcpy(&tmp_i, code, sizeof(tmp_i)); code += sizeof(tmp_i);
+ exogenous.push_back(tmp_i);
+ }
+ for (unsigned int i = 0; i < other_endo_size; i++)
+ {
+ unsigned int tmp_i;
+ memcpy(&tmp_i, code, sizeof(tmp_i)); code += sizeof(tmp_i);
+ other_endogenous.push_back(tmp_i);
+ }
return code;
};
#endif
@@ -1036,6 +1224,7 @@ public:
CompiledCode.close();
nb_blocks = 0;
bool done = false;
+ int instruction = 0;
while (!done)
{
switch (*code)
@@ -1181,6 +1370,20 @@ public:
tags_liste.push_back(make_pair(FSTPG, code));
code += sizeof(FSTPG_);
break;
+ case FSTPG2:
+# ifdef DEBUGL
+ mexPrintf("FSTPG2\n");
+# endif
+ tags_liste.push_back(make_pair(FSTPG2, code));
+ code += sizeof(FSTPG2_);
+ break;
+ case FSTPG3:
+# ifdef DEBUGL
+ mexPrintf("FSTPG3\n");
+# endif
+ tags_liste.push_back(make_pair(FSTPG3, code));
+ code += sizeof(FSTPG3_);
+ break;
case FUNARY:
# ifdef DEBUGL
mexPrintf("FUNARY\n");
@@ -1257,10 +1460,25 @@ public:
nb_blocks++;
}
break;
+ case FJMPIFEVAL:
+# ifdef DEBUGL
+ mexPrintf("FJMPIFEVAL\n");
+# endif
+ tags_liste.push_back(make_pair(FJMPIFEVAL, code));
+ code += sizeof(FJMPIFEVAL_);
+ break;
+ case FJMP:
+# ifdef DEBUGL
+ mexPrintf("FJMP\n");
+# endif
+ tags_liste.push_back(make_pair(FJMP, code));
+ code += sizeof(FJMP_);
+ break;
default:
mexPrintf("Unknown Tag value=%d code=%x\n", *code, code);
done = true;
}
+ instruction++;
}
return tags_liste;
};
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 4817e0c3d..1850f1150 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -57,28 +57,28 @@ DynamicModel::AddVariable(int symb_id, int lag)
}
void
-DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, const map_idx_t &map_idx) const
+DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const
{
first_derivatives_t::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag)));
if (it != first_derivatives.end())
- (it->second)->compile(code_file, false, temporary_terms, map_idx, true, false);
+ (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
else
{
FLDZ_ fldz;
- fldz.write(code_file);
+ fldz.write(code_file, instruction_number);
}
}
void
-DynamicModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, const map_idx_t &map_idx) const
+DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, const map_idx_t &map_idx) const
{
map >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
if (it != first_chain_rule_derivatives.end())
- (it->second)->compile(code_file, false, temporary_terms, map_idx, true, false);
+ (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
else
{
FLDZ_ fldz;
- fldz.write(code_file);
+ fldz.write(code_file, instruction_number);
}
}
@@ -140,7 +140,6 @@ DynamicModel::computeTemporaryTermsOrdered()
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
-
set temporary_terms_in_use;
temporary_terms_in_use.clear();
v_temporary_terms_inuse[block] = temporary_terms_in_use;
@@ -471,7 +470,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
goto evaluation;
feedback_variables.push_back(variable_ID);
output << " % equation " << equation_ID+1 << " variable : " << sModel
- << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
+ << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(eEndogenous, variable_ID) << endl;
output << " " << "residual(" << i+1-block_recursive << ") = (";
goto end;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
@@ -480,7 +479,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
goto evaluation;
feedback_variables.push_back(variable_ID);
output << " % equation " << equation_ID+1 << " variable : " << sModel
- << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
+ << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(eEndogenous, variable_ID) << endl;
Uf[equation_ID] << " b(" << i+1-block_recursive << "+Per_J_) = -residual(" << i+1-block_recursive << ", it_)";
output << " residual(" << i+1-block_recursive << ", it_) = (";
goto end;
@@ -526,6 +525,34 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
+ for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != derivative_exo[block].end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.first;
+ int var = it->first.second.second;
+ int eqr = getBlockInitialEquationID(block, eq);
+ expr_t id = it->second;
+ output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
+ << "(" << lag
+ << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != derivative_exo_det[block].end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.first;
+ int var = it->first.second.second;
+ int eqr = getBlockInitialEquationID(block, eq);
+ expr_t id = it->second;
+ output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
+ << "(" << lag
+ << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
int lag = it->first.first;
@@ -564,7 +591,36 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
+ for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != derivative_exo[block].end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.first;
+ int var = it->first.second.second;
+ int eqr = getBlockInitialEquationID(block, eq);
+ expr_t id = it->second;
+ output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
+ << "(" << lag
+ << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != derivative_exo_det[block].end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.first;
+ int var = it->first.second.second;
+ int eqr = getBlockInitialEquationID(block, eq);
+ expr_t id = it->second;
+
+ output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
+ << "(" << lag
+ << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
int lag = it->first.first;
@@ -695,6 +751,36 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
+ for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != derivative_exo[block].end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.first;
+ int var = it->first.second.second;
+ int eqr = getBlockInitialEquationID(block, eq);
+ expr_t id = it->second;
+
+ output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
+ << "(" << lag
+ << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != derivative_exo_det[block].end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.first;
+ int var = it->first.second.second;
+ int eqr = getBlockInitialEquationID(block, eq);
+ expr_t id = it->second;
+
+ output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
+ << "(" << lag
+ << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
int lag = it->first.first;
@@ -726,6 +812,7 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
{
ostringstream tmp_output;
ofstream code_file;
+ unsigned int instruction_number = 0;
bool file_open = false;
string main_name = file_name;
@@ -737,7 +824,6 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
exit(EXIT_FAILURE);
}
-
int count_u;
int u_count_int = 0;
BlockSimulationType simulation_type;
@@ -753,8 +839,15 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
//Temporary variables declaration
FDIMT_ fdimt(temporary_terms.size());
- fdimt.write(code_file);
+ fdimt.write(code_file, instruction_number);
+ int other_endo_size = 0;
+ vector exo, exo_det, other_endo;
+
+ for(int i = 0; i < symbol_table.exo_det_nbr(); i++)
+ exo_det.push_back(i);
+ for(int i = 0; i < symbol_table.exo_nbr(); i++)
+ exo.push_back(i);
FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
simulation_type,
0,
@@ -765,16 +858,24 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
symbol_table.endo_nbr(),
0,
0,
- u_count_int
+ u_count_int,
+ 0,
+ symbol_table.exo_det_nbr(),
+ symbol_table.exo_nbr(),
+ other_endo_size,
+ 0,
+ exo_det,
+ exo,
+ other_endo
);
- fbeginblock.write(code_file);
+ fbeginblock.write(code_file, instruction_number);
- compileTemporaryTerms(code_file, temporary_terms, map_idx, true, false);
+ compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
- compileModelEquations(code_file, temporary_terms, map_idx, true, false);
+ compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, true, false);
FENDEQU_ fendequ;
- fendequ.write(code_file);
+ fendequ.write(code_file, instruction_number);
vector, int > > > derivatives;
derivatives.resize(symbol_table.endo_nbr());
count_u = symbol_table.endo_nbr();
@@ -790,45 +891,45 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
unsigned int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var, lag);
- fnumexpr.write(code_file);
+ fnumexpr.write(code_file, instruction_number);
if (!derivatives[eq].size())
derivatives[eq].clear();
derivatives[eq].push_back(make_pair(make_pair(var, lag), count_u));
- d1->compile(code_file, false, temporary_terms, map_idx, true, false);
+ d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
FSTPU_ fstpu(count_u);
- fstpu.write(code_file);
+ fstpu.write(code_file, instruction_number);
count_u++;
}
}
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
FLDR_ fldr(i);
- fldr.write(code_file);
+ fldr.write(code_file, instruction_number);
for(vector, int> >::const_iterator it = derivatives[i].begin();
it != derivatives[i].end(); it++)
{
FLDU_ fldu(it->second);
- fldu.write(code_file);
+ fldu.write(code_file, instruction_number);
FLDV_ fldv(eEndogenous, it->first.first, it->first.second);
- fldv.write(code_file);
+ fldv.write(code_file, instruction_number);
FBINARY_ fbinary(oTimes);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
if (it != derivatives[i].begin())
{
FBINARY_ fbinary(oPlus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
}
}
FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FSTPU_ fstpu(i);
- fstpu.write(code_file);
+ fstpu.write(code_file, instruction_number);
}
FENDBLOCK_ fendblock;
- fendblock.write(code_file);
+ fendblock.write(code_file, instruction_number);
FEND_ fend;
- fend.write(code_file);
+ fend.write(code_file, instruction_number);
code_file.close();
}
@@ -852,6 +953,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
string tmp_s;
ostringstream tmp_output;
ofstream code_file;
+ unsigned int instruction_number = 0;
expr_t lhs = NULL, rhs = NULL;
BinaryOpNode *eq_node;
Uff Uf[symbol_table.endo_nbr()];
@@ -870,7 +972,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
//Temporary variables declaration
FDIMT_ fdimt(temporary_terms.size());
- fdimt.write(code_file);
+ fdimt.write(code_file, instruction_number);
for (unsigned int block = 0; block < getNbBlocks(); block++)
{
@@ -878,7 +980,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
if (block > 0)
{
FENDBLOCK_ fendblock;
- fendblock.write(code_file);
+ fendblock.write(code_file, instruction_number);
}
int count_u;
int u_count_int = 0;
@@ -886,6 +988,8 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
unsigned int block_size = getBlockSize(block);
unsigned int block_mfs = getBlockMfs(block);
unsigned int block_recursive = block_size - block_mfs;
+ unsigned int block_exo_det_size = exo_det_block[block].size();
+ unsigned int block_other_endo_size = other_endo_block[block].size();
int block_max_lag = max_leadlag_block[block].first;
if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE
@@ -895,6 +999,45 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE);
file_open = true;
}
+ map >, expr_t> tmp_block_endo_derivative;
+ for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+ tmp_block_endo_derivative[make_pair(it->second.first, make_pair(it->first.second, it->first.first) )] = it->second.second ;
+ map >, expr_t> tmp_exo_derivative;
+ for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != (derivative_exo[block]).end(); it++)
+ tmp_exo_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first) )] = it->second ;
+ map >, expr_t> tmp_exo_det_derivative;
+ for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != (derivative_exo_det[block]).end(); it++)
+ tmp_exo_det_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first) )] = it->second;
+ map >, expr_t> tmp_other_endo_derivative;
+ for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != (derivative_other_endo[block]).end(); it++)
+ tmp_other_endo_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first) )] = it->second;
+ int prev_var = -1;
+ int prev_lag = -999999999;
+ int count_col_endo = 0;
+ for (map >, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
+ {
+ int lag = it->first.first;
+ int var = it->first.second.first;
+ if(prev_var != var || prev_lag != lag)
+ {
+ prev_var = var;
+ prev_lag = lag;
+ count_col_endo++;
+ }
+ }
+ vector exo_det;
+ for (lag_var_t::const_iterator it = exo_det_block[block].begin(); it != exo_det_block[block].end(); it++)
+ exo_det.push_back(it->first);
+ vector exo;
+ for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
+ exo.push_back(it->first);
+ vector other_endo;
+ unsigned int count_col_other_endo = 0;
+ for (lag_var_t::const_iterator it = other_endo_block[block].begin(); it != other_endo_block[block].end(); it++)
+ {
+ other_endo.push_back(it->first);
+ count_col_other_endo += it->second.size();
+ }
FBEGINBLOCK_ fbeginblock(block_mfs,
simulation_type,
getBlockFirstEquation(block),
@@ -905,9 +1048,17 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
symbol_table.endo_nbr(),
block_max_lag,
block_max_lag,
- u_count_int
+ u_count_int,
+ count_col_endo,
+ block_exo_det_size,
+ getBlockExoColSize(block),
+ block_other_endo_size,
+ count_col_other_endo,
+ exo_det,
+ exo,
+ other_endo
);
- fbeginblock.write(code_file);
+ fbeginblock.write(code_file, instruction_number);
// The equations
for (i = 0; i < (int) block_size; i++)
@@ -921,14 +1072,15 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
it != v_temporary_terms[block][i].end(); it++)
{
FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
- fnumexpr.write(code_file);
- (*it)->compile(code_file, false, tt2, map_idx, true, false);
+ fnumexpr.write(code_file, instruction_number);
+ (*it)->compile(code_file, instruction_number, false, tt2, map_idx, true, false);
FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
- fstpt.write(code_file);
+ fstpt.write(code_file, instruction_number);
// Insert current node into tt2
tt2.insert(*it);
#ifdef DEBUGC
cout << "FSTPT " << v << "\n";
+ instruction_number++;
code_file.write(&FOK, sizeof(FOK));
code_file.write(reinterpret_cast(&k), sizeof(k));
ki++;
@@ -956,23 +1108,23 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
equ_type = getBlockEquationType(block, i);
{
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
- fnumexpr.write(code_file);
+ fnumexpr.write(code_file, instruction_number);
}
if (equ_type == E_EVALUATE)
{
eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
- lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+ lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, true, false);
}
else if (equ_type == E_EVALUATE_S)
{
eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
- lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+ lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, true, false);
}
break;
case SOLVE_BACKWARD_COMPLETE:
@@ -989,22 +1141,28 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
default:
end:
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
- fnumexpr.write(code_file);
+ fnumexpr.write(code_file, instruction_number);
eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- lhs->compile(code_file, false, temporary_terms, map_idx, true, false);
- rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
+ lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FSTPR_ fstpr(i - block_recursive);
- fstpr.write(code_file);
+ fstpr.write(code_file, instruction_number);
}
}
FENDEQU_ fendequ;
- fendequ.write(code_file);
- // The Jacobian if we have to solve the block
+ fendequ.write(code_file, instruction_number);
+
+ // Get the current code_file position and jump if eval = true
+ streampos pos1 = code_file.tellp();
+ FJMPIFEVAL_ fjmp_if_eval(0);
+ fjmp_if_eval.write(code_file, instruction_number);
+ int prev_instruction_number = instruction_number;
+ // The Jacobian if we have to solve the block determinsitic bloc
if (simulation_type != EVALUATE_BACKWARD
&& simulation_type != EVALUATE_FORWARD)
{
@@ -1013,13 +1171,13 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
{
- FNUMEXPR_ fnumexpr(FirstEndoDerivative, getBlockEquationID(block, i), getBlockVariableID(block, 0), 0);
- fnumexpr.write(code_file);
+ FNUMEXPR_ fnumexpr(FirstEndoDerivative, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0);
+ fnumexpr.write(code_file, instruction_number);
}
- compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
+ compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
{
FSTPG_ fstpg(0);
- fstpg.write(code_file);
+ fstpg.write(code_file, instruction_number);
}
break;
@@ -1030,12 +1188,12 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
count_u = feedback_variables.size();
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
{
+ int lag = it->second.first;
unsigned int eq = it->first.first;
unsigned int var = it->first.second;
unsigned int eqr = getBlockEquationID(block, eq);
unsigned int varr = getBlockVariableID(block, var);
- int lag = it->second.first;
- if (eq >= block_recursive && var >= block_recursive)
+ if (eq >= block_recursive and var >= block_recursive)
{
if (!Uf[eqr].Ufl)
{
@@ -1052,10 +1210,10 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
Uf[eqr].Ufl->var = varr;
Uf[eqr].Ufl->lag = lag;
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, lag);
- fnumexpr.write(code_file);
- compileChainRuleDerivative(code_file, eqr, varr, lag, map_idx);
+ fnumexpr.write(code_file, instruction_number);
+ compileChainRuleDerivative(code_file, instruction_number, eqr, varr, lag, map_idx);
FSTPU_ fstpu(count_u);
- fstpu.write(code_file);
+ fstpu.write(code_file, instruction_number);
count_u++;
}
}
@@ -1064,24 +1222,24 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
if (i >= (int) block_recursive)
{
FLDR_ fldr(i-block_recursive);
- fldr.write(code_file);
+ fldr.write(code_file, instruction_number);
FLDZ_ fldz;
- fldz.write(code_file);
+ fldz.write(code_file, instruction_number);
v = getBlockEquationID(block, i);
for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
{
FLDU_ fldu(Uf[v].Ufl->u);
- fldu.write(code_file);
+ fldu.write(code_file, instruction_number);
FLDV_ fldv(eEndogenous, Uf[v].Ufl->var, Uf[v].Ufl->lag);
- fldv.write(code_file);
+ fldv.write(code_file, instruction_number);
FBINARY_ fbinary(oTimes);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FCUML_ fcuml;
- fcuml.write(code_file);
+ fcuml.write(code_file, instruction_number);
}
Uf[v].Ufl = Uf[v].Ufl_First;
while (Uf[v].Ufl)
@@ -1091,10 +1249,10 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
Uf[v].Ufl = Uf[v].Ufl_First;
}
FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FSTPU_ fstpu(i - block_recursive);
- fstpu.write(code_file);
+ fstpu.write(code_file, instruction_number);
}
}
break;
@@ -1102,11 +1260,125 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
break;
}
}
+ // Get the current code_file position and jump = true
+ streampos pos2 = code_file.tellp();
+ FJMP_ fjmp(0);
+ fjmp.write(code_file, instruction_number);
+ // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
+ streampos pos3 = code_file.tellp();
+ code_file.seekp(pos1);
+ FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
+ fjmp_if_eval1.write(code_file, instruction_number);
+ code_file.seekp(pos3);
+ prev_instruction_number = instruction_number ;
+ // The Jacobian if we have to solve the block determinsitic bloc
+
+ prev_var = -1;
+ prev_lag = -999999999;
+ count_col_endo = 0;
+ for (map >, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
+ {
+ int lag = it->first.first;
+ unsigned int eq = it->first.second.second;
+ int var = it->first.second.first;
+ unsigned int eqr = getBlockEquationID(block, eq);
+ unsigned int varr = getBlockVariableID(block, var);
+ if(prev_var != var || prev_lag != lag)
+ {
+ prev_var = var;
+ prev_lag = lag;
+ count_col_endo++;
+ }
+ FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, lag);
+ fnumexpr.write(code_file, instruction_number);
+ compileDerivative(code_file, instruction_number, eqr, varr, lag, map_idx);
+ FSTPG3_ fstpg3(eq, var, lag, count_col_endo-1);
+ fstpg3.write(code_file, instruction_number);
+ }
+ prev_var = -1;
+ prev_lag = -999999999;
+ int count_col_exo = 0;
+ for (map >, expr_t>::const_iterator it = tmp_exo_derivative.begin(); it != tmp_exo_derivative.end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.second;
+ int var = it->first.second.first;
+ int eqr = getBlockInitialEquationID(block, eq);
+ int varr = getBlockInitialExogenousID(block, var);
+ if(prev_var != var || prev_lag != lag)
+ {
+ prev_var = var;
+ prev_lag = lag;
+ count_col_exo++;
+ }
+ expr_t id = it->second;
+
+ FNUMEXPR_ fnumexpr(FirstExoDerivative, eqr, varr, lag);
+ fnumexpr.write(code_file, instruction_number);
+ id->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+ FSTPG3_ fstpg3(eq, var, lag, /*var*/count_col_exo-1);
+ fstpg3.write(code_file, instruction_number);
+ }
+ prev_var = -1;
+ prev_lag = -999999999;
+ int count_col_exo_det = 0;
+ for (map >, expr_t>::const_iterator it = tmp_exo_det_derivative.begin(); it != tmp_exo_det_derivative.end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.second;
+ int var = it->first.second.first;
+ int eqr = getBlockInitialEquationID(block, eq);
+ int varr = getBlockInitialDetExogenousID(block, var);
+ if(prev_var != var || prev_lag != lag)
+ {
+ prev_var = var;
+ prev_lag = lag;
+ count_col_exo_det++;
+ }
+ expr_t id = it->second;
+
+ FNUMEXPR_ fnumexpr(FirstExodetDerivative, eqr, varr, lag);
+ fnumexpr.write(code_file, instruction_number);
+ id->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+ FSTPG3_ fstpg3(eq, var, lag, count_col_exo_det-1);
+ fstpg3.write(code_file, instruction_number);
+ }
+ prev_var = -1;
+ prev_lag = -999999999;
+ count_col_other_endo = 0;
+ for (map >, expr_t>::const_iterator it = tmp_other_endo_derivative.begin(); it != tmp_other_endo_derivative.end(); it++)
+ {
+ int lag = it->first.first;
+ int eq = it->first.second.second;
+ int var = it->first.second.first;
+ int eqr = getBlockInitialEquationID(block, eq);
+ int varr = getBlockInitialOtherEndogenousID(block, var);;
+ if(prev_var != var || prev_lag != lag)
+ {
+ prev_var = var;
+ prev_lag = lag;
+ count_col_other_endo++;
+ }
+ expr_t id = it->second;
+
+ FNUMEXPR_ fnumexpr(FirstOtherEndoDerivative, eqr, varr, lag);
+ fnumexpr.write(code_file, instruction_number);
+ id->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+ FSTPG3_ fstpg3(eq, var, lag, count_col_other_endo-1);
+ fstpg3.write(code_file, instruction_number);
+ }
+
+ // Set codefile position to previous JMP_ and set the number of instructions to jump
+ pos1 = code_file.tellp();
+ code_file.seekp(pos2);
+ FJMP_ fjmp1(instruction_number - prev_instruction_number);
+ fjmp1.write(code_file, instruction_number);
+ code_file.seekp(pos1);
}
FENDBLOCK_ fendblock;
- fendblock.write(code_file);
+ fendblock.write(code_file, instruction_number);
FEND_ fend;
- fend.write(code_file);
+ fend.write(code_file, instruction_number);
code_file.close();
}
@@ -1909,7 +2181,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
if (block_decomposition)
{
int count_lead_lag_incidence = 0;
- int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo;
+ int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo, max_lag_exo_det, max_lead_exo_det;
unsigned int nb_blocks = getNbBlocks();
for (unsigned int block = 0; block < nb_blocks; block++)
{
@@ -1921,11 +2193,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
max_lead = max_leadlag_block[block].second;
max_lag_endo = endo_max_leadlag_block[block].first;
max_lead_endo = endo_max_leadlag_block[block].second;
- max_lag_exo = max(exo_max_leadlag_block[block].first, exo_det_max_leadlag_block[block].first);
- max_lead_exo = max(exo_max_leadlag_block[block].second, exo_det_max_leadlag_block[block].second);
- vector exogenous(symbol_table.exo_nbr(), -1);
- vector::iterator it_exogenous;
- exogenous.clear();
+ max_lag_exo = exo_max_leadlag_block[block].first;
+ max_lead_exo = exo_max_leadlag_block[block].second;
+ max_lag_exo_det = exo_det_max_leadlag_block[block].first;
+ max_lead_exo_det = exo_det_max_leadlag_block[block].second;
ostringstream tmp_s, tmp_s_eq;
tmp_s.str("");
tmp_s_eq.str("");
@@ -1934,10 +2205,21 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
tmp_s << " " << getBlockVariableID(block, i)+1;
tmp_s_eq << " " << getBlockEquationID(block, i)+1;
}
- it_exogenous = exogenous.begin();
+ set exogenous;
+ exogenous.clear();
for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
- it_exogenous = set_union(it->second.begin(), it->second.end(), exogenous.begin(), exogenous.end(), it_exogenous);
- output << "M_.block_structure.block(" << block+1 << ").num = " << block+1 << ";\n";
+ for(var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+ exogenous.insert(*it1);
+ set exogenous_det;
+ exogenous_det.clear();
+ for (lag_var_t::const_iterator it = exo_det_block[block].begin(); it != exo_det_block[block].end(); it++)
+ for(var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+ exogenous_det.insert(*it1);
+ set other_endogenous;
+ other_endogenous.clear();
+ for (lag_var_t::const_iterator it = other_endo_block[block].begin(); it != other_endo_block[block].end(); it++)
+ for(var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+ other_endogenous.insert(*it1);
output << "M_.block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_lead = " << max_lead << ";\n";
@@ -1945,21 +2227,43 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
output << "M_.block_structure.block(" << block+1 << ").maximum_endo_lead = " << max_lead_endo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_lag = " << max_lag_exo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_lead = " << max_lead_exo << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").maximum_exo_det_lag = " << max_lag_exo_det << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").maximum_exo_det_lead = " << max_lead_exo_det << ";\n";
output << "M_.block_structure.block(" << block+1 << ").endo_nbr = " << block_size << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").mfs = " << getBlockMfs(block) << ";\n";
output << "M_.block_structure.block(" << block+1 << ").equation = [" << tmp_s_eq.str() << "];\n";
output << "M_.block_structure.block(" << block+1 << ").variable = [" << tmp_s.str() << "];\n";
+ output << "M_.block_structure.block(" << block+1 << ").exo_nbr = " << getBlockExoSize(block) << ";\n";
output << "M_.block_structure.block(" << block+1 << ").exogenous = [";
int i = 0;
- for (it_exogenous = exogenous.begin(); it_exogenous != exogenous.end(); it_exogenous++)
+ for (set::iterator it_exogenous = exogenous.begin(); it_exogenous != exogenous.end(); it_exogenous++)
if (*it_exogenous >= 0)
{
output << " " << *it_exogenous+1;
i++;
}
output << "];\n";
- output << "M_.block_structure.block(" << block+1 << ").exo_nbr = " << i << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").exo_det_nbr = " << i << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").exogenous_det = [";
+ i = 0;
+ for (set::iterator it_exogenous_det = exogenous_det.begin(); it_exogenous_det != exogenous_det.end(); it_exogenous_det++)
+ if (*it_exogenous_det >= 0)
+ {
+ output << " " << *it_exogenous_det+1;
+ i++;
+ }
+ output << "];\n";
- output << "M_.block_structure.block(" << block+1 << ").exo_det_nbr = " << exo_det_block.size() << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").other_endo_nbr = " << i << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").other_endogenous = [";
+ i = 0;
+ for (set::iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++)
+ if (*it_other_endogenous >= 0)
+ {
+ output << " " << *it_other_endogenous+1;
+ i++;
+ }
+ output << "];\n";
tmp_s.str("");
count_lead_lag_incidence = 0;
@@ -1970,7 +2274,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
int last_var = -1;
for (int lag = -max_lag_endo; lag < max_lead_endo+1; lag++)
{
- last_var = 0;
+ last_var = -1;
for (dynamic_jacob_map_t::const_iterator it = reordered_dynamic_jacobian.begin(); it != reordered_dynamic_jacobian.end(); it++)
{
if (lag == it->first.first && last_var != it->first.second.first)
@@ -1989,6 +2293,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
output << "M_.block_structure.block(" << block+1 << ").lead_lag_incidence = [ M_.block_structure.block(" << block+1 << ").lead_lag_incidence; " << tmp_s.str() << "]; %lag = " << lag << "\n";
tmp_s.str("");
}
+ output << "M_.block_structure.block(" << block+1 << ").n_static = " << block_col_type[block].first.first << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").n_forward = " << block_col_type[block].first.second << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").n_backward = " << block_col_type[block].second.first << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").n_mixed = " << block_col_type[block].second.second << ";\n";
}
string cst_s;
int nb_endo = symbol_table.endo_nbr();
@@ -2027,7 +2335,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
output << "M_.block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").lead_lag = " << prev_lag << ";\n";
output << "M_.block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").sparse_IM = [";
}
- output << it->first.second.first << " " << it->first.second.second << ";\n";
+ output << it->first.second.first+1 << " " << it->first.second.second+1 << ";\n";
}
output << "];\n";
}
@@ -2142,6 +2450,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
if (block)
{
+ vector n_static, n_forward, n_backward, n_mixed;
jacob_map_t contemporaneous_jacobian, static_jacobian;
// for each block contains pair
@@ -2159,10 +2468,11 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
cout << "Finding the optimal block decomposition of the model ...\n";
- if (prologue+epilogue < (unsigned int) equation_number())
- computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered);
+ lag_lead_vector_t equation_lag_lead, variable_lag_lead;
- block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered);
+ computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed);
+
+ block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type);
printBlockDecomposition(blocks);
@@ -2172,10 +2482,24 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
collect_block_first_order_derivatives();
+ collectBlockVariables();
+
global_temporary_terms = true;
if (!no_tmp_terms)
computeTemporaryTermsOrdered();
-
+ int k = 0;
+ equation_block = vector(equation_number());
+ variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equation_number());
+ for (unsigned int i = 0; i < getNbBlocks(); i++)
+ {
+ for(unsigned int j = 0; j< getBlockSize(i); j++)
+ {
+ equation_block[equation_reordered[k]] = i;
+ int l = variable_reordered[k];
+ variable_block_lead_lag[l] = make_pair(i, make_pair(variable_lag_lead[l].first, variable_lag_lead[l].second));
+ k++;
+ }
+ }
}
else
if (!no_tmp_terms)
@@ -2241,82 +2565,52 @@ DynamicModel::get_Derivatives(int block)
}
void
-DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
+DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_endo_derivatives)
{
map recursive_variables;
unsigned int nb_blocks = getNbBlocks();
- blocks_derivatives = blocks_derivatives_t(nb_blocks);
+ blocks_endo_derivatives = blocks_derivatives_t(nb_blocks);
for (unsigned int block = 0; block < nb_blocks; block++)
{
block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
recursive_variables.clear();
- BlockSimulationType simulation_type = getBlockSimulationType(block);
int block_size = getBlockSize(block);
int block_nb_mfs = getBlockMfs(block);
int block_nb_recursives = block_size - block_nb_mfs;
- if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
+ blocks_endo_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
+ for (int i = 0; i < block_nb_recursives; i++)
{
- blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
- for (int i = 0; i < block_nb_recursives; i++)
+ if (getBlockEquationType(block, i) == E_EVALUATE_S)
+ recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
+ else
+ recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
+ }
+ map >, pair >, int> Derivatives = get_Derivatives(block);
+ map >, pair >, int>::const_iterator it = Derivatives.begin();
+ for (int i = 0; i < (int) Derivatives.size(); i++)
+ {
+ int Deriv_type = it->second;
+ pair >, pair > it_l(it->first);
+ it++;
+ int lag = it_l.first.first;
+ int eq = it_l.first.second.first;
+ int var = it_l.first.second.second;
+ int eqr = it_l.second.first;
+ int varr = it_l.second.second;
+ if (Deriv_type == 0)
+ first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))];
+ else if (Deriv_type == 1)
+ first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
+ else if (Deriv_type == 2)
{
- if (getBlockEquationType(block, i) == E_EVALUATE_S)
- recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
- else
- recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
- }
- map >, pair >, int> Derivatives = get_Derivatives(block);
- map >, pair >, int>::const_iterator it = Derivatives.begin();
- for (int i = 0; i < (int) Derivatives.size(); i++)
- {
- int Deriv_type = it->second;
- pair >, pair > it_l(it->first);
- it++;
- int lag = it_l.first.first;
- int eq = it_l.first.second.first;
- int var = it_l.first.second.second;
- int eqr = it_l.second.first;
- int varr = it_l.second.second;
- if (Deriv_type == 0)
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))];
- else if (Deriv_type == 1)
+ if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursives)
first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
- else if (Deriv_type == 2)
- {
- if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursives)
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
- else
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
- }
- tmp_derivatives.push_back(make_pair(make_pair(eq, var), make_pair(lag, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))])));
- }
- }
- else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE
- || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
- {
- blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
- for (int i = 0; i < block_nb_recursives; i++)
- {
- if (getBlockEquationType(block, i) == E_EVALUATE_S)
- recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
else
- recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
- }
- for (int eq = block_nb_recursives; eq < block_size; eq++)
- {
- int eqr = getBlockEquationID(block, eq);
- for (int var = block_nb_recursives; var < block_size; var++)
- {
- int varr = getBlockVariableID(block, var);
- expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
- if (d1 == Zero)
- continue;
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1;
- tmp_derivatives.push_back(
- make_pair(make_pair(eq, var), make_pair(0, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))])));
- }
+ first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
}
+ tmp_derivatives.push_back(make_pair(make_pair(eq, var), make_pair(lag, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))])));
}
- blocks_derivatives[block] = tmp_derivatives;
+ blocks_endo_derivatives[block] = tmp_derivatives;
}
}
@@ -2379,6 +2673,20 @@ DynamicModel::collect_block_first_order_derivatives()
if (lag > 0 && lag > other_endo_max_leadlag_block[block_eq].second)
other_endo_max_leadlag_block[block_eq] = make_pair(other_endo_max_leadlag_block[block_eq].first, lag);
tmp_derivative = derivative_other_endo[block_eq];
+ {
+ map< int, map >::const_iterator it = block_other_endo_index.find(block_eq);
+ if (it == block_other_endo_index.end())
+ block_other_endo_index[block_eq][var] = 0;
+ else
+ {
+ map::const_iterator it1 = it->second.find(var);
+ if (it1 == it->second.end())
+ {
+ int size = block_other_endo_index[block_eq].size();
+ block_other_endo_index[block_eq][var] = size;
+ }
+ }
+ }
tmp_derivative[make_pair(lag, make_pair(eq, var))] = first_derivatives[make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, var), lag))];
derivative_other_endo[block_eq] = tmp_derivative;
lag_var = other_endo_block[block_eq];
@@ -2394,6 +2702,20 @@ DynamicModel::collect_block_first_order_derivatives()
if (lag > 0 && lag > exo_max_leadlag_block[block_eq].second)
exo_max_leadlag_block[block_eq] = make_pair(exo_max_leadlag_block[block_eq].first, lag);
tmp_derivative = derivative_exo[block_eq];
+ {
+ map< int, map >::const_iterator it = block_exo_index.find(block_eq);
+ if (it == block_exo_index.end())
+ block_exo_index[block_eq][var] = 0;
+ else
+ {
+ map::const_iterator it1 = it->second.find(var);
+ if (it1 == it->second.end())
+ {
+ int size = block_exo_index[block_eq].size();
+ block_exo_index[block_eq][var] = size;
+ }
+ }
+ }
tmp_derivative[make_pair(lag, make_pair(eq, var))] = first_derivatives[make_pair(eq, getDerivID(symbol_table.getID(eExogenous, var), lag))];
derivative_exo[block_eq] = tmp_derivative;
lag_var = exo_block[block_eq];
@@ -2408,6 +2730,20 @@ DynamicModel::collect_block_first_order_derivatives()
if (lag > 0 && lag > exo_det_max_leadlag_block[block_eq].second)
exo_det_max_leadlag_block[block_eq] = make_pair(exo_det_max_leadlag_block[block_eq].first, lag);
tmp_derivative = derivative_exo_det[block_eq];
+ {
+ map< int, map >::const_iterator it = block_det_exo_index.find(block_eq);
+ if (it == block_det_exo_index.end())
+ block_det_exo_index[block_eq][var] = 0;
+ else
+ {
+ map::const_iterator it1 = it->second.find(var);
+ if (it1 == it->second.end())
+ {
+ int size = block_det_exo_index[block_eq].size();
+ block_det_exo_index[block_eq][var] = size;
+ }
+ }
+ }
tmp_derivative[make_pair(lag, make_pair(eq, var))] = first_derivatives[make_pair(eq, getDerivID(symbol_table.getID(eExogenous, var), lag))];
derivative_exo_det[block_eq] = tmp_derivative;
lag_var = exo_det_block[block_eq];
@@ -2427,6 +2763,35 @@ DynamicModel::collect_block_first_order_derivatives()
}
+
+void DynamicModel::collectBlockVariables()
+{
+ for (unsigned int block = 0; block < getNbBlocks(); block++)
+ {
+ int prev_var = -1;
+ int prev_lag = -999999999;
+ int count_col_exo = 0;
+ var_t tmp_var_exo;
+ for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
+ {
+ int lag = it->first;
+ for(var_t::const_iterator it2 = it->second.begin(); it2 != it->second.end(); it2++)
+ {
+ int var = *it2;
+ tmp_var_exo.insert(var);
+ if(prev_var != var || prev_lag != lag)
+ {
+ prev_var = var;
+ prev_lag = lag;
+ count_col_exo++;
+ }
+ }
+ }
+ block_var_exo.push_back(make_pair(tmp_var_exo, count_col_exo));
+ }
+}
+
+
void
DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order) const
{
diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh
index 32e2296af..1d42f281d 100644
--- a/preprocessor/DynamicModel.hh
+++ b/preprocessor/DynamicModel.hh
@@ -142,9 +142,9 @@ private:
//! creates a mapping from the index of temporary terms to a natural index
void computeTemporaryTermsMapping();
//! Write derivative code of an equation w.r. to a variable
- void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, const map_idx_t &map_idx) const;
+ void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const;
//! Write chain rule derivative code of an equation w.r. to a variable
- void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, const map_idx_t &map_idx) const;
+ void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, const map_idx_t &map_idx) const;
//! Get the type corresponding to a derivation ID
virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
@@ -180,6 +180,9 @@ private:
//! Collecte the derivatives w.r. to endogenous of the block, to endogenous of previouys blocks and to exogenous
void collect_block_first_order_derivatives();
+ //! Collecte the informations about exogenous, deterministic exogenous and endogenous from the previous block for each block
+ void collectBlockVariables();
+
//! Factorized code for substitutions of leads/lags
/*! \param[in] type determines which type of variables is concerned */
void substituteLeadLagInternal(aux_var_t type);
@@ -211,11 +214,26 @@ private:
//! Vector of derivative for each blocks
vector derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
- //!List for each block and for each lag-leag all the other endogenous variables and exogenous variables
+ //!List for each block and for each lag-lead all the other endogenous variables and exogenous variables
typedef set var_t;
typedef map lag_var_t;
vector other_endo_block, exo_block, exo_det_block;
+ //!List for each block the exogenous variables
+ vector > block_var_exo;
+
+ map< int, map > block_exo_index, block_det_exo_index, block_other_endo_index;
+
+ //! for each block described the number of static, forward, backward and mixed variables in the block
+ /*! pair< pair, pair > */
+ vector, pair > > block_col_type;
+
+ //! List for each variable its block number and its maximum lag and lead inside the block
+ vector > > variable_block_lead_lag;
+ //! List for each equation its block number
+ vector equation_block;
+
+
//!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
vector > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
@@ -323,6 +341,18 @@ public:
{
return (block_type_firstequation_size_mfs[block_number].second.first);
};
+ //! Return the number of exogenous variable in the block block_number
+ virtual unsigned int
+ getBlockExoSize(int block_number) const
+ {
+ return (block_var_exo[block_number].first.size());
+ };
+ //! Return the number of colums in the jacobian matrix for exogenous variable in the block block_number
+ virtual unsigned int
+ getBlockExoColSize(int block_number) const
+ {
+ return (block_var_exo[block_number].second);
+ };
//! Return the number of feedback variable of the block block_number
virtual unsigned int
getBlockMfs(int block_number) const
@@ -377,6 +407,13 @@ public:
{
return (variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);
};
+ //! Return the original number of the exogenous variable varexo_number belonging to the block block_number
+ virtual int
+ getBlockVariableExoID(int block_number, int variable_number) const
+ {
+ map::const_iterator it = exo_block[block_number].find(variable_number);
+ return (it->first);
+ };
//! Return the position of equation_number in the block number belonging to the block block_number
virtual int
getBlockInitialEquationID(int block_number, int equation_number) const
@@ -389,7 +426,60 @@ public:
{
return ((int) inv_variable_reordered[variable_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
};
-
+ //! Return the block number containing the endogenous variable variable_number
+ int
+ getBlockVariableID(int variable_number) const
+ {
+ return(variable_block_lead_lag[variable_number].first);
+ };
+ //! Return the position of the exogenous variable_number in the block number belonging to the block block_number
+ virtual int
+ getBlockInitialExogenousID(int block_number, int variable_number) const
+ {
+ map< int, map >::const_iterator it = block_exo_index.find(block_number);
+ if (it != block_exo_index.end())
+ {
+ map::const_iterator it1 = it->second.find(variable_number);
+ if( it1 != it->second.end())
+ return it1->second;
+ else
+ return -1;
+ }
+ else
+ return (-1);
+ };
+ //! Return the position of the deterministic exogenous variable_number in the block number belonging to the block block_number
+ virtual int
+ getBlockInitialDetExogenousID(int block_number, int variable_number) const
+ {
+ map< int, map >::const_iterator it = block_det_exo_index.find(block_number);
+ if (it != block_det_exo_index.end())
+ {
+ map::const_iterator it1 = it->second.find(variable_number);
+ if( it1 != it->second.end())
+ return it1->second;
+ else
+ return -1;
+ }
+ else
+ return (-1);
+ };
+ //! Return the position of the other endogenous variable_number in the block number belonging to the block block_number
+ virtual int
+ getBlockInitialOtherEndogenousID(int block_number, int variable_number) const
+ {
+ map< int, map >::const_iterator it = block_other_endo_index.find(block_number);
+ if (it != block_other_endo_index.end())
+ {
+ map::const_iterator it1 = it->second.find(variable_number);
+ if( it1 != it->second.end())
+ return it1->second;
+ else
+ return -1;
+ }
+ else
+ return (-1);
+ };
};
#endif
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index 6dedb0c67..9e0ea561a 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -300,10 +300,10 @@ NumConstNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
}
void
-NumConstNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+NumConstNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
FLDC_ fldc(datatree.num_constants.getDouble(id));
- fldc.write(CompileCode);
+ fldc.write(CompileCode, instruction_number);
}
void
@@ -689,10 +689,10 @@ VariableNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
}
void
-VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
if (type == eModelLocalVariable || type == eModFileLocalVariable)
- datatree.local_variables_table[symb_id]->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+ datatree.local_variables_table[symb_id]->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
else
{
int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
@@ -705,26 +705,26 @@ VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
if (steady_dynamic) // steady state values in a dynamic model
{
FLDVS_ fldvs(type, tsid);
- fldvs.write(CompileCode);
+ fldvs.write(CompileCode, instruction_number);
}
else
{
if (type == eParameter)
{
FLDV_ fldv(type, tsid);
- fldv.write(CompileCode);
+ fldv.write(CompileCode, instruction_number);
}
else
{
FLDV_ fldv(type, tsid, lag);
- fldv.write(CompileCode);
+ fldv.write(CompileCode, instruction_number);
}
}
}
else
{
FLDSV_ fldsv(type, tsid);
- fldsv.write(CompileCode);
+ fldsv.write(CompileCode, instruction_number);
}
}
else
@@ -741,19 +741,19 @@ VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
if (type == eParameter)
{
FSTPV_ fstpv(type, tsid);
- fstpv.write(CompileCode);
+ fstpv.write(CompileCode, instruction_number);
}
else
{
FSTPV_ fstpv(type, tsid, lag);
- fstpv.write(CompileCode);
+ fstpv.write(CompileCode, instruction_number);
}
}
}
else
{
FSTPSV_ fstpsv(type, tsid);
- fstpsv.write(CompileCode);
+ fstpsv.write(CompileCode, instruction_number);
}
}
}
@@ -1567,7 +1567,7 @@ UnaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException
}
void
-UnaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+UnaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
temporary_terms_t::const_iterator it = temporary_terms.find(const_cast(this));
if (it != temporary_terms.end())
@@ -1576,23 +1576,23 @@ UnaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDT_ fldt(ii->second);
- fldt.write(CompileCode);
+ fldt.write(CompileCode, instruction_number);
}
else
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDST_ fldst(ii->second);
- fldst.write(CompileCode);
+ fldst.write(CompileCode, instruction_number);
}
return;
}
if (op_code == oSteadyState)
- arg->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, true);
+ arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, true);
else
{
- arg->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+ arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
FUNARY_ funary(op_code);
- funary.write(CompileCode);
+ funary.write(CompileCode, instruction_number);
}
}
@@ -2271,7 +2271,7 @@ BinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
}
void
-BinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
// If current node is a temporary term
temporary_terms_t::const_iterator it = temporary_terms.find(const_cast(this));
@@ -2281,20 +2281,20 @@ BinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDT_ fldt(ii->second);
- fldt.write(CompileCode);
+ fldt.write(CompileCode, instruction_number);
}
else
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDST_ fldst(ii->second);
- fldst.write(CompileCode);
+ fldst.write(CompileCode, instruction_number);
}
return;
}
- arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
- arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+ arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+ arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
FBINARY_ fbinary(op_code);
- fbinary.write(CompileCode);
+ fbinary.write(CompileCode, instruction_number);
}
void
@@ -3205,7 +3205,7 @@ TrinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalExcepti
}
void
-TrinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+TrinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
// If current node is a temporary term
temporary_terms_t::const_iterator it = temporary_terms.find(const_cast(this));
@@ -3215,21 +3215,21 @@ TrinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDT_ fldt(ii->second);
- fldt.write(CompileCode);
+ fldt.write(CompileCode, instruction_number);
}
else
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDST_ fldst(ii->second);
- fldst.write(CompileCode);
+ fldst.write(CompileCode, instruction_number);
}
return;
}
- arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
- arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
- arg3->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+ arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+ arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+ arg3->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
FTRINARY_ ftrinary(op_code);
- ftrinary.write(CompileCode);
+ ftrinary.write(CompileCode, instruction_number);
}
void
@@ -3665,7 +3665,7 @@ ExternalFunctionNode::eval(const eval_context_t &eval_context) const throw (Eval
}
void
-ExternalFunctionNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+ExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
cerr << "ExternalFunctionNode::compile: operation impossible!" << endl;
exit(EXIT_FAILURE);
diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh
index 79773a6be..825e43439 100644
--- a/preprocessor/ExprNode.hh
+++ b/preprocessor/ExprNode.hh
@@ -241,7 +241,7 @@ public:
};
virtual double eval(const eval_context_t &eval_context) const throw (EvalException) = 0;
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const = 0;
+ virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const = 0;
//! Creates a static version of this node
/*!
This method duplicates the current node by creating a similar node from which all leads/lags have been stripped,
@@ -387,7 +387,7 @@ public:
virtual void collectVariables(SymbolType type_arg, set > &result) const;
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+ virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t toStatic(DataTree &static_datatree) const;
virtual pair normalizeEquation(int symb_id_endo, vector > > &List_of_Op_RHS) const;
virtual expr_t getChainRuleDerivative(int deriv_id, const map &recursive_variables);
@@ -429,8 +429,13 @@ public:
int equation) const;
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+ virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t toStatic(DataTree &static_datatree) const;
+ SymbolType
+ get_type() const
+ {
+ return type;
+ };
int
get_symb_id() const
{
@@ -486,7 +491,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException);
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+ virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
//! Returns operand
expr_t
get_arg() const
@@ -549,7 +554,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+ virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
//! Returns first operand
expr_t
@@ -620,7 +625,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+ virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t toStatic(DataTree &static_datatree) const;
virtual pair normalizeEquation(int symb_id_endo, vector > > &List_of_Op_RHS) const;
virtual expr_t getChainRuleDerivative(int deriv_id, const map &recursive_variables);
@@ -678,7 +683,7 @@ public:
virtual void collectVariables(SymbolType type_arg, set > &result) const;
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+ virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t toStatic(DataTree &static_datatree) const;
virtual pair normalizeEquation(int symb_id_endo, vector > > &List_of_Op_RHS) const;
virtual expr_t getChainRuleDerivative(int deriv_id, const map &recursive_variables);
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 5a5a52c70..50a62e616 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -302,7 +302,7 @@ ModFile::computingPass(bool no_tmp_terms)
if (mod_file_struct.simul_present)
{
dynamic_model.initializeVariablesAndEquations();
- dynamic_model.computingPass(false, false, false, false, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
+ dynamic_model.computingPass(true, false, false, false, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
}
else
{
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 005c47937..ede21b364 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -207,6 +207,8 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = 0;
if (contemporaneous_jacobian.find(make_pair(it->first.first, it->first.second)) == contemporaneous_jacobian.end())
contemporaneous_jacobian[make_pair(it->first.first, it->first.second)] = 0;
+ if (first_derivatives.find(make_pair(it->first.first, getDerivID(symbol_table.getID(eEndogenous, it->first.second), 0))) == first_derivatives.end())
+ first_derivatives[make_pair(it->first.first, getDerivID(symbol_table.getID(eEndogenous, it->first.second), 0))] = Zero;
}
}
}
@@ -486,7 +488,7 @@ ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian
for (dynamic_jacob_map_t::const_iterator it = dynamic_jacobian.begin(); it != dynamic_jacobian.end(); it++)
{
int lag = it->first.first;
- int j_1 = it->first.second.second;
+ int j_1 = it->first.second.first;
int i_1 = it->first.second.second;
if (variable_blck[i_1] == equation_blck[j_1])
{
@@ -503,7 +505,7 @@ ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian
}
void
-ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector &equation_reordered, vector &variable_reordered, vector > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector &inv_equation_reordered, vector &inv_variable_reordered) const
+ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector &equation_reordered, vector &variable_reordered, vector > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector &inv_equation_reordered, vector &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed) const
{
int nb_var = variable_reordered.size();
int n = nb_var - prologue - epilogue;
@@ -570,8 +572,6 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
components_set[endo2block[i]].first.insert(i);
}
- lag_lead_vector_t equation_lag_lead, variable_lag_lead;
-
getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
vector tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
@@ -594,6 +594,23 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
add_edge(i, i, G2);
}
+ //Determines the dynamic structure of each equation
+ n_static = vector(prologue+num+epilogue, 0);
+ n_forward = vector(prologue+num+epilogue, 0);
+ n_backward = vector(prologue+num+epilogue, 0);
+ n_mixed = vector(prologue+num+epilogue, 0);
+
+ for (int i = 0; i < prologue; i++)
+ {
+ if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
+ n_mixed[i]++;
+ else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
+ n_forward[i]++;
+ else if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
+ n_backward[i]++;
+ else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
+ n_static[i]++;
+ }
//For each block, the minimum set of feedback variable is computed
// and the non-feedback variables are reordered to get
// a sub-recursive block without feedback variables
@@ -611,21 +628,88 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
//First we have the recursive equations conditional on feedback variables
- for (vector::iterator its = Reordered_Vertice.begin(); its != Reordered_Vertice.end(); its++)
+ for (int j = 0; j < 4; j++)
{
- equation_reordered[order] = tmp_equation_reordered[*its+prologue];
- variable_reordered[order] = tmp_variable_reordered[*its+prologue];
- order++;
+ for (vector::iterator its = Reordered_Vertice.begin(); its != Reordered_Vertice.end(); its++)
+ {
+ bool something_done = false;
+ if (j == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
+ {
+ n_mixed[prologue+i]++;
+ something_done = true;
+ }
+ else if (j == 1 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
+ {
+ n_forward[prologue+i]++;
+ something_done = true;
+ }
+ else if (j == 2 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
+ {
+ n_backward[prologue+i]++;
+ something_done = true;
+ }
+ else if (j == 3 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
+ {
+ n_static[prologue+i]++;
+ something_done = true;
+ }
+ if (something_done)
+ {
+ equation_reordered[order] = tmp_equation_reordered[*its+prologue];
+ variable_reordered[order] = tmp_variable_reordered[*its+prologue];
+ order++;
+ }
+ }
}
components_set[i].second.second = Reordered_Vertice;
//Second we have the equations related to the feedback variables
- for (set::iterator its = feed_back_vertices.begin(); its != feed_back_vertices.end(); its++)
+ for (int j = 0; j < 4; j++)
{
- equation_reordered[order] = tmp_equation_reordered[v_index[vertex(*its, G)]+prologue];
- variable_reordered[order] = tmp_variable_reordered[v_index[vertex(*its, G)]+prologue];
- order++;
+ for (set::iterator its = feed_back_vertices.begin(); its != feed_back_vertices.end(); its++)
+ {
+ bool something_done = false;
+ if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
+ {
+ n_mixed[prologue+i]++;
+ something_done = true;
+ }
+ else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
+ {
+ n_forward[prologue+i]++;
+ something_done = true;
+ }
+ else if (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
+ {
+ n_backward[prologue+i]++;
+ something_done = true;
+ }
+ else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
+ {
+ n_static[prologue+i]++;
+ something_done = true;
+ }
+ if (something_done)
+ {
+ equation_reordered[order] = tmp_equation_reordered[v_index[vertex(*its, G)]+prologue];
+ variable_reordered[order] = tmp_variable_reordered[v_index[vertex(*its, G)]+prologue];
+ order++;
+ }
+ }
}
}
+
+ for (int i = 0; i < epilogue; i++)
+ {
+ if (variable_lag_lead[tmp_variable_reordered[prologue+num+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+num+i]].second != 0)
+ n_mixed[prologue+num+i]++;
+ else if (variable_lag_lead[tmp_variable_reordered[prologue+num+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+num+i]].second != 0)
+ n_forward[prologue+num+i]++;
+ else if (variable_lag_lead[tmp_variable_reordered[prologue+num+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+num+i]].second == 0)
+ n_backward[prologue+num+i]++;
+ else if (variable_lag_lead[tmp_variable_reordered[prologue+num+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+num+i]].second == 0)
+ n_static[prologue+num+i]++;
+ }
+
inv_equation_reordered = vector(nb_var);
inv_variable_reordered = vector(nb_var);
for (int i = 0; i < nb_var; i++)
@@ -665,7 +749,7 @@ ModelTree::printBlockDecomposition(const vector > &blocks) const
}
block_type_firstequation_size_mfs_t
-ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, const vector > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector &variable_reordered, const vector &equation_reordered)
+ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector &variable_reordered, const vector &equation_reordered, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed, vector, pair > > &block_col_type)
{
int i = 0;
int count_equ = 0, blck_count_simult = 0;
@@ -674,6 +758,10 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
block_type_firstequation_size_mfs_t block_type_size_mfs;
BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
int eq = 0;
+ unsigned int l_n_static = 0;
+ unsigned int l_n_forward = 0;
+ unsigned int l_n_backward = 0;
+ unsigned int l_n_mixed = 0;
for (i = 0; i < prologue+(int) blocks.size()+epilogue; i++)
{
int first_count_equ = count_equ;
@@ -736,6 +824,10 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
else
Simulation_Type = SOLVE_FORWARD_SIMPLE;
}
+ l_n_static = n_static[i];
+ l_n_forward = n_forward[i];
+ l_n_backward = n_backward[i];
+ l_n_mixed = n_mixed[i];
if (Blck_Size == 1)
{
if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
@@ -754,29 +846,35 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
BlockSimulationType c_Type = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.first;
int c_Size = (block_type_size_mfs[block_type_size_mfs.size()-1]).second.first;
int first_equation = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.second;
- block_type_size_mfs[block_type_size_mfs.size()-1] = make_pair(make_pair(c_Type, first_equation), make_pair(++c_Size, block_type_size_mfs[block_type_size_mfs.size()-1].second.second));
+ c_Size++;
+ block_type_size_mfs[block_type_size_mfs.size()-1] = make_pair(make_pair(c_Type, first_equation), make_pair(c_Size, c_Size));
if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
block_lag_lead[block_type_size_mfs.size()-1] = make_pair(Lag, Lead);
+ pair< pair< unsigned int, unsigned int>, pair > tmp = block_col_type[block_col_type.size()-1];
+ block_col_type[block_col_type.size()-1] = make_pair( make_pair(tmp.first.first+l_n_static, tmp.first.second+l_n_forward), make_pair(tmp.second.first+l_n_backward, tmp.second.second+l_n_mixed) );
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
+ block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
}
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
+ block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
}
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
+ block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
}
prev_Type = Simulation_Type;
eq += Blck_Size;
@@ -1015,7 +1113,7 @@ ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, ostream &output,
}
void
-ModelTree::compileTemporaryTerms(ostream &code_file, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const
+ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const
{
// Local var used to keep track of temp nodes already written
temporary_terms_t tt2;
@@ -1023,17 +1121,17 @@ ModelTree::compileTemporaryTerms(ostream &code_file, const temporary_terms_t &tt
it != tt.end(); it++)
{
FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
- fnumexpr.write(code_file);
- (*it)->compile(code_file, false, tt2, map_idx, dynamic, steady_dynamic);
+ fnumexpr.write(code_file, instruction_number);
+ (*it)->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic);
if (dynamic)
{
FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
- fstpt.write(code_file);
+ fstpt.write(code_file, instruction_number);
}
else
{
FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second));
- fstpst.write(code_file);
+ fstpst.write(code_file, instruction_number);
}
// Insert current node into tt2
tt2.insert(*it);
@@ -1108,7 +1206,7 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type)
}
void
-ModelTree::compileModelEquations(ostream &code_file, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
for (int eq = 0; eq < (int) equations.size(); eq++)
{
@@ -1116,7 +1214,7 @@ ModelTree::compileModelEquations(ostream &code_file, const temporary_terms_t &tt
expr_t lhs = eq_node->get_arg1();
expr_t rhs = eq_node->get_arg2();
FNUMEXPR_ fnumexpr(ModelEquation, eq);
- fnumexpr.write(code_file);
+ fnumexpr.write(code_file, instruction_number);
// Test if the right hand side of the equation is empty.
double vrhs = 1.0;
try
@@ -1129,20 +1227,20 @@ ModelTree::compileModelEquations(ostream &code_file, const temporary_terms_t &tt
if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
{
- lhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic);
- rhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic);
+ lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
+ rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FSTPR_ fstpr(eq);
- fstpr.write(code_file);
+ fstpr.write(code_file, instruction_number);
}
else // The right hand side of the equation is empty ==> residual=lhs;
{
- lhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic);
+ lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
FSTPR_ fstpr(eq);
- fstpr.write(code_file);
+ fstpr.write(code_file, instruction_number);
}
}
}
diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh
index 9297a1b22..47ad302d7 100644
--- a/preprocessor/ModelTree.hh
+++ b/preprocessor/ModelTree.hh
@@ -110,7 +110,7 @@ protected:
//! Writes temporary terms
void writeTemporaryTerms(const temporary_terms_t &tt, ostream &output, ExprNodeOutputType output_type) const;
//! Compiles temporary terms
- void compileTemporaryTerms(ostream &code_file, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const;
+ void compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const;
//! Adds informations for simulation in a binary file
void Write_Inf_To_Bin_File(const string &basename, int &u_count_int, bool &file_open, bool is_two_boundaries, int block_mfs) const;
@@ -120,7 +120,7 @@ protected:
//! Writes model equations
void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const;
//! Compiles model equations
- void compileModelEquations(ostream &code_file, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+ void compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
//! Writes LaTeX model file
void writeLatexModelFile(const string &filename, ExprNodeOutputType output_type) const;
@@ -167,9 +167,9 @@ protected:
//! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
equation_type_and_normalized_equation_t equationTypeDetermination(const map >, expr_t> &first_order_endo_derivatives, const vector &Index_Var_IM, const vector &Index_Equ_IM, int mfs) const;
//! Compute the block decomposition and for a non-recusive block find the minimum feedback set
- void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector &equation_reordered, vector &variable_reordered, vector > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector &inv_equation_reordered, vector &inv_variable_reordered) const;
+ void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector &equation_reordered, vector &variable_reordered, vector > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector &inv_equation_reordered, vector &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead_t, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed) const;
//! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
- block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, const vector > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector &variable_reordered, const vector &equation_reordered);
+ block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector &variable_reordered, const vector &equation_reordered, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed, vector, pair > > &block_col_type);
//! Determine the maximum number of lead and lag for the endogenous variable in a bloc
void getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector &equation_reordered, const vector &variable_reordered) const;
//! Print an abstract of the block structure of the model
@@ -189,6 +189,10 @@ protected:
virtual unsigned int getBlockFirstEquation(int block_number) const = 0;
//! Return the size of the block block_number
virtual unsigned int getBlockSize(int block_number) const = 0;
+ //! Return the number of exogenous variable in the block block_number
+ virtual unsigned int getBlockExoSize(int block_number) const = 0;
+ //! Return the number of colums in the jacobian matrix for exogenous variable in the block block_number
+ virtual unsigned int getBlockExoColSize(int block_number) const = 0;
//! Return the number of feedback variable of the block block_number
virtual unsigned int getBlockMfs(int block_number) const = 0;
//! Return the maximum lag in a block
@@ -207,11 +211,18 @@ protected:
virtual int getBlockEquationID(int block_number, int equation_number) const = 0;
//! Return the original number of variable variable_number belonging to the block block_number
virtual int getBlockVariableID(int block_number, int variable_number) const = 0;
+ //! Return the original number of the exogenous variable varexo_number belonging to the block block_number
+ virtual int getBlockVariableExoID(int block_number, int variable_number) const = 0;
//! Return the position of equation_number in the block number belonging to the block block_number
virtual int getBlockInitialEquationID(int block_number, int equation_number) const = 0;
//! Return the position of variable_number in the block number belonging to the block block_number
virtual int getBlockInitialVariableID(int block_number, int variable_number) const = 0;
-
+ //! Return the position of variable_number in the block number belonging to the block block_number
+ virtual int getBlockInitialExogenousID(int block_number, int variable_number) const = 0;
+ //! Return the position of the deterministic exogenous variable_number in the block number belonging to the block block_number
+ virtual int getBlockInitialDetExogenousID(int block_number, int variable_number) const = 0;
+ //! Return the position of the other endogenous variable_number in the block number belonging to the block block_number
+ virtual int getBlockInitialOtherEndogenousID(int block_number, int variable_number) const = 0;
public:
ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_arg);
//! Declare a node as an equation of the model
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 58a7328f5..d060556c5 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -46,28 +46,28 @@ StaticModel::StaticModel(SymbolTable &symbol_table_arg,
}
void
-StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_t &map_idx) const
+StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx) const
{
first_derivatives_t::const_iterator it = first_derivatives.find(make_pair(eq, symbol_table.getID(eEndogenous, symb_id)));
if (it != first_derivatives.end())
- (it->second)->compile(code_file, false, temporary_terms, map_idx, false, false);
+ (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
else
{
FLDZ_ fldz;
- fldz.write(code_file);
+ fldz.write(code_file, instruction_number);
}
}
void
-StaticModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, map_idx_t &map_idx) const
+StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx) const
{
map >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
if (it != first_chain_rule_derivatives.end())
- (it->second)->compile(code_file, false, temporary_terms, map_idx, false, false);
+ (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
else
{
FLDZ_ fldz;
- fldz.write(code_file);
+ fldz.write(code_file, instruction_number);
}
}
@@ -406,6 +406,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
ostringstream tmp_output;
ofstream code_file;
+ unsigned int instruction_number = 0;
bool file_open = false;
string main_name = file_name;
@@ -424,8 +425,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
//Temporary variables declaration
FDIMT_ fdimt(temporary_terms.size());
- fdimt.write(code_file);
-
+ fdimt.write(code_file, instruction_number);
FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
SOLVE_FORWARD_COMPLETE,
0,
@@ -436,22 +436,22 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
symbol_table.endo_nbr(),
0,
0,
- u_count_int
+ u_count_int,
+ symbol_table.endo_nbr()
);
- fbeginblock.write(code_file);
-
+ fbeginblock.write(code_file, instruction_number);
// Add a mapping form node ID to temporary terms order
int j = 0;
for (temporary_terms_t::const_iterator it = temporary_terms.begin();
it != temporary_terms.end(); it++)
map_idx[(*it)->idx] = j++;
- compileTemporaryTerms(code_file, temporary_terms, map_idx, false, false);
+ compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, false, false);
- compileModelEquations(code_file, temporary_terms, map_idx, false, false);
+ compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, false, false);
FENDEQU_ fendequ;
- fendequ.write(code_file);
+ fendequ.write(code_file, instruction_number);
vector > > derivatives;
derivatives.resize(symbol_table.endo_nbr());
@@ -467,46 +467,46 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
int symb = getSymbIDByDerivID(deriv_id);
unsigned int var = symbol_table.getTypeSpecificID(symb);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
- fnumexpr.write(code_file);
+ fnumexpr.write(code_file, instruction_number);
if (!derivatives[eq].size())
derivatives[eq].clear();
derivatives[eq].push_back(make_pair(var, count_u));
- d1->compile(code_file, false, temporary_terms, map_idx, false, false);
+ d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
FSTPSU_ fstpsu(count_u);
- fstpsu.write(code_file);
+ fstpsu.write(code_file, instruction_number);
count_u++;
}
}
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
FLDR_ fldr(i);
- fldr.write(code_file);
+ fldr.write(code_file, instruction_number);
for(vector >::const_iterator it = derivatives[i].begin();
it != derivatives[i].end(); it++)
{
FLDSU_ fldsu(it->second);
- fldsu.write(code_file);
+ fldsu.write(code_file, instruction_number);
FLDSV_ fldsv(eEndogenous, it->first);
- fldsv.write(code_file);
+ fldsv.write(code_file, instruction_number);
FBINARY_ fbinary(oTimes);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
if (it != derivatives[i].begin())
{
FBINARY_ fbinary(oPlus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
}
}
FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FSTPSU_ fstpsu(i);
- fstpsu.write(code_file);
+ fstpsu.write(code_file, instruction_number);
}
FENDBLOCK_ fendblock;
- fendblock.write(code_file);
+ fendblock.write(code_file, instruction_number);
FEND_ fend;
- fend.write(code_file);
+ fend.write(code_file, instruction_number);
code_file.close();
}
@@ -528,6 +528,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
string tmp_s;
ostringstream tmp_output;
ofstream code_file;
+ unsigned int instruction_number = 0;
expr_t lhs = NULL, rhs = NULL;
BinaryOpNode *eq_node;
Uff Uf[symbol_table.endo_nbr()];
@@ -546,7 +547,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
//Temporary variables declaration
FDIMT_ fdimt(temporary_terms.size());
- fdimt.write(code_file);
+ fdimt.write(code_file, instruction_number);
for (unsigned int block = 0; block < getNbBlocks(); block++)
{
@@ -554,7 +555,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
if (block > 0)
{
FENDBLOCK_ fendblock;
- fendblock.write(code_file);
+ fendblock.write(code_file, instruction_number);
}
int count_u;
int u_count_int = 0;
@@ -580,9 +581,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
symbol_table.endo_nbr(),
0,
0,
- u_count_int
+ u_count_int,
+ symbol_table.endo_nbr()
);
- fbeginblock.write(code_file);
+ fbeginblock.write(code_file, instruction_number);
// The equations
for (i = 0; i < (int) block_size; i++)
@@ -596,10 +598,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
it != v_temporary_terms[block][i].end(); it++)
{
FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
- fnumexpr.write(code_file);
- (*it)->compile(code_file, false, tt2, map_idx, false, false);
+ fnumexpr.write(code_file, instruction_number);
+ (*it)->compile(code_file, instruction_number, false, tt2, map_idx, false, false);
FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second));
- fstpst.write(code_file);
+ fstpst.write(code_file, instruction_number);
// Insert current node into tt2
tt2.insert(*it);
}
@@ -615,23 +617,23 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
equ_type = getBlockEquationType(block, i);
{
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
- fnumexpr.write(code_file);
+ fnumexpr.write(code_file, instruction_number);
}
if (equ_type == E_EVALUATE)
{
eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- rhs->compile(code_file, false, temporary_terms, map_idx, false, false);
- lhs->compile(code_file, true, temporary_terms, map_idx, false, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
+ lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
}
else if (equ_type == E_EVALUATE_S)
{
eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- rhs->compile(code_file, false, temporary_terms, map_idx, false, false);
- lhs->compile(code_file, true, temporary_terms, map_idx, false, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
+ lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
}
break;
case SOLVE_BACKWARD_COMPLETE:
@@ -646,22 +648,22 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
default:
end:
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
- fnumexpr.write(code_file);
+ fnumexpr.write(code_file, instruction_number);
eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- lhs->compile(code_file, false, temporary_terms, map_idx, false, false);
- rhs->compile(code_file, false, temporary_terms, map_idx, false, false);
+ lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FSTPR_ fstpr(i - block_recursive);
- fstpr.write(code_file);
+ fstpr.write(code_file, instruction_number);
}
}
FENDEQU_ fendequ;
- fendequ.write(code_file);
+ fendequ.write(code_file, instruction_number);
// The Jacobian if we have to solve the block
if (simulation_type != EVALUATE_BACKWARD
&& simulation_type != EVALUATE_FORWARD)
@@ -672,12 +674,12 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
case SOLVE_FORWARD_SIMPLE:
{
FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
- fnumexpr.write(code_file);
+ fnumexpr.write(code_file, instruction_number);
}
- compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx);
+ compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx);
{
FSTPG_ fstpg(0);
- fstpg.write(code_file);
+ fstpg.write(code_file, instruction_number);
}
break;
@@ -706,10 +708,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
Uf[eqr].Ufl->u = count_u;
Uf[eqr].Ufl->var = varr;
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
- fnumexpr.write(code_file);
- compileChainRuleDerivative(code_file, eqr, varr, 0, map_idx);
+ fnumexpr.write(code_file, instruction_number);
+ compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx);
FSTPSU_ fstpsu(count_u);
- fstpsu.write(code_file);
+ fstpsu.write(code_file, instruction_number);
count_u++;
}
}
@@ -718,24 +720,24 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
if (i >= (int) block_recursive)
{
FLDR_ fldr(i-block_recursive);
- fldr.write(code_file);
+ fldr.write(code_file, instruction_number);
FLDZ_ fldz;
- fldz.write(code_file);
+ fldz.write(code_file, instruction_number);
v = getBlockEquationID(block, i);
for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
{
FLDSU_ fldsu(Uf[v].Ufl->u);
- fldsu.write(code_file);
+ fldsu.write(code_file, instruction_number);
FLDSV_ fldsv(eEndogenous, Uf[v].Ufl->var);
- fldsv.write(code_file);
+ fldsv.write(code_file, instruction_number);
FBINARY_ fbinary(oTimes);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FCUML_ fcuml;
- fcuml.write(code_file);
+ fcuml.write(code_file, instruction_number);
}
Uf[v].Ufl = Uf[v].Ufl_First;
while (Uf[v].Ufl)
@@ -745,10 +747,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
Uf[v].Ufl = Uf[v].Ufl_First;
}
FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
+ fbinary.write(code_file, instruction_number);
FSTPSU_ fstpsu(i - block_recursive);
- fstpsu.write(code_file);
+ fstpsu.write(code_file, instruction_number);
}
}
@@ -759,9 +761,9 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
}
}
FENDBLOCK_ fendblock;
- fendblock.write(code_file);
+ fendblock.write(code_file, instruction_number);
FEND_ fend;
- fend.write(code_file);
+ fend.write(code_file, instruction_number);
code_file.close();
}
@@ -858,6 +860,7 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
if (block)
{
jacob_map_t contemporaneous_jacobian, static_jacobian;
+ vector