From 623bde1b56612f16d5aeec6059107a6082c06aa5 Mon Sep 17 00:00:00 2001
From: sebastien <sebastien@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Fri, 17 Oct 2008 17:21:22 +0000
Subject: [PATCH] trunk preprocessor: removed LCC_COMPILER and GCC_COMPILER
options of SPARSE_DLL mode
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2177 ac1d8469-bf42-47a9-8791-bf33cf982152
---
ComputingTasks.cc | 22 +-
DynareBison.yy | 30 +-
DynareFlex.ll | 3 -
ModelTree.cc | 1871 ++++++++++---------------------------
ParsingDriver.cc | 16 +-
include/ComputingTasks.hh | 4 +-
include/ModelTree.hh | 17 +-
include/ParsingDriver.hh | 7 -
8 files changed, 490 insertions(+), 1480 deletions(-)
diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index c357d32f..612c56ec 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -104,9 +104,8 @@ SimulStatement::writeOutput(ostream &output, const string &basename) const
}
SimulSparseStatement::SimulSparseStatement(const OptionsList &options_list_arg,
- int compiler_arg, int mode_arg) :
+ int mode_arg) :
options_list(options_list_arg),
- compiler(compiler_arg),
mode(mode_arg)
{
}
@@ -130,24 +129,7 @@ SimulSparseStatement::writeOutput(ostream &output, const string &basename) const
output << " end\n";
output << "end\n";
if(mode==eSparseDLLMode)
- {
- if(compiler!=NO_COMPILE)
- {
- output << "disp('compiling...');\n";
- output << "t0=clock;\n";
- if (compiler == 0)
- output << "mex " << basename << "_dynamic.c;\n";
- else
- output << "mex " << basename << "_dynamic.cc;\n";
- output << "disp(['compiling time: ' num2str(etime(clock,t0))]);\n";
- output << "oo_.endo_simul=" << basename << "_dynamic;\n";
- }
- else
- {
- output << "oo_.endo_simul=simulate;\n";
- output << "clear simulate.dll;\n";
- }
- }
+ output << "oo_.endo_simul=simulate;\n";
else
{
//output << "oo_.endo_simul=" << basename << "_dynamic();\n";
diff --git a/DynareBison.yy b/DynareBison.yy
index 89fcf57d..ccdcaf96 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -90,19 +90,19 @@ class ParsingDriver;
%token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS
%token <string_val> FLOAT_NUMBER
%token FORECAST
-%token GAMMA_PDF GAUSSIAN_ELIMINATION GCC_COMPILER GMRES GRAPH
+%token GAMMA_PDF GAUSSIAN_ELIMINATION GMRES GRAPH
%token HISTVAL HP_FILTER HP_NGRID
%token INITVAL INITVAL_FILE
%token <string_val> INT_NUMBER
%token INV_GAMMA1_PDF INV_GAMMA2_PDF IRF
%token KALMAN_ALGO KALMAN_TOL
-%token LAPLACE LCC_COMPILER LIK_ALGO LIK_INIT LINEAR LOAD_MH_FILE LOGLINEAR LU
+%token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_MH_FILE LOGLINEAR LU
%token MARKOWITZ MARGINAL_DENSITY MAX
%token METHOD MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN
%token MODE_CHECK MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS
%token MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER
%token <string_val> NAME
-%token NO_COMPILER NOBS NOCONSTANT NOCORR NODIAGNOSTIC NOFUNCTIONS
+%token NOBS NOCONSTANT NOCORR NODIAGNOSTIC NOFUNCTIONS
%token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF
%token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS
%token PARAMETERS PERIODS PLANNER_OBJECTIVE PREFILTER PRESAMPLE
@@ -393,27 +393,11 @@ histval_list : histval_list histval_elem
histval_elem : NAME '(' signed_integer ')' EQUAL expression ';' { driver.hist_val($1, $3, $6); };
-model_sparse_dll_options_list : model_sparse_dll_options_list COMMA model_sparse_dll_options
- | model_sparse_dll_options
+model_sparse_options_list : model_sparse_options_list COMMA model_sparse_options
+ | model_sparse_options
;
-model_sparse_options_list : model_sparse_options_list COMMA model_sparse_common_options
- | model_sparse_common_options
- ;
-
-model_sparse_dll_options : model_compiler_options
- | model_sparse_common_options
- ;
-
-model_compiler_options : LCC_COMPILER
- { driver.init_compiler(0); }
- | GCC_COMPILER
- { driver.init_compiler(1); }
- | NO_COMPILER
- { driver.init_compiler(2); }
- ;
-
-model_sparse_common_options : o_cutoff
+model_sparse_options : o_cutoff
| o_markowitz
;
@@ -423,7 +407,7 @@ model : MODEL ';' { driver.begin_model(); }
equation_list END { driver.reset_data_tree(); }
| MODEL '(' USE_DLL ')' ';' { driver.begin_model(); driver.use_dll(); }
equation_list END { driver.reset_data_tree(); }
- | MODEL '(' SPARSE_DLL COMMA model_sparse_dll_options_list ')'
+ | MODEL '(' SPARSE_DLL COMMA model_sparse_options_list ')'
{ driver.begin_model(); driver.sparse_dll(); } ';'
equation_list END { driver.reset_data_tree(); }
| MODEL '(' SPARSE_DLL ')' { driver.begin_model(); driver.sparse_dll(); } ';'
diff --git a/DynareFlex.ll b/DynareFlex.ll
index 30ee679c..a7df257f 100644
--- a/DynareFlex.ll
+++ b/DynareFlex.ll
@@ -283,9 +283,6 @@ int sigma_e = 0;
<DYNARE_STATEMENT,DYNARE_BLOCK>bicgstab {return token::BICGSTAB;}
<DYNARE_STATEMENT,DYNARE_BLOCK>sparse {return token::SPARSE;}
<DYNARE_STATEMENT,DYNARE_BLOCK>sparse_dll {return token::SPARSE_DLL;}
-<DYNARE_STATEMENT,DYNARE_BLOCK>gcc_compiler {return token::GCC_COMPILER;}
-<DYNARE_STATEMENT,DYNARE_BLOCK>lcc_compiler {return token::LCC_COMPILER;}
-<DYNARE_STATEMENT,DYNARE_BLOCK>no_compiler {return token::NO_COMPILER;}
<DYNARE_STATEMENT,DYNARE_BLOCK>linear {return token::LINEAR;}
<DYNARE_STATEMENT,DYNARE_BLOCK>[,] {return token::COMMA;}
<DYNARE_STATEMENT,DYNARE_BLOCK>[:] {return Dynare::parser::token_type (yytext[0]);}
diff --git a/ModelTree.cc b/ModelTree.cc
index 86a0fd93..becc64be 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -32,7 +32,6 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg) :
DataTree(symbol_table_arg, num_constants_arg),
mode(eStandardMode),
- compiler(NO_COMPILE),
cutoff(1e-12),
markowitz(0.7),
new_SGE(true),
@@ -344,244 +343,6 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
/*EndNew*/
}
-void
-ModelTree::writeModelEquationsOrdered_C(ostream &output, Model_Block *ModelBlock) const
-{
- int i,j,k,m;
- string tmp_s;
- ostringstream tmp_output;
- NodeID lhs=NULL, rhs=NULL;
- BinaryOpNode *eq_node;
- bool OK, lhs_rhs_done, skip_the_head;
- ostringstream Uf[symbol_table.endo_nbr];
- map<NodeID, int> reference_count;
- int prev_Simulation_Type=-1;
- temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
- //----------------------------------------------------------------------
- //Temporary variables declaration
- OK=true;
- for(temporary_terms_type::const_iterator it = temporary_terms.begin();
- it != temporary_terms.end(); it++)
- {
- if (OK)
- OK=false;
- else
- tmp_output << ", ";
-
- (*it)->writeOutput(tmp_output, oCDynamicModel, temporary_terms);
-
- tmp_output << "[" << block_triangular.periods + variable_table.max_lag+variable_table.max_lead << "]";
- }
- if (tmp_output.str().length()>0)
- {
- output << "double " << tmp_output.str() << ";\n\n";
- }
- //For each block
- for(j = 0;j < ModelBlock->Size;j++)
- {
- //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
- if (ModelBlock->Block_List[j].Size==1)
- {
- lhs_rhs_done=true;
- eq_node = equations[ModelBlock->Block_List[j].Equation[0]];
- lhs = eq_node->arg1;
- rhs = eq_node->arg2;
- tmp_output.str("");
- lhs->writeOutput(tmp_output, oCDynamicModelSparseDLL, temporary_terms);
- }
- else
- lhs_rhs_done=false;
- if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type)
- && (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
- ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD
- ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
- ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FOREWARD_R ))
- skip_the_head=true;
- else
- skip_the_head=false;
- if (!skip_the_head)
- {
- if (j>0)
- output << "}\n\n";
- output << "void Dynamic" << j+1 << "(double *y, double *x, double *residual, double *g1, double *g2)\n";
- output << "{\n";
- output << " ////////////////////////////////////////////////////////////////////////\n" <<
- " //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) <<
- " //\n" <<
- " // Simulation type ";
- output << BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) << " //\n" <<
- " ////////////////////////////////////////////////////////////////////////\n";
-#ifdef CONDITION
- if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
- output << " longd condition[" << ModelBlock->Block_List[j].Size << "]; /*to improve condition*/\n";
-#endif
- }
- //The Temporary terms
- temporary_terms_type tt2;
- if (ModelBlock->Block_List[j].Temporary_terms->size())
- output << " //Temporary variables\n";
- i=0;
- for(temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin();
- it != ModelBlock->Block_List[j].Temporary_terms->end(); it++)
- {
- output << " ";
- (*it)->writeOutput(output, oCDynamicModelSparseDLL, temporary_terms);
- output << " = ";
- (*it)->writeOutput(output, oCDynamicModelSparseDLL, tt2);
- // Insert current node into tt2
- tt2.insert(*it);
- output << ";" << endl;
- i++;
- }
- // The equations
- for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0);
- string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ;
- output << " //equation " << ModelBlock->Block_List[j].Equation[i] << " variable : " <<
- sModel << " (" << ModelBlock->Block_List[j].Variable[i] << ")\n";
- if (!lhs_rhs_done)
- {
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
- lhs = eq_node->arg1;
- rhs = eq_node->arg2;
- tmp_output.str("");
- lhs->writeOutput(tmp_output, oCDynamicModelSparseDLL, temporary_terms);
- }
- output << " ";
- switch(ModelBlock->Block_List[j].Simulation_Type)
- {
- case EVALUATE_BACKWARD:
- case EVALUATE_FOREWARD:
- output << tmp_output.str();
- output << " = ";
- rhs->writeOutput(output, oCDynamicModelSparseDLL, temporary_terms);
- output << ";\n";
- break;
- case EVALUATE_BACKWARD_R:
- case EVALUATE_FOREWARD_R:
- rhs->writeOutput(output, oCDynamicModelSparseDLL, temporary_terms);
- output << " = ";
- lhs->writeOutput(output, oCDynamicModelSparseDLL, temporary_terms);
- output << ";\n";
- break;
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_FOREWARD_COMPLETE:
- Uf[ModelBlock->Block_List[j].Equation[i]] << " u[" << i << "] = residual[" << i << "]";
- goto end;
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- Uf[ModelBlock->Block_List[j].Equation[i]] << " u[" << i << "+Per_u_] = residual[" << i << "]";
- goto end;
- default:
- end:
- output << "residual[" << i << "] = (";
- output << tmp_output.str();
- output << ") - (";
- rhs->writeOutput(output, oCDynamicModelSparseDLL, temporary_terms);
- output << ");\n";
-#ifdef CONDITION
- if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
- output << " condition[" << i << "]=0;\n";
-#endif
- }
- }
- // The Jacobian if we have to solve the block
- if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD
- && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD
- && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R
- && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FOREWARD_R)
- {
- output << " /* Jacobian */\n";
- switch(ModelBlock->Block_List[j].Simulation_Type)
- {
- case SOLVE_BACKWARD_SIMPLE:
- case SOLVE_FOREWARD_SIMPLE:
- output << " g1[0]=";
- writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oCDynamicModelSparseDLL, temporary_terms);
- output << "; /* variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0])
- <<"(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) << ") " << ModelBlock->Block_List[j].Variable[0]
- << ", equation=" << ModelBlock->Block_List[j].Equation[0] << "*/\n";
- break;
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_FOREWARD_COMPLETE:
- m=ModelBlock->Block_List[j].Max_Lag;
- for(i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
- int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u[" << u << "]*y[Per_y_+" << var << "]";
- output << " u[" << u << "] = ";
- writeDerivative(output, eq, var, 0, oCDynamicModelSparseDLL, temporary_terms);
- output << "; // variable=" << symbol_table.getNameByID(eEndogenous, var)
- <<"(" << variable_table.getLag(variable_table.getSymbolID(var))<< ") " << var
- << ", equation=" << eq << "\n";
- }
- for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
- output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
- break;
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for(i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
- int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- if (k==0)
- Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u[" << u << "+Per_u_]*y[Per_y_+" << var << "]";
- else if (k>0)
- Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u[" << u << "+Per_u_]*y[(it_+" << k << ")*y_size+" << var << "]";
- else if (k<0)
- Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u[" << u << "+Per_u_]*y[(it_" << k << ")*y_size+" << var << "]";
- output << " u[" << u << "+Per_u_] = ";
- writeDerivative(output, eq, var, k, oCDynamicModelSparseDLL, temporary_terms);
- output << "; // variable=" << symbol_table.getNameByID(eEndogenous, var)
- <<"(" << k << ") " << var
- << ", equation=" << eq << "\n";
-#ifdef CONDITION
- output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
- output << " condition[" << eqr << "]=u[" << u << "+Per_u_];\n";
-#endif
- }
- }
- for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
-#ifdef CONDITION
- output << " if (fabs(condition[" << i << "])<fabs(u[" << i << "+Per_u_]))\n";
- output << " condition[" << i << "]=u[" << i << "+Per_u_];\n";
-#endif
- }
-#ifdef CONDITION
- for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for(i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
- int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- output << " u[" << u << "+Per_u_] /= condition[" << eqr << "];\n";
- }
- }
- for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
- output << " u[" << i << "+Per_u_] /= condition[" << i << "];\n";
-#endif
- break;
- }
- }
- prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
- }
- output << "}\n\n";
-}
-
-
-
void
ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock, const string &dynamic_basename) const
{
@@ -1923,97 +1684,6 @@ ModelTree::reform(const string name1) const
return (name);
}
-void
-ModelTree::writeSparseDLLDynamicHFile(const string &dynamic_basename) const
-{
- string filename;
- ofstream mDynamicModelFile;
- string tmp_s;
- int i, j;
-
- if (compiler == LCC_COMPILE)
- filename = dynamic_basename + ".h";
- else
- filename = dynamic_basename + ".hh";
- mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
- if (!mDynamicModelFile.is_open())
- {
- cout << "ModelTree::Open : Error : Can't open file " << filename
- << ".h for writing\n";
- exit(-1);
- }
- filename.erase(filename.end() - 2, filename.end());
- tmp_s = filename;
- j = tmp_s.size();
- for(i = 0;i < j;i++)
- if ((tmp_s[i] == '\\') || (tmp_s[i] == '.') || (tmp_s[i] == ':'))
- tmp_s[i] = '_';
- mDynamicModelFile << "#ifndef " << tmp_s << "\n";
- mDynamicModelFile << "#define " << tmp_s << "\n";
- if (compiler==LCC_COMPILE)
- {
- mDynamicModelFile << "typedef struct IM_compact\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " int size, u_init, u_finish, nb_endo;\n";
- mDynamicModelFile << " int *u, *Var, *Equ, *Var_Index, *Equ_Index, *Var_dyn_Index;\n";
- mDynamicModelFile << "} IM_compact;\n";
- mDynamicModelFile << "typedef struct Variable_l\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " int* Index;\n";
- mDynamicModelFile << "} Variable_l;\n";
- mDynamicModelFile << "typedef struct tBlock\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " int Size, Sized, Type, Max_Lead, Max_Lag, Simulation_Type, /*icc1_size,*/ Nb_Lead_Lag_Endo;\n";
- mDynamicModelFile << " int *Variable, *dVariable, *Equation/*, *icc1, *ics*/;\n";
- mDynamicModelFile << " int *variable_dyn_index, *variable_dyn_leadlag;\n";
- mDynamicModelFile << " IM_compact *IM_lead_lag;\n";
- mDynamicModelFile << "} tBlock;\n";
- mDynamicModelFile << "\n";
- mDynamicModelFile << "typedef struct tModel_Block\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " int Size;\n";
- mDynamicModelFile << " tBlock * List;\n";
- mDynamicModelFile << "} tModel_Block;\n";
- mDynamicModelFile << "\n";
- mDynamicModelFile << "double *u, slowc, max_res, res2, res1;\n";
- mDynamicModelFile << "double *params;\n";
- mDynamicModelFile << "int it_,Per_u_;\n";
- mDynamicModelFile << "bool cvg;\n";
- mDynamicModelFile << "int nb_row_x;\n";
- mDynamicModelFile << "int y_kmin, y_kmax,periods, x_size, y_size, u_size, maxit_;\n";
- mDynamicModelFile << "double *y=NULL, *x=NULL, *r=NULL, *g1=NULL, *g2=NULL, solve_tolf, dynaretol;\n";
- mDynamicModelFile << "pctimer_t t0, t1;\n";
- }
- else
- {
- mDynamicModelFile << "typedef struct IM_compact\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " int size, u_init, u_finish, nb_endo;\n";
- mDynamicModelFile << " int *u, *Var, *Equ, *Var_Index, *Equ_Index, *Var_dyn_Index;\n";
- mDynamicModelFile << "};\n";
- mDynamicModelFile << "typedef struct Variable_l\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " int* Index;\n";
- mDynamicModelFile << "};\n";
- mDynamicModelFile << "typedef struct tBlock\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " int Size, Sized, Type, Max_Lead, Max_Lag, Simulation_Type, /*icc1_size,*/ Nb_Lead_Lag_Endo;\n";
- mDynamicModelFile << " int *Variable, *dVariable, *Equation/*, *icc1, *ics*/;\n";
- mDynamicModelFile << " int *variable_dyn_index, *variable_dyn_leadlag;\n";
- mDynamicModelFile << " IM_compact *IM_lead_lag;\n";
- mDynamicModelFile << "};\n";
- mDynamicModelFile << "\n";
- mDynamicModelFile << "typedef struct tModel_Block\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " int Size;\n";
- mDynamicModelFile << " tBlock * List;\n";
- mDynamicModelFile << "};\n";
- mDynamicModelFile << "\n";
- }
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile.close();
-}
-
void
ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, const int &num,
int &u_count_int, bool &file_open) const
@@ -2291,1107 +1961,518 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
mStaticModelFile.close();
}
-
-/*void
-ModelTree::writeSparseDLLDynamicCFileAndBinFile(const string &dynamic_basename, const string &bin_basename, ExprNodeOutputType output_type) const*/
void
-ModelTree::writeSparseDynamicFileAndBinFile(const string &dynamic_basename, const string &bin_basename, ExprNodeOutputType output_type, const int mode) const
+ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string &basename, const int mode) const
{
- string filename, sp;
+ string sp;
ofstream mDynamicModelFile;
ostringstream tmp, tmp1, tmp_eq;
int prev_Simulation_Type;
SymbolicGaussElimination SGE;
bool OK;
- if(mode==eSparseDLLMode)
+ string filename = dynamic_basename + ".m";
+ mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
+ if (!mDynamicModelFile.is_open())
{
- if (compiler == LCC_COMPILE || compiler == GCC_COMPILE)
+ cerr << "Error: Can't open file " << filename << " for writing" << endl;
+ exit(-1);
+ }
+ mDynamicModelFile << "%\n";
+ mDynamicModelFile << "% " << filename << " : Computes dynamic model for Dynare\n";
+ mDynamicModelFile << "%\n";
+ mDynamicModelFile << "% Warning : this file is generated automatically by Dynare\n";
+ mDynamicModelFile << "% from model file (.mod)\n\n";
+ mDynamicModelFile << "%/\n";
+
+ int i, k, Nb_SGE=0;
+ bool printed = false, skip_head, open_par=false;
+ if (computeJacobian || computeJacobianExo || computeHessian)
+ {
+ mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n";
+ mDynamicModelFile << " global oo_ options_ M_ ;\n";
+ mDynamicModelFile << " g2=[];g3=[];\n";
+ //Temporary variables declaration
+ OK=true;
+ ostringstream tmp_output;
+ for(temporary_terms_type::const_iterator it = temporary_terms.begin();
+ it != temporary_terms.end(); it++)
{
- if (compiler == LCC_COMPILE)
- filename = dynamic_basename + ".c";
+ if (OK)
+ OK=false;
else
- filename = dynamic_basename + ".cc";
- mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
- if (!mDynamicModelFile.is_open())
+ tmp_output << " ";
+ (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
+ }
+ if (tmp_output.str().length()>0)
+ mDynamicModelFile << " global " << tmp_output.str() << " M_ ;\n";
+
+ mDynamicModelFile << " T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n";
+ tmp_output.clear();
+ OK=true;
+ for(temporary_terms_type::const_iterator it = temporary_terms.begin();
+ it != temporary_terms.end(); it++)
+ {
+ if (OK)
+ OK=false;
+ else
+ tmp_output << "=T_init;\n ";
+ (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
+ }
+ if (tmp_output.str().length()>0)
+ mDynamicModelFile << tmp_output.str() << "=T_init;\n";
+
+ mDynamicModelFile << " y_kmin=M_.maximum_lag;\n";
+ mDynamicModelFile << " y_kmax=M_.maximum_lead;\n";
+ mDynamicModelFile << " y_size=M_.endo_nbr;\n";
+ mDynamicModelFile << " if(length(varargin)>0)\n";
+ mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n";
+ //mDynamicModelFile << " global it_;\n";
+ mDynamicModelFile << " it_=varargin{3};\n";
+ //mDynamicModelFile << " g=zeros(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
+ mDynamicModelFile << " g1=spalloc(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1),y_size*y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
+ mDynamicModelFile << " Per_u_=0;\n";
+ mDynamicModelFile << " Per_y_=it_*y_size;\n";
+ mDynamicModelFile << " y=varargin{1};\n";
+ mDynamicModelFile << " ys=y(it_,:);\n";
+ /*mDynamicModelFile << " y1=varargin{1};\n";
+ mDynamicModelFile << " cnb_nz_elem=1;\n";
+ mDynamicModelFile << " for i = -y_kmin:y_kmax\n";
+ mDynamicModelFile << " nz_elem=find(M_.lead_lag_incidence(:,1+i+y_kmin));\n";
+ mDynamicModelFile << " nb_nz_elem=length(nz_elem);\n";
+ mDynamicModelFile << " y(it_+i, nz_elem)=y1(cnb_nz_elem:(cnb_nz_elem+nb_nz_elem));\n";
+ mDynamicModelFile << " if(i==0)\n";
+ mDynamicModelFile << " ys(nz_elem)=y(it_, nz_elem);\n";
+ mDynamicModelFile << " nz_elem_s=nz_elem;\n";
+ mDynamicModelFile << " end;\n";
+ mDynamicModelFile << " cnb_nz_elem=cnb_nz_elem+nb_nz_elem;\n";
+ mDynamicModelFile << " end;\n";*/
+ mDynamicModelFile << " x=varargin{2};\n";
+ prev_Simulation_Type=-1;
+ tmp.str("");
+ tmp_eq.str("");
+ for(i = 0;i < block_triangular.ModelBlock->Size;i++)
+ {
+ k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
+ if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) &&
+ ((prev_Simulation_Type==EVALUATE_FOREWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FOREWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R)
+ || (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)))
{
- cerr << "Error: Can't open file " << filename << " for writing" << endl;
- exit(-1);
+ mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
+ tmp_eq.str("");
+ mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
+ tmp.str("");
+ mDynamicModelFile << tmp1.str();
+ tmp1.str("");
}
- mDynamicModelFile << "/*\n";
- mDynamicModelFile << " * " << filename << " : Computes dynamic model for Dynare\n";
- mDynamicModelFile << " *\n";
- mDynamicModelFile << " * Warning : this file is generated automatically by Dynare\n";
- mDynamicModelFile << " * from model file (.mod)\n\n";
- mDynamicModelFile << " */\n";
- if (compiler==LCC_COMPILE)
+ //mDynamicModelFile << " y_index=[";
+ for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
- mDynamicModelFile << "#include <math.h>\n";
- mDynamicModelFile << "#include <stdio.h>\n";
- mDynamicModelFile << "#include <string.h>\n";
- mDynamicModelFile << "#include \"pctimer_h.h\"\n";
- mDynamicModelFile << "#include \"mex.h\" /* The Last include file*/\n";
- mDynamicModelFile << "#include \"" << dynamic_basename.c_str() << ".h\"\n";
- mDynamicModelFile << "#include \"simulate.h\"\n";
+ tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
}
+ //mDynamicModelFile << " ];\n";
+ if(k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)
+ {
+ if(i==block_triangular.ModelBlock->Size-1)
+ {
+ mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
+ tmp_eq.str("");
+ mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
+ tmp.str("");
+ mDynamicModelFile << tmp1.str();
+ tmp1.str("");
+ }
+ }
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
+ (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
+ skip_head=true;
else
+ skip_head=false;
+ switch(k)
{
- mDynamicModelFile << "#include \"" << dynamic_basename.c_str() << ".hh\"\n";
- mDynamicModelFile << "#include \"simulate.cc\"\n";
+ case EVALUATE_FOREWARD:
+ case EVALUATE_BACKWARD:
+ case EVALUATE_FOREWARD_R:
+ case EVALUATE_BACKWARD_R:
+ if(!skip_head)
+ {
+ tmp1 << " [y, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n";
+ tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n";
+ }
+ break;
+ case SOLVE_FOREWARD_SIMPLE:
+ case SOLVE_BACKWARD_SIMPLE:
+ mDynamicModelFile << " y_index=" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ";\n";
+ mDynamicModelFile << " [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, y_index, 1);\n";
+ mDynamicModelFile << " residual(y_index_eq)=r;\n";
+ break;
+ case SOLVE_FOREWARD_COMPLETE:
+ case SOLVE_BACKWARD_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_COMPLETE:
+ //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_size, 1);\n";
+ mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
+ mDynamicModelFile << " y_index = [" << tmp.str() << "];\n";
+ tmp.str("");
+ tmp_eq.str("");
+ int tmp_i=variable_table.max_lag+variable_table.max_lead+1;
+ mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n";
+ mDynamicModelFile << " y_index_c = y_index;\n";
+ tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1;
+ mDynamicModelFile << " for i=1:" << tmp_i-1 << ",\n";
+ mDynamicModelFile << " y_index_c = [y_index_c (y_index+i*y_size)];\n";
+ mDynamicModelFile << " end;\n";
+ mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n";
+ if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead)
+ mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga;\n";
+ else
+ mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";
+ mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
+ break;
}
- mDynamicModelFile << "//#define DEBUG\n";
+ prev_Simulation_Type=k;
}
- writeModelLocalVariables(mDynamicModelFile, oCDynamicModelSparseDLL);
- if (compiler==NO_COMPILE)
- writeModelEquationsCodeOrdered(dynamic_basename, block_triangular.ModelBlock, bin_basename, oCDynamicModelSparseDLL);
- else
- writeModelEquationsOrdered_C(mDynamicModelFile, block_triangular.ModelBlock);
- }
- else
- {
- filename = dynamic_basename + ".m";
- mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
- if (!mDynamicModelFile.is_open())
+ if(tmp1.str().length())
{
- cerr << "Error: Can't open file " << filename << " for writing" << endl;
- exit(-1);
+ mDynamicModelFile << tmp1.str();
+ tmp1.str("");
}
- mDynamicModelFile << "%\n";
- mDynamicModelFile << "% " << filename << " : Computes dynamic model for Dynare\n";
- mDynamicModelFile << "%\n";
- mDynamicModelFile << "% Warning : this file is generated automatically by Dynare\n";
- mDynamicModelFile << "% from model file (.mod)\n\n";
- mDynamicModelFile << "%/\n";
+ mDynamicModelFile << " varargout{1}=residual;\n";
+ mDynamicModelFile << " varargout{2}=g1;\n";
+ mDynamicModelFile << " return;\n";
+ mDynamicModelFile << " end;\n";
+ mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n";
+ mDynamicModelFile << " if(options_.simulation_method==0)\n";
+ mDynamicModelFile << " mthd='Sparse LU';\n";
+ mDynamicModelFile << " elseif(options_.simulation_method==2)\n";
+ mDynamicModelFile << " mthd='GMRES';\n";
+ mDynamicModelFile << " elseif(options_.simulation_method==3)\n";
+ mDynamicModelFile << " mthd='BICGSTAB';\n";
+ mDynamicModelFile << " else\n";
+ mDynamicModelFile << " mthd='UNKNOWN';\n";
+ mDynamicModelFile << " end;\n";
+ mDynamicModelFile << " disp (['-----------------------------------------------------']) ;\n";
+ mDynamicModelFile << " disp (['MODEL SIMULATION: (method=' mthd ')']) ;\n";
+ mDynamicModelFile << " fprintf('\\n') ;\n";
+ mDynamicModelFile << " periods=options_.periods;\n";
+ mDynamicModelFile << " maxit_=options_.maxit_;\n";
+ mDynamicModelFile << " solve_tolf=options_.solve_tolf;\n";
+ mDynamicModelFile << " y=oo_.endo_simul';\n";
+ mDynamicModelFile << " x=oo_.exo_simul;\n";
}
- int i, j, k, Nb_SGE=0;
- bool printed = false, skip_head, open_par=false;
- if (computeJacobian || computeJacobianExo || computeHessian)
+ prev_Simulation_Type=-1;
+ for(i = 0;i < block_triangular.ModelBlock->Size;i++)
{
- if (compiler!=NO_COMPILE || mode==eSparseMode)
+ k = block_triangular.ModelBlock->Block_List[i].Simulation_Type;
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
+ (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
+ skip_head=true;
+ else
+ skip_head=false;
+ if ((k == EVALUATE_FOREWARD || k == EVALUATE_FOREWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size))
{
- //mDynamicModelFile << "void Dynamic_Init(tModel_Block *Model_Block)\n";
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << "void Dynamic_Init()\n";
- mDynamicModelFile << " {\n";
- }
- else
+ if (!skip_head)
{
- mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n";
- mDynamicModelFile << " global oo_ options_ M_ ;\n";
- mDynamicModelFile << " g2=[];g3=[];\n";
- //Temporary variables declaration
- {
- OK=true;
- ostringstream tmp_output;
- for(temporary_terms_type::const_iterator it = temporary_terms.begin();
- it != temporary_terms.end(); it++)
- {
- if (OK)
- OK=false;
- else
- tmp_output << " ";
- (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
- }
- if (tmp_output.str().length()>0)
- {
- mDynamicModelFile << " global " << tmp_output.str() << " M_ ;\n";
- }
- }
- mDynamicModelFile << " T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n";
- {
- ostringstream tmp_output;
- OK=true;
- for(temporary_terms_type::const_iterator it = temporary_terms.begin();
- it != temporary_terms.end(); it++)
- {
- if (OK)
- OK=false;
- else
- tmp_output << "=T_init;\n ";
- (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
- }
- if (tmp_output.str().length()>0)
- {
- mDynamicModelFile << tmp_output.str() << "=T_init;\n";
- }
- }
- mDynamicModelFile << " y_kmin=M_.maximum_lag;\n";
- mDynamicModelFile << " y_kmax=M_.maximum_lead;\n";
- mDynamicModelFile << " y_size=M_.endo_nbr;\n";
- mDynamicModelFile << " if(length(varargin)>0)\n";
- mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n";
- //mDynamicModelFile << " global it_;\n";
- mDynamicModelFile << " it_=varargin{3};\n";
- //mDynamicModelFile << " g=zeros(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
- mDynamicModelFile << " g1=spalloc(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1),y_size*y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
- mDynamicModelFile << " Per_u_=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " y=varargin{1};\n";
- mDynamicModelFile << " ys=y(it_,:);\n";
- /*mDynamicModelFile << " y1=varargin{1};\n";
- mDynamicModelFile << " cnb_nz_elem=1;\n";
- mDynamicModelFile << " for i = -y_kmin:y_kmax\n";
- mDynamicModelFile << " nz_elem=find(M_.lead_lag_incidence(:,1+i+y_kmin));\n";
- mDynamicModelFile << " nb_nz_elem=length(nz_elem);\n";
- mDynamicModelFile << " y(it_+i, nz_elem)=y1(cnb_nz_elem:(cnb_nz_elem+nb_nz_elem));\n";
- mDynamicModelFile << " if(i==0)\n";
- mDynamicModelFile << " ys(nz_elem)=y(it_, nz_elem);\n";
- mDynamicModelFile << " nz_elem_s=nz_elem;\n";
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " cnb_nz_elem=cnb_nz_elem+nb_nz_elem;\n";
- mDynamicModelFile << " end;\n";*/
- mDynamicModelFile << " x=varargin{2};\n";
- prev_Simulation_Type=-1;
- tmp.str("");
- tmp_eq.str("");
- for(i = 0;i < block_triangular.ModelBlock->Size;i++)
- {
-
- k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) &&
- ((prev_Simulation_Type==EVALUATE_FOREWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FOREWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R)
- || (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R)))
- {
- mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
- tmp_eq.str("");
- mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
- tmp.str("");
- mDynamicModelFile << tmp1.str();
- tmp1.str("");
- }
- //mDynamicModelFile << " y_index=[";
- for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
- {
- tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
- tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
- }
- //mDynamicModelFile << " ];\n";
- if(k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R
- )
- {
- if(i==block_triangular.ModelBlock->Size-1)
- {
- mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
- tmp_eq.str("");
- mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
- tmp.str("");
- mDynamicModelFile << tmp1.str();
- tmp1.str("");
- }
- }
- if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
- (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
- skip_head=true;
- else
- skip_head=false;
- switch(k)
- {
- case EVALUATE_FOREWARD:
- case EVALUATE_BACKWARD:
- case EVALUATE_FOREWARD_R:
- case EVALUATE_BACKWARD_R:
- if(!skip_head)
- {
- tmp1 << " [y, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n";
- tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n";
- }
- break;
- case SOLVE_FOREWARD_SIMPLE:
- case SOLVE_BACKWARD_SIMPLE:
- mDynamicModelFile << " y_index=" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ";\n";
- mDynamicModelFile << " [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, y_index, 1);\n";
- mDynamicModelFile << " residual(y_index_eq)=r;\n";
- break;
- case SOLVE_FOREWARD_COMPLETE:
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_size, 1);\n";
- mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
- mDynamicModelFile << " y_index = [" << tmp.str() << "];\n";
- tmp.str("");
- tmp_eq.str("");
- int tmp_i=variable_table.max_lag+variable_table.max_lead+1;
- mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n";
- mDynamicModelFile << " y_index_c = y_index;\n";
- tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1;
- mDynamicModelFile << " for i=1:" << tmp_i-1 << ",\n";
- mDynamicModelFile << " y_index_c = [y_index_c (y_index+i*y_size)];\n";
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n";
- if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead)
- mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga;\n";
- else
- mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";
- mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
- break;
- }
- prev_Simulation_Type=k;
- }
- if(tmp1.str().length())
+ if (open_par)
{
- mDynamicModelFile << tmp1.str();
- tmp1.str("");
+ mDynamicModelFile << " end\n";
}
- mDynamicModelFile << " varargout{1}=residual;\n";
- mDynamicModelFile << " varargout{2}=g1;\n";
- mDynamicModelFile << " return;\n";
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n";
- mDynamicModelFile << " if(options_.simulation_method==0)\n";
- mDynamicModelFile << " mthd='Sparse LU';\n";
- mDynamicModelFile << " elseif(options_.simulation_method==2)\n";
- mDynamicModelFile << " mthd='GMRES';\n";
- mDynamicModelFile << " elseif(options_.simulation_method==3)\n";
- mDynamicModelFile << " mthd='BICGSTAB';\n";
- mDynamicModelFile << " else\n";
- mDynamicModelFile << " mthd='UNKNOWN';\n";
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " disp (['-----------------------------------------------------']) ;\n";
- mDynamicModelFile << " disp (['MODEL SIMULATION: (method=' mthd ')']) ;\n";
- mDynamicModelFile << " fprintf('\\n') ;\n";
- mDynamicModelFile << " periods=options_.periods;\n";
- mDynamicModelFile << " maxit_=options_.maxit_;\n";
- mDynamicModelFile << " solve_tolf=options_.solve_tolf;\n";
- mDynamicModelFile << " y=oo_.endo_simul';\n";
- mDynamicModelFile << " x=oo_.exo_simul;\n";
+ mDynamicModelFile << " Per_u_=0;\n";
+ mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
+ mDynamicModelFile << " Per_y_=it_*y_size;\n";
+ mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
+ mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n";
}
- prev_Simulation_Type=-1;
- for(i = 0;i < block_triangular.ModelBlock->Size;i++)
+ open_par=true;
+ }
+ else if ((k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (!skip_head)
{
- k = block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
- (k==EVALUATE_FOREWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FOREWARD_R || k==EVALUATE_BACKWARD_R))
- skip_head=true;
- else
- skip_head=false;
- if ((k == EVALUATE_FOREWARD || k == EVALUATE_FOREWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (!skip_head)
- {
- if(mode==eSparseDLLMode)
- {
- if (open_par)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- }
- mDynamicModelFile << " for(it_=y_kmin;it_<periods+y_kmin;it_++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << "#ifdef DEBUG\n";
- }
- else
- {
- if (open_par)
- {
- mDynamicModelFile << " end\n";
- }
- mDynamicModelFile << " Per_u_=0;\n";
- mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
- mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n";
- }
- }
- if(mode==eSparseDLLMode)
- for(j = 0;j < block_triangular.ModelBlock->Block_List[i].Size;j++)
- mDynamicModelFile << " mexPrintf(\"y[%d, %d]=%f \\n\",it_," << block_triangular.ModelBlock->Block_List[i].Variable[j] << ",double(y[it_," << block_triangular.ModelBlock->Block_List[i].Variable[j] << "]));\n";
- open_par=true;
- }
- else if ((k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size))
+ if (open_par)
{
- if (!skip_head)
- {
- if(mode==eSparseDLLMode)
- {
- if (open_par)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- }
- mDynamicModelFile << " for(it_=periods+y_kmin;it_>y_kmin;it_--)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << "#ifdef DEBUG\n";
- }
- else
- {
- if (open_par)
- {
- mDynamicModelFile << " end\n";
- }
- mDynamicModelFile << " Per_u_=0;\n";
- mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3);\n";
- }
- }
- if(mode==eSparseDLLMode)
- for(j = 0;j < block_triangular.ModelBlock->Block_List[i].Size;j++)
- mDynamicModelFile << " mexPrintf(\"y[%d, %d]=%f \\n\",it_," << block_triangular.ModelBlock->Block_List[i].Variable[j] << ",double(y[it_," << block_triangular.ModelBlock->Block_List[i].Variable[j] << "]));\n";
- open_par=true;
+ mDynamicModelFile << " end\n";
}
- else if ((k == SOLVE_FOREWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (open_par)
- {
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- //if (!skip_head)
- mDynamicModelFile << " end\n";
- }
- }
- open_par=false;
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << " g1=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " r=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " for(it_=y_kmin;it_<periods+y_kmin;it_++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " cvg=false;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " while(!((cvg)||(iter>maxit_)))\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << " y[Per_y_+" << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] += -r[0]/g1[0];\n";
- mDynamicModelFile << " cvg=((r[0]*r[0])<solve_tolf);\n";
- mDynamicModelFile << " iter++;\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " if (!cvg)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Convergence not achieved in block " << i << ", at time %d after %d iterations\\n\",it_,iter);\n";
- mDynamicModelFile << " mexErrMsgTxt(\"End of simulate\");\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << "#ifdef DEBUG\n";
- mDynamicModelFile << " mexPrintf(\"y[%d, %d]=%f \\n\",it_," << block_triangular.ModelBlock->Block_List[i].Variable[0] << ",y[it_," << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]);\n";
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mxFree(g1);\n";
- mDynamicModelFile << " mxFree(r);\n";
- }
- else
- {
- mDynamicModelFile << " g1=0;\n";
- mDynamicModelFile << " r=0;\n";
- mDynamicModelFile << " for it_=y_kmin+1:periods+y_kmin\n";
- mDynamicModelFile << " cvg=0;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n";
- mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n";
- mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n";
- mDynamicModelFile << " cvg=((r*r)<solve_tolf);\n";
- mDynamicModelFile << " iter=iter+1;\n";
- mDynamicModelFile << " end\n";
- mDynamicModelFile << " if cvg==0\n";
- mDynamicModelFile << " fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
- mDynamicModelFile << " return;\n";
- mDynamicModelFile << " end\n";
- mDynamicModelFile << " end\n";
- }
- }
- else if ((k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (open_par)
- {
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " end\n";
- }
- }
- open_par=false;
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << " g1=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " r=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " for(it_=periods+y_kmin;it_>y_kmin;it_--)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " cvg=false;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " while(!((cvg)||(iter>maxit_)))\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << " y[Per_y_+" << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] += -r[0]/g1[0];\n";
- mDynamicModelFile << " cvg=((r[0]*r[0])<solve_tolf);\n";
- mDynamicModelFile << " iter++;\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " if (!cvg)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Convergence not achieved in block " << i << ", at time %d after %d iterations\\n\",it_,iter);\n";
- mDynamicModelFile << " mexErrMsgTxt(\"End of simulate\");\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << "#ifdef DEBUG\n";
- mDynamicModelFile << " mexPrintf(\"y[%d, %d]=%f \\n\",it_," << block_triangular.ModelBlock->Block_List[i].Variable[0] << ",y[it_," << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]);\n";
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mxFree(g1);\n";
- mDynamicModelFile << " mxFree(r);\n";
- }
- else
- {
- mDynamicModelFile << " g1=0;\n";
- mDynamicModelFile << " r=0;\n";
- mDynamicModelFile << " for it_=periods+y_kmin:-1:y_kmin+1\n";
- mDynamicModelFile << " cvg=0;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n";
- mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n";
- mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n";
- mDynamicModelFile << " cvg=((r*r)<solve_tolf);\n";
- mDynamicModelFile << " iter=iter+1;\n";
- mDynamicModelFile << " end\n";
- mDynamicModelFile << " if cvg==0\n";
- mDynamicModelFile << " fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
- mDynamicModelFile << " return;\n";
- mDynamicModelFile << " end\n";
- mDynamicModelFile << " end\n";
- }
- }
- else if ((k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (open_par)
- {
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " end\n";
- }
- }
- open_par=false;
- if (!printed)
- {
- printed = true;
- }
- SGE.SGE_compute(block_triangular.ModelBlock, i, true, bin_basename, symbol_table.endo_nbr);
- Nb_SGE++;
+ mDynamicModelFile << " Per_u_=0;\n";
+ mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
+ mDynamicModelFile << " Per_y_=it_*y_size;\n";
+ mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3);\n";
+ }
+ open_par=true;
+ }
+ else if ((k == SOLVE_FOREWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (open_par)
+ mDynamicModelFile << " end\n";
+ open_par=false;
+ mDynamicModelFile << " g1=0;\n";
+ mDynamicModelFile << " r=0;\n";
+ mDynamicModelFile << " for it_=y_kmin+1:periods+y_kmin\n";
+ mDynamicModelFile << " cvg=0;\n";
+ mDynamicModelFile << " iter=0;\n";
+ mDynamicModelFile << " Per_y_=it_*y_size;\n";
+ mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n";
+ mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n";
+ mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n";
+ mDynamicModelFile << " cvg=((r*r)<solve_tolf);\n";
+ mDynamicModelFile << " iter=iter+1;\n";
+ mDynamicModelFile << " end\n";
+ mDynamicModelFile << " if cvg==0\n";
+ mDynamicModelFile << " fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
+ mDynamicModelFile << " return;\n";
+ mDynamicModelFile << " end\n";
+ mDynamicModelFile << " end\n";
+ }
+ else if ((k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (open_par)
+ mDynamicModelFile << " end\n";
+ open_par=false;
+ mDynamicModelFile << " g1=0;\n";
+ mDynamicModelFile << " r=0;\n";
+ mDynamicModelFile << " for it_=periods+y_kmin:-1:y_kmin+1\n";
+ mDynamicModelFile << " cvg=0;\n";
+ mDynamicModelFile << " iter=0;\n";
+ mDynamicModelFile << " Per_y_=it_*y_size;\n";
+ mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n";
+ mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n";
+ mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n";
+ mDynamicModelFile << " cvg=((r*r)<solve_tolf);\n";
+ mDynamicModelFile << " iter=iter+1;\n";
+ mDynamicModelFile << " end\n";
+ mDynamicModelFile << " if cvg==0\n";
+ mDynamicModelFile << " fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
+ mDynamicModelFile << " return;\n";
+ mDynamicModelFile << " end\n";
+ mDynamicModelFile << " end\n";
+ }
+ else if ((k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (open_par)
+ mDynamicModelFile << " end\n";
+ open_par=false;
+ if (!printed)
+ {
+ printed = true;
+ }
+ SGE.SGE_compute(block_triangular.ModelBlock, i, true, basename, symbol_table.endo_nbr);
+ Nb_SGE++;
#ifdef PRINT_OUT
- cout << "end of Gaussian elimination\n";
+ cout << "end of Gaussian elimination\n";
#endif
- mDynamicModelFile << " Read_file(\"" << reform(bin_basename) << "\",periods," <<
- block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr <<
- ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << ");\n";
- mDynamicModelFile << " g1=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " r=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- if (!block_triangular.ModelBlock->Block_List[i].is_linear)
- {
- mDynamicModelFile << " cvg=false;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " while(!((cvg)||(iter>maxit_)))\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " res2=0;\n";
- mDynamicModelFile << " res1=0;\n";
- mDynamicModelFile << " max_res=0;\n";
- mDynamicModelFile << " for(it_=y_kmin;it_<periods+y_kmin;it_++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Per_u_=(it_-y_kmin)*" << block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << " for(i=0;i<" << block_triangular.ModelBlock->Block_List[i].Size << ";i++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " if (max_res<fabs(r[i]))\n";
- mDynamicModelFile << " max_res=fabs(r[i]);\n";
- mDynamicModelFile << " res2+=r[i]*r[i];\n";
- mDynamicModelFile << " res1+=fabs(r[i]);\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " iter++;\n";
- mDynamicModelFile << " cvg=(max_res<solve_tolf);\n";
- mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", periods, true);\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " if (!cvg)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Convergence not achieved in block " << i << ", after %d iterations\\n\",iter);\n";
- mDynamicModelFile << " mexErrMsgTxt(\"End of simulate\");\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " for(it_=y_kmin;it_<periods+y_kmin;it_++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Per_u_=(it_-y_kmin)*" << block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, r, g1, g2, g3);\n";
- mDynamicModelFile << "#ifdef PRINT_OUT\n";
- mDynamicModelFile << " for(j=0;j<" << block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";j++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\" %f\",u[Per_u_+j]);\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mexPrintf(\"\\n\");\n";
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", periods, true);\n";
- }
- mDynamicModelFile << " mxFree(g1);\n";
- mDynamicModelFile << " mxFree(r);\n";
- mDynamicModelFile << " mxFree(u);\n";
- mDynamicModelFile << " //mexErrMsgTxt(\"Exit from Dynare\");\n";
- }
- else if ((k == SOLVE_FOREWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (open_par)
- {
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " end\n";
- }
- }
- open_par=false;
- if (!printed)
- {
- printed = true;
- }
- SGE.SGE_compute(block_triangular.ModelBlock, i, false, bin_basename, /*mod_param.endo_nbr*/symbol_table.endo_nbr);
- Nb_SGE++;
- mDynamicModelFile << " Read_file(\"" << reform(bin_basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";
- mDynamicModelFile << " g1=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " r=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " for(it_=y_kmin;it_<periods+y_kmin;it_++)\n";
- mDynamicModelFile << " {\n";
- if (!block_triangular.ModelBlock->Block_List[i].is_linear)
- {
- mDynamicModelFile << " cvg=false;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " while(!((cvg)||(iter>maxit_)))\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", 0, false);\n";
- mDynamicModelFile << " res2=0;\n";
- mDynamicModelFile << " res1=0;\n";
- mDynamicModelFile << " max_res=0;\n";
- mDynamicModelFile << " for(i=0;i<" << block_triangular.ModelBlock->Block_List[i].Size << ";i++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " if (max_res<fabs(r[i]))\n";
- mDynamicModelFile << " max_res=fabs(r[i]);\n";
- mDynamicModelFile << " res2+=r[i]*r[i];\n";
- mDynamicModelFile << " res1+=fabs(r[i]);\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " cvg=(max_res<solve_tolf);\n";
- mDynamicModelFile << " iter++;\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " if (!cvg)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Convergence not achieved in block " << i << ", at time %d after %d iterations\\n\",it_,iter);\n";
- mDynamicModelFile << " mexErrMsgTxt(\"End of simulate\");\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2, g3);\n";
- mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", 0, false);\n";
- }
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mxFree(g1);\n";
- mDynamicModelFile << " mxFree(r);\n";
- mDynamicModelFile << " mxFree(u);\n";
- }
- else if ((k == SOLVE_BACKWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (open_par)
- {
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " end\n";
- }
- }
- open_par=false;
- SGE.SGE_compute(block_triangular.ModelBlock, i, false, bin_basename, /*mod_param.endo_nbr*/symbol_table.endo_nbr);
- Nb_SGE++;
- mDynamicModelFile << " Read_file(\"" << reform(bin_basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";
- mDynamicModelFile << " g1=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " r=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " for(it_=periods+y_kmin;it_>y_kmin;it_--)\n";
- mDynamicModelFile << " {\n";
- if (!block_triangular.ModelBlock->Block_List[i].is_linear)
- {
- mDynamicModelFile << " cvg=false;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " while(!((cvg)||(iter>maxit_)))\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", 0, false);\n";
- mDynamicModelFile << " res2=0;\n";
- mDynamicModelFile << " for(i=0;i<" << block_triangular.ModelBlock->Block_List[i].Size << ";i++)\n";
- mDynamicModelFile << " res2+=r[i]*r[i];\n";
- mDynamicModelFile << " cvg=(res2<solve_tolf);\n";
- mDynamicModelFile << " iter++;\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " if (!cvg)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Convergence not achieved in block " << i << ", at time %d after %d iterations\\n\",it_,iter);\n";
- mDynamicModelFile << " mexErrMsgTxt(\"End of simulate\");\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", 0, false);\n";
- }
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mxFree(g1);\n";
- mDynamicModelFile << " mxFree(r);\n";
- mDynamicModelFile << " mxFree(u);\n";
- }
- else if ((k == SOLVE_TWO_BOUNDARIES_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
- {
- if (open_par)
- {
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " end\n";
- }
- }
- open_par=false;
- if (!printed)
- {
- printed = true;
- }
- Nb_SGE++;
- //cout << "new_SGE=" << new_SGE << "\n";
- if(mode==eSparseDLLMode)
- {
- if (new_SGE)
- {
- int u_count_int=0;
- Write_Inf_To_Bin_File(dynamic_basename, bin_basename, i, u_count_int,SGE.file_open);
- SGE.file_is_open();
- mDynamicModelFile << " u_count=" << u_count_int << "*periods;\n";
- mDynamicModelFile << " u_count_alloc = 2*u_count;\n";
- mDynamicModelFile << " u=(longd*)mxMalloc(u_count_alloc*sizeof(longd));\n";
- mDynamicModelFile << " memset(u, 0, u_count_alloc*sizeof(longd));\n";
- mDynamicModelFile << " u_count_init=" <<
- block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag +
- block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";\n";
- mDynamicModelFile << " Read_SparseMatrix(\"" << reform(bin_basename) << "\","
- << block_triangular.ModelBlock->Block_List[i].Size << ", periods, y_kmin, y_kmax"
- << ");\n";
- mDynamicModelFile << " u_count=" << u_count_int << "*(periods+y_kmax+y_kmin);\n";
- }
- else
- {
- SGE.SGE_compute(block_triangular.ModelBlock, i, true, bin_basename, /*mod_param.endo_nbr*/symbol_table.endo_nbr);
- mDynamicModelFile << " Read_file(\"" << reform(bin_basename) << "\",periods," <<
- block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr <<
- ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << ");\n";
- }
- mDynamicModelFile << " g1=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- mDynamicModelFile << " r=(double*)mxMalloc(" << block_triangular.ModelBlock->Block_List[i].Size << "*sizeof(double));\n";
- if (!block_triangular.ModelBlock->Block_List[i].is_linear)
- {
- mDynamicModelFile << " cvg=false;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " while(!((cvg)||(iter>maxit_)))\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " res2=0;\n";
- mDynamicModelFile << " res1=0;\n";
- mDynamicModelFile << " max_res=0;\n";
- mDynamicModelFile << " for(it_=y_kmin;it_<periods+y_kmin;it_++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Per_u_=(it_-y_kmin)*" << block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << " if (isnan(res1)||isinf(res1))\n";
- mDynamicModelFile << " break;\n";
- mDynamicModelFile << " for(i=0;i<" << block_triangular.ModelBlock->Block_List[i].Size << ";i++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " if (max_res<fabs(r[i]))\n";
- mDynamicModelFile << " max_res=fabs(r[i]);\n";
- mDynamicModelFile << " res2+=r[i]*r[i];\n";
- mDynamicModelFile << " res1+=fabs(r[i]);\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " cvg=(max_res<solve_tolf);\n";
- if (new_SGE)
- mDynamicModelFile << " simulate_NG1(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", periods, true, cvg);\n";
- else
- mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", periods, true);\n";
- mDynamicModelFile << " iter++;\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " if (!cvg)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Convergence not achieved in block " << i << ", after %d iterations\\n\",iter);\n";
- mDynamicModelFile << " mexErrMsgTxt(\"End of simulate\");\n";
- mDynamicModelFile << " }\n";
- }
- else
- {
- mDynamicModelFile << " for(it_=y_kmin;it_<periods+y_kmin;it_++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " Per_u_=(it_-y_kmin)*" << block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " Dynamic" << i + 1 << "(y, x, r, g1, g2);\n";
- mDynamicModelFile << "#ifdef PRINT_OUT\n";
- mDynamicModelFile << " for(j=0;j<" << block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ";j++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\" %f\",u[Per_u_+j]);\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mexPrintf(\"\\n\");\n";
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- if (new_SGE)
- mDynamicModelFile << " simulate_NG1(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", periods, true, cvg);\n";
- else
- mDynamicModelFile << " simulate(" << i << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr << ", it_, y_kmin, y_kmax," << block_triangular.ModelBlock->Block_List[i].Size << ", periods, true);\n";
- }
- mDynamicModelFile << " mxFree(g1);\n";
- mDynamicModelFile << " mxFree(r);\n";
- mDynamicModelFile << " mxFree(u);\n";
- mDynamicModelFile << " mxFree(index_vara);\n";
- mDynamicModelFile << " memset(direction,0,size_of_direction);\n";
- mDynamicModelFile << " //mexErrMsgTxt(\"Exit from Dynare\");\n";
- }
- else
- {
- mDynamicModelFile << " cvg=0;\n";
- mDynamicModelFile << " iter=0;\n";
- mDynamicModelFile << " Per_u_=0;\n";
- mDynamicModelFile << " y_index=[";
- for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
- {
- mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
- }
- mDynamicModelFile << " ];\n";
- mDynamicModelFile << " Blck_size=" << block_triangular.ModelBlock->Block_List[i].Size << ";\n";
- /*mDynamicModelFile << " if(options_.simulation_method==2 | options_.simulation_method==3),\n";
- mDynamicModelFile << " [r, g1]= " << bin_basename << "_static(y, x);\n";
- mDynamicModelFile << " [L1,U1] = lu(g1,1e-5);\n";
- mDynamicModelFile << " I = speye(periods);\n";
- mDynamicModelFile << " L1=kron(I,L1);\n";
- mDynamicModelFile << " U1=kron(I,U1);\n";
- mDynamicModelFile << " end;\n";*/
- mDynamicModelFile << " y_kmin_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lag << ";\n";
- mDynamicModelFile << " y_kmax_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lead << ";\n";
- mDynamicModelFile << " lambda=options_.slowc;\n";
- mDynamicModelFile << " correcting_factor=0.01;\n";
- mDynamicModelFile << " luinc_tol=1e-10;\n";
- mDynamicModelFile << " max_resa=1e100;\n";
- int nze, m;
- for(nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++)
- nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size;
- mDynamicModelFile << " Jacobian_Size=" << block_triangular.ModelBlock->Block_List[i].Size << "*(y_kmin+" << block_triangular.ModelBlock->Block_List[i].Max_Lead << " +periods);\n";
- mDynamicModelFile << " g1=spalloc( length(y_index)*periods, Jacobian_Size, " << nze << "*periods" << ");\n";
- /*mDynamicModelFile << " cpath=path;\n";
- mDynamicModelFile << " addpath(fullfile(matlabroot,'toolbox','matlab','sparfun'));\n";*/
- mDynamicModelFile << " bicgstabh=@bicgstab;\n";
- //mDynamicModelFile << " path(cpath);\n";
- mDynamicModelFile << sp << " reduced = 0;\n";
- //mDynamicModelFile << " functions(bicgstabh)\n";
- if (!block_triangular.ModelBlock->Block_List[i].is_linear)
- {
- sp=" ";
- mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n";
- }
- else
- {
- sp="";
- }
- mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, y_kmin, Blck_size, periods, g1, g2, g3);\n";
- mDynamicModelFile << sp << " g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);\n";
- mDynamicModelFile << sp << " b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';\n";
- mDynamicModelFile << sp << " if(~isreal(r))\n";
- mDynamicModelFile << sp << " max_res=(-(max(max(abs(r))))^2)^0.5;\n";
- mDynamicModelFile << sp << " else\n";
- mDynamicModelFile << sp << " max_res=max(max(abs(r)));\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " if(iter>0)\n";
- mDynamicModelFile << sp << " if(~isreal(max_res) | isnan(max_res) | (max_resa<max_res && iter>1))\n";
- mDynamicModelFile << sp << " if(isnan(max_res))\n";
- mDynamicModelFile << sp << " detJ=det(g1aa);\n";
- mDynamicModelFile << sp << " if(abs(detJ)<1e-7)\n";
- mDynamicModelFile << sp << " max_factor=max(max(abs(g1aa)));\n";
- mDynamicModelFile << sp << " ze_elem=sum(diag(g1aa)<options_.cutoff);\n";
- mDynamicModelFile << sp << " disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(options_.cutoff,'%f') ')']);\n";
- mDynamicModelFile << sp << " if(correcting_factor<max_factor)\n";
- mDynamicModelFile << sp << " correcting_factor=correcting_factor*4;\n";
- mDynamicModelFile << sp << " disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.']);\n";
- mDynamicModelFile << sp << " disp([' trying to correct the Jacobian matrix:']);\n";
- mDynamicModelFile << sp << " disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]);\n";
- mDynamicModelFile << sp << " dx = (g1aa+correcting_factor*speye(periods*Blck_size))\\ba- ya;\n";
- mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';\n";
- if (!block_triangular.ModelBlock->Block_List[i].is_linear)
- mDynamicModelFile << sp << " continue;\n";
- mDynamicModelFile << sp << " else\n";
- mDynamicModelFile << sp << " disp('The singularity of the jacobian matrix could not be corrected');\n";
- mDynamicModelFile << sp << " return;\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " elseif(lambda>1e-8)\n";
- mDynamicModelFile << sp << " lambda=lambda/2;\n";
- mDynamicModelFile << sp << " reduced = 1;\n";
- mDynamicModelFile << sp << " disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);\n";
- mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';\n";
- if (!block_triangular.ModelBlock->Block_List[i].is_linear)
- mDynamicModelFile << sp << " continue;\n";
- mDynamicModelFile << sp << " else\n";
- mDynamicModelFile << sp << " disp(['No convergence after ' num2str(iter,'%d') ' iterations']);\n";
- mDynamicModelFile << sp << " return;\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " else\n";
- mDynamicModelFile << sp << " if(lambda<1)\n";
- mDynamicModelFile << sp << " lambda=max(lambda*2, 1);\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " end;\n";
-
- mDynamicModelFile << sp << " ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)';\n";
- mDynamicModelFile << sp << " ya_save=ya;\n";
- mDynamicModelFile << sp << " g1aa=g1a;\n";
- mDynamicModelFile << sp << " ba=b;\n";
- mDynamicModelFile << sp << " max_resa=max_res;\n";
- mDynamicModelFile << sp << " if(options_.simulation_method==0),\n";
- mDynamicModelFile << sp << " dx = g1a\\b- ya;\n";
- mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
- mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
- mDynamicModelFile << sp << " elseif(options_.simulation_method==2),\n";
- mDynamicModelFile << sp << " flag1=1;\n";
- mDynamicModelFile << sp << " while(flag1>0)\n";
- mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n";
- mDynamicModelFile << sp << " [za,flag1] = gmres(g1a,b," << block_triangular.ModelBlock->Block_List[i].Size << ",1e-6," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n";
- mDynamicModelFile << sp << " if (flag1>0 | reduced)\n";
- mDynamicModelFile << sp << " if(flag1==1)\n";
- mDynamicModelFile << sp << " disp(['No convergence inside GMRES after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n";
- mDynamicModelFile << sp << " elseif(flag1==2)\n";
- mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n";
- mDynamicModelFile << sp << " elseif(flag1==3)\n";
- mDynamicModelFile << sp << " disp(['GMRES stagnated. (Two consecutive iterates were the same.)']);\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " luinc_tol = luinc_tol/10;\n";
- mDynamicModelFile << sp << " reduced = 0;\n";
- mDynamicModelFile << sp << " else\n";
- mDynamicModelFile << sp << " dx = za - ya;\n";
- mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
- mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " elseif(options_.simulation_method==3),\n";
- mDynamicModelFile << sp << " flag1=1;\n";
- mDynamicModelFile << sp << " while(flag1>0)\n";
- mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n";
- mDynamicModelFile << sp << " [za,flag1] = bicgstabh(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n";
- mDynamicModelFile << sp << " if (flag1>0 | reduced)\n";
- mDynamicModelFile << sp << " if(flag1==1)\n";
- mDynamicModelFile << sp << " disp(['No convergence inside BICGSTAB after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n";
- mDynamicModelFile << sp << " elseif(flag1==2)\n";
- mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n";
- mDynamicModelFile << sp << " elseif(flag1==3)\n";
- mDynamicModelFile << sp << " disp(['BICGSTAB stagnated. (Two consecutive iterates were the same.)']);\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " luinc_tol = luinc_tol/10;\n";
- mDynamicModelFile << sp << " reduced = 0;\n";
- mDynamicModelFile << sp << " else\n";
- mDynamicModelFile << sp << " dx = za - ya;\n";
- mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
- mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " end;\n";
- mDynamicModelFile << sp << " end;\n";
- if(!block_triangular.ModelBlock->Block_List[i].is_linear)
- {
- mDynamicModelFile << " cvg=(max_res<solve_tolf);\n";
- mDynamicModelFile << " iter=iter+1;\n";
- }
- mDynamicModelFile << " disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]);\n";
- if(!block_triangular.ModelBlock->Block_List[i].is_linear)
- {
- mDynamicModelFile << " end\n";
- mDynamicModelFile << " if (iter>maxit_)\n";
- mDynamicModelFile << " disp(['No convergence after ' num2str(iter,'%4d') ' iterations']);\n";
- mDynamicModelFile << " return;\n";
- mDynamicModelFile << " end;\n";
- }
- }
+ mDynamicModelFile << " Read_file(\"" << reform(basename) << "\",periods," <<
+ block_triangular.ModelBlock->Block_List[i].IM_lead_lag[block_triangular.ModelBlock->Block_List[i].Max_Lag + block_triangular.ModelBlock->Block_List[i].Max_Lead].u_finish + 1 << ", " << /*mod_param.endo_nbr*/symbol_table.endo_nbr <<
+ ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << ");\n";
+ cerr << "Not implemented block SOLVE_TWO_BOUNDARIES_COMPLETE" << endl;
+ exit(-1);
+ }
+ else if ((k == SOLVE_FOREWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (open_par)
+ mDynamicModelFile << " end\n";
+ open_par=false;
+ if (!printed)
+ {
+ printed = true;
+ }
+ SGE.SGE_compute(block_triangular.ModelBlock, i, false, basename, /*mod_param.endo_nbr*/symbol_table.endo_nbr);
+ Nb_SGE++;
+ mDynamicModelFile << " Read_file(\"" << reform(basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";
+ cerr << "Not implemented block SOLVE_FOREWARD_COMPLETE" << endl;
+ exit(-1);
+ }
+ else if ((k == SOLVE_BACKWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (open_par)
+ mDynamicModelFile << " end\n";
+ open_par=false;
+ SGE.SGE_compute(block_triangular.ModelBlock, i, false, basename, /*mod_param.endo_nbr*/symbol_table.endo_nbr);
+ Nb_SGE++;
+ mDynamicModelFile << " Read_file(\"" << reform(basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";
+ cerr << "Not implemented block SOLVE_BACKWARD_COMPLETE" << endl;
+ exit(-1);
+ }
+ else if ((k == SOLVE_TWO_BOUNDARIES_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ {
+ if (open_par)
+ mDynamicModelFile << " end\n";
+ open_par=false;
+ if (!printed)
+ {
+ printed = true;
+ }
+ Nb_SGE++;
+ //cout << "new_SGE=" << new_SGE << "\n";
- }
- prev_Simulation_Type=k;
+ mDynamicModelFile << " cvg=0;\n";
+ mDynamicModelFile << " iter=0;\n";
+ mDynamicModelFile << " Per_u_=0;\n";
+ mDynamicModelFile << " y_index=[";
+ for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
+ {
+ mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
}
- // Writing the gateway routine
- if(mode==eSparseDLLMode)
- {
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " }\n";
- }
- if(mode==eSparseMode)
+ mDynamicModelFile << " ];\n";
+ mDynamicModelFile << " Blck_size=" << block_triangular.ModelBlock->Block_List[i].Size << ";\n";
+ /*mDynamicModelFile << " if(options_.simulation_method==2 | options_.simulation_method==3),\n";
+ mDynamicModelFile << " [r, g1]= " << basename << "_static(y, x);\n";
+ mDynamicModelFile << " [L1,U1] = lu(g1,1e-5);\n";
+ mDynamicModelFile << " I = speye(periods);\n";
+ mDynamicModelFile << " L1=kron(I,L1);\n";
+ mDynamicModelFile << " U1=kron(I,U1);\n";
+ mDynamicModelFile << " end;\n";*/
+ mDynamicModelFile << " y_kmin_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lag << ";\n";
+ mDynamicModelFile << " y_kmax_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lead << ";\n";
+ mDynamicModelFile << " lambda=options_.slowc;\n";
+ mDynamicModelFile << " correcting_factor=0.01;\n";
+ mDynamicModelFile << " luinc_tol=1e-10;\n";
+ mDynamicModelFile << " max_resa=1e100;\n";
+ int nze, m;
+ for(nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++)
+ nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size;
+ mDynamicModelFile << " Jacobian_Size=" << block_triangular.ModelBlock->Block_List[i].Size << "*(y_kmin+" << block_triangular.ModelBlock->Block_List[i].Max_Lead << " +periods);\n";
+ mDynamicModelFile << " g1=spalloc( length(y_index)*periods, Jacobian_Size, " << nze << "*periods" << ");\n";
+ /*mDynamicModelFile << " cpath=path;\n";
+ mDynamicModelFile << " addpath(fullfile(matlabroot,'toolbox','matlab','sparfun'));\n";*/
+ mDynamicModelFile << " bicgstabh=@bicgstab;\n";
+ //mDynamicModelFile << " path(cpath);\n";
+ mDynamicModelFile << sp << " reduced = 0;\n";
+ //mDynamicModelFile << " functions(bicgstabh)\n";
+ if (!block_triangular.ModelBlock->Block_List[i].is_linear)
{
- if(open_par)
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " oo_.endo_simul = y';\n";
- mDynamicModelFile << "return;\n";
+ sp=" ";
+ mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n";
}
- if(mode==eSparseDLLMode)
+ else
{
- mDynamicModelFile << "/* The gateway routine */\n";
- mDynamicModelFile << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\n";
- mDynamicModelFile << "{\n";
- mDynamicModelFile << " mxArray *M_, *oo_, *options_;\n";
- mDynamicModelFile << " int i, row_y, col_y, row_x, col_x, nb_row_xd;\n";
- mDynamicModelFile << " double * pind ;\n";
- mDynamicModelFile << "\n";
- mDynamicModelFile << " /* Gets model parameters from global workspace of Matlab */\n";
- mDynamicModelFile << " M_ = mexGetVariable(\"global\",\"M_\");\n";
- mDynamicModelFile << " if (M_ == NULL )\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Global variable not found : \");\n";
- mDynamicModelFile << " mexErrMsgTxt(\"M_ \\n\");\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " /* Gets variables and parameters from global workspace of Matlab */\n";
- mDynamicModelFile << " oo_ = mexGetVariable(\"global\",\"oo_\");\n";
- mDynamicModelFile << " if (oo_ == NULL )\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Global variable not found : \");\n";
- mDynamicModelFile << " mexErrMsgTxt(\"oo_ \\n\");\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " options_ = mexGetVariable(\"global\",\"options_\");\n";
- mDynamicModelFile << " if (options_ == NULL )\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " mexPrintf(\"Global variable not found : \");\n";
- mDynamicModelFile << " mexErrMsgTxt(\"options_ \\n\");\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"params\")));\n";
- mDynamicModelFile << " double *yd, *xd;\n";
- mDynamicModelFile << " yd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,\"endo_simul\")));\n";
- mDynamicModelFile << " row_y=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,\"endo_simul\")));\n";
- mDynamicModelFile << " xd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,\"exo_simul\")));\n";
- mDynamicModelFile << " row_x=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,\"exo_simul\")));\n";
- mDynamicModelFile << " col_x=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,\"exo_simul\")));\n";
- if (compiler==GCC_COMPILE)
- {
- mDynamicModelFile << " y_kmin=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"maximum_lag\"))))));\n";
- mDynamicModelFile << " y_kmax=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"maximum_lead\"))))));\n";
- mDynamicModelFile << " y_decal=max(0,y_kmin-int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"maximum_endo_lag\")))))));\n";
- mDynamicModelFile << " periods=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"periods\"))))));\n";
- mDynamicModelFile << " maxit_=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"maxit_\"))))));\n";
- mDynamicModelFile << " slowc=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"slowc\")))));\n";
- mDynamicModelFile << " markowitz_c=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"markowitz\")))));\n";
- mDynamicModelFile << " nb_row_xd=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"exo_det_nbr\"))))));\n";
- }
- else
- {
- mDynamicModelFile << " y_kmin=(int)floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"maximum_lag\")))));\n";
- mDynamicModelFile << " y_kmax=(int)floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"maximum_lead\")))));\n";
- mDynamicModelFile << " y_decal=max(0,y_kmin-int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"maximum_endo_lag\")))))));\n";
- mDynamicModelFile << " periods=(int)floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"periods\")))));\n";
- mDynamicModelFile << " maxit_=(int)floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"maxit_\")))));\n";
- mDynamicModelFile << " slowc=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"slowc\")))));\n";
- mDynamicModelFile << " markowitz_c=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"markowitz\")))));\n";
- mDynamicModelFile << " nb_row_xd=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"exo_det_nbr\"))))));\n";
- }
- mDynamicModelFile << " mxArray *mxa=mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"fname\"));\n";
- mDynamicModelFile << " int buflen=mxGetM(mxa) * mxGetN(mxa) + 1;\n";
- mDynamicModelFile << " char *fname;\n";
- mDynamicModelFile << " fname=(char*)mxCalloc(buflen, sizeof(char));\n";
- mDynamicModelFile << " int status = mxGetString(mxa, fname, buflen);\n";
- mDynamicModelFile << " if (status != 0)\n";
- mDynamicModelFile << " mexWarnMsgTxt(\"Not enough space. Filename is truncated.\");\n";
- mDynamicModelFile << " mexPrintf(\"fname=%s\\n\",fname);\n";
- mDynamicModelFile << " col_y=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,\"endo_simul\")));;\n";
- mDynamicModelFile << " if (col_y<row_x)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " row_y=row_y/row_x;\n";
- mDynamicModelFile << " col_y=row_x;\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " solve_tolf=*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"dynatol\"))));\n";
- mDynamicModelFile << " size_of_direction=col_y*row_y*sizeof(longd);\n";
- mDynamicModelFile << " y=(longd*)mxMalloc(size_of_direction);\n";
- mDynamicModelFile << " ya=(longd*)mxMalloc(size_of_direction);\n";
- mDynamicModelFile << " direction=(longd*)mxMalloc(size_of_direction);\n";
- mDynamicModelFile << " memset(direction,0,size_of_direction);\n";
- mDynamicModelFile << " x=(longd*)mxMalloc(col_x*row_x*sizeof(longd));\n";
- mDynamicModelFile << " for(i=0;i<row_x*col_x;i++)\n";
- mDynamicModelFile << " x[i]=longd(xd[i]);\n";
- mDynamicModelFile << " for(i=0;i<row_y*col_y;i++)\n";
- mDynamicModelFile << " y[i]=longd(yd[i]);\n";
- mDynamicModelFile << " \n";
- mDynamicModelFile << " y_size=row_y;\n";
- mDynamicModelFile << " x_size=row_x;\n";
- mDynamicModelFile << " nb_row_x=row_x;\n";
- mDynamicModelFile << "#ifdef DEBUG\n";
- mDynamicModelFile << " for(j=0;j<periods+y_kmin+y_kmax;j++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " for(i=0;i<row_y;i++)\n";
- mDynamicModelFile << " mexPrintf(\"y[%d,%d]=%f \",j,i,y[j*y_size+i]);\n";
- mDynamicModelFile << " mexPrintf(\"\\n\");\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mexPrintf(\"\\n\");\n";
- mDynamicModelFile << " mexPrintf(\"x=%x\\n\",x);\n";
- mDynamicModelFile << " for(j=0;j<periods+y_kmin+y_kmax;j++)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " for(i=0;i<col_x;i++)\n";
- mDynamicModelFile << " mexPrintf(\"x[%d,%d]=%f \",j,i,x[i*x_size+j]);\n";
- mDynamicModelFile << " mexPrintf(\"\\n\");\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mexPrintf(\"x[1]=%f\\n\",x[1]);\n";
- mDynamicModelFile << "#endif\n";
- mDynamicModelFile << " /* Gets it_ from global workspace of Matlab */\n";
- mDynamicModelFile << " //it_ = (int) floor(mxGetScalar(mexGetVariable(\"global\", \"it_\")))-1;\n";
- mDynamicModelFile << " /* Call the C subroutines. */\n";
- mDynamicModelFile << " t0= pctimer();\n";
- mDynamicModelFile << " Dynamic_Init();\n";
- mDynamicModelFile << " t1= pctimer();\n";
- mDynamicModelFile << " mexPrintf(\"Simulation Time=%f milliseconds\\n\",1000*(t1-t0));\n";
- if (compiler==LCC_COMPILE )
- {
- mDynamicModelFile << " if (SaveCode)\n";
- mDynamicModelFile << " fclose(SaveCode);\n";
- }
- else
- {
- mDynamicModelFile << " if (SaveCode.is_open())\n";
- mDynamicModelFile << " SaveCode.close();\n";
- }
- mDynamicModelFile << " if (nlhs>0)\n";
- mDynamicModelFile << " {\n";
- mDynamicModelFile << " plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);\n";
- mDynamicModelFile << " pind = mxGetPr(plhs[0]);\n";
- mDynamicModelFile << " for(i=0;i<row_y*col_y;i++)\n";
- mDynamicModelFile << " pind[i]=y[i];\n";
- mDynamicModelFile << " }\n";
- mDynamicModelFile << " mxFree(x);\n";
- mDynamicModelFile << " mxFree(y);\n";
- mDynamicModelFile << " mxFree(ya);\n";
- mDynamicModelFile << " mxFree(direction);\n";
- mDynamicModelFile << "}\n";
+ sp="";
+ }
+ mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, y_kmin, Blck_size, periods, g1, g2, g3);\n";
+ mDynamicModelFile << sp << " g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);\n";
+ mDynamicModelFile << sp << " b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';\n";
+ mDynamicModelFile << sp << " if(~isreal(r))\n";
+ mDynamicModelFile << sp << " max_res=(-(max(max(abs(r))))^2)^0.5;\n";
+ mDynamicModelFile << sp << " else\n";
+ mDynamicModelFile << sp << " max_res=max(max(abs(r)));\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " if(iter>0)\n";
+ mDynamicModelFile << sp << " if(~isreal(max_res) | isnan(max_res) | (max_resa<max_res && iter>1))\n";
+ mDynamicModelFile << sp << " if(isnan(max_res))\n";
+ mDynamicModelFile << sp << " detJ=det(g1aa);\n";
+ mDynamicModelFile << sp << " if(abs(detJ)<1e-7)\n";
+ mDynamicModelFile << sp << " max_factor=max(max(abs(g1aa)));\n";
+ mDynamicModelFile << sp << " ze_elem=sum(diag(g1aa)<options_.cutoff);\n";
+ mDynamicModelFile << sp << " disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(options_.cutoff,'%f') ')']);\n";
+ mDynamicModelFile << sp << " if(correcting_factor<max_factor)\n";
+ mDynamicModelFile << sp << " correcting_factor=correcting_factor*4;\n";
+ mDynamicModelFile << sp << " disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.']);\n";
+ mDynamicModelFile << sp << " disp([' trying to correct the Jacobian matrix:']);\n";
+ mDynamicModelFile << sp << " disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]);\n";
+ mDynamicModelFile << sp << " dx = (g1aa+correcting_factor*speye(periods*Blck_size))\\ba- ya;\n";
+ mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';\n";
+ if (!block_triangular.ModelBlock->Block_List[i].is_linear)
+ mDynamicModelFile << sp << " continue;\n";
+ mDynamicModelFile << sp << " else\n";
+ mDynamicModelFile << sp << " disp('The singularity of the jacobian matrix could not be corrected');\n";
+ mDynamicModelFile << sp << " return;\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " elseif(lambda>1e-8)\n";
+ mDynamicModelFile << sp << " lambda=lambda/2;\n";
+ mDynamicModelFile << sp << " reduced = 1;\n";
+ mDynamicModelFile << sp << " disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);\n";
+ mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';\n";
+ if (!block_triangular.ModelBlock->Block_List[i].is_linear)
+ mDynamicModelFile << sp << " continue;\n";
+ mDynamicModelFile << sp << " else\n";
+ mDynamicModelFile << sp << " disp(['No convergence after ' num2str(iter,'%d') ' iterations']);\n";
+ mDynamicModelFile << sp << " return;\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " else\n";
+ mDynamicModelFile << sp << " if(lambda<1)\n";
+ mDynamicModelFile << sp << " lambda=max(lambda*2, 1);\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " end;\n";
+
+ mDynamicModelFile << sp << " ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)';\n";
+ mDynamicModelFile << sp << " ya_save=ya;\n";
+ mDynamicModelFile << sp << " g1aa=g1a;\n";
+ mDynamicModelFile << sp << " ba=b;\n";
+ mDynamicModelFile << sp << " max_resa=max_res;\n";
+ mDynamicModelFile << sp << " if(options_.simulation_method==0),\n";
+ mDynamicModelFile << sp << " dx = g1a\\b- ya;\n";
+ mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
+ mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
+ mDynamicModelFile << sp << " elseif(options_.simulation_method==2),\n";
+ mDynamicModelFile << sp << " flag1=1;\n";
+ mDynamicModelFile << sp << " while(flag1>0)\n";
+ mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n";
+ mDynamicModelFile << sp << " [za,flag1] = gmres(g1a,b," << block_triangular.ModelBlock->Block_List[i].Size << ",1e-6," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n";
+ mDynamicModelFile << sp << " if (flag1>0 | reduced)\n";
+ mDynamicModelFile << sp << " if(flag1==1)\n";
+ mDynamicModelFile << sp << " disp(['No convergence inside GMRES after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n";
+ mDynamicModelFile << sp << " elseif(flag1==2)\n";
+ mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n";
+ mDynamicModelFile << sp << " elseif(flag1==3)\n";
+ mDynamicModelFile << sp << " disp(['GMRES stagnated. (Two consecutive iterates were the same.)']);\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " luinc_tol = luinc_tol/10;\n";
+ mDynamicModelFile << sp << " reduced = 0;\n";
+ mDynamicModelFile << sp << " else\n";
+ mDynamicModelFile << sp << " dx = za - ya;\n";
+ mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
+ mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " elseif(options_.simulation_method==3),\n";
+ mDynamicModelFile << sp << " flag1=1;\n";
+ mDynamicModelFile << sp << " while(flag1>0)\n";
+ mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n";
+ mDynamicModelFile << sp << " [za,flag1] = bicgstabh(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n";
+ mDynamicModelFile << sp << " if (flag1>0 | reduced)\n";
+ mDynamicModelFile << sp << " if(flag1==1)\n";
+ mDynamicModelFile << sp << " disp(['No convergence inside BICGSTAB after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n";
+ mDynamicModelFile << sp << " elseif(flag1==2)\n";
+ mDynamicModelFile << sp << " disp(['Preconditioner is ill-conditioned ']);\n";
+ mDynamicModelFile << sp << " elseif(flag1==3)\n";
+ mDynamicModelFile << sp << " disp(['BICGSTAB stagnated. (Two consecutive iterates were the same.)']);\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " luinc_tol = luinc_tol/10;\n";
+ mDynamicModelFile << sp << " reduced = 0;\n";
+ mDynamicModelFile << sp << " else\n";
+ mDynamicModelFile << sp << " dx = za - ya;\n";
+ mDynamicModelFile << sp << " ya = ya + lambda*dx;\n";
+ mDynamicModelFile << sp << " y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " end;\n";
+ mDynamicModelFile << sp << " end;\n";
+ if(!block_triangular.ModelBlock->Block_List[i].is_linear)
+ {
+ mDynamicModelFile << " cvg=(max_res<solve_tolf);\n";
+ mDynamicModelFile << " iter=iter+1;\n";
+ }
+ mDynamicModelFile << " disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]);\n";
+ if(!block_triangular.ModelBlock->Block_List[i].is_linear)
+ {
+ mDynamicModelFile << " end\n";
+ mDynamicModelFile << " if (iter>maxit_)\n";
+ mDynamicModelFile << " disp(['No convergence after ' num2str(iter,'%4d') ' iterations']);\n";
+ mDynamicModelFile << " return;\n";
+ mDynamicModelFile << " end;\n";
}
}
- if(mode==eSparseMode)
- writeModelEquationsOrdered_M(mDynamicModelFile, block_triangular.ModelBlock, dynamic_basename);
- mDynamicModelFile.close();
+
+ prev_Simulation_Type=k;
+ if(open_par)
+ mDynamicModelFile << " end;\n";
+ mDynamicModelFile << " oo_.endo_simul = y';\n";
+ mDynamicModelFile << "return;\n";
}
+
+ writeModelEquationsOrdered_M(mDynamicModelFile, block_triangular.ModelBlock, dynamic_basename);
+
+ mDynamicModelFile.close();
+
if (printed)
cout << "done\n";
}
@@ -4007,6 +3088,7 @@ ModelTree::writeStaticFile(const string &basename) const
break;
case eSparseMode:
writeSparseStaticMFile(basename + "_static", basename, mode);
+ break;
case eDLLMode:
writeStaticCFile(basename + "_static");
break;
@@ -4016,14 +3098,13 @@ ModelTree::writeStaticFile(const string &basename) const
void
ModelTree::writeDynamicFile(const string &basename) const
{
- ExprNodeOutputType output_type = (mode == eDLLMode ? oCStaticModel : oMatlabStaticModel);
switch(mode)
{
case eStandardMode:
writeDynamicMFile(basename + "_dynamic");
break;
case eSparseMode:
- writeSparseDynamicFileAndBinFile(basename + "_dynamic", basename, output_type, mode);
+ writeSparseDynamicMFile(basename + "_dynamic", basename, mode);
block_triangular.Free_Block(block_triangular.ModelBlock);
block_triangular.Free_IM(block_triangular.First_IM);
break;
@@ -4031,9 +3112,7 @@ ModelTree::writeDynamicFile(const string &basename) const
writeDynamicCFile(basename + "_dynamic");
break;
case eSparseDLLMode:
- writeSparseDynamicFileAndBinFile(basename + "_dynamic", basename, output_type, mode);
- if (compiler==GCC_COMPILE || compiler==LCC_COMPILE )
- writeSparseDLLDynamicHFile(basename + "_dynamic");
+ writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL);
block_triangular.Free_Block(block_triangular.ModelBlock);
block_triangular.Free_IM(block_triangular.First_IM);
break;
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index 52e885ec..cd2cd313 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -24,14 +24,6 @@
#include "ParsingDriver.hh"
#include "Statement.hh"
-ParsingDriver::ParsingDriver()
-{
-}
-
-ParsingDriver::~ParsingDriver()
-{
-}
-
bool
ParsingDriver::symbol_exists_and_is_not_modfile_local_variable(const char *s)
{
@@ -705,16 +697,10 @@ void ParsingDriver::simulate()
void
ParsingDriver::simul_sparse()
{
- mod_file->addStatement(new SimulSparseStatement(options_list, mod_file->model_tree.compiler, mod_file->model_tree.mode));
+ mod_file->addStatement(new SimulSparseStatement(options_list, mod_file->model_tree.mode));
options_list.clear();
}
-void
-ParsingDriver::init_compiler(int compiler_type)
-{
- mod_file->model_tree.compiler = compiler_type;
-}
-
void
ParsingDriver::simul()
{
diff --git a/include/ComputingTasks.hh b/include/ComputingTasks.hh
index 08d24a1d..96e547e2 100644
--- a/include/ComputingTasks.hh
+++ b/include/ComputingTasks.hh
@@ -69,9 +69,9 @@ class SimulSparseStatement : public Statement
{
private:
const OptionsList options_list;
- const int compiler, mode;
+ const int mode;
public:
- SimulSparseStatement(const OptionsList &options_list_arg, int compiler_arg, int mode_arg);
+ SimulSparseStatement(const OptionsList &options_list_arg, int mode_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
diff --git a/include/ModelTree.hh b/include/ModelTree.hh
index 1e4fbcea..7cb392e6 100644
--- a/include/ModelTree.hh
+++ b/include/ModelTree.hh
@@ -33,11 +33,6 @@ using namespace std;
#include "BlockTriangular.hh"
#include "SymbolGaussElim.hh"
-#define LCC_COMPILE 0
-#define GCC_COMPILE 1
-#define NO_COMPILE 2
-//#define CONDITION
-
//! The three in which ModelTree can work
enum ModelTreeMode
{
@@ -106,13 +101,11 @@ private:
//! Writes the dynamic model equations and its derivatives
/*! \todo add third derivatives handling in C output */
void writeDynamicModel(ostream &DynamicOutput) const;
- //! Writes the Block reordred structure of the model in C output
- void writeModelEquationsOrdered_C(ostream &output, Model_Block *ModelBlock) const;
//! Writes the Block reordred structure of the model in M output
void writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock, const string &dynamic_basename) const;
//! Writes the Block reordred structure of the static model in M output
void writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *ModelBlock, const string &static_basename) const;
- //! Writes the code of the Block reordred structure of the model in C output
+ //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, ExprNodeOutputType output_type) const;
//! Writes static model file (Matlab version)
void writeStaticMFile(const string &static_basename) const;
@@ -125,10 +118,8 @@ private:
//! Writes dynamic model file (C version)
/*! \todo add third derivatives handling */
void writeDynamicCFile(const string &dynamic_basename) const;
- //! Writes dynamic model header file when SparseDLL option is on
- void writeSparseDLLDynamicHFile(const string &dynamic_basename) const;
//! Writes dynamic model file when SparseDLL option is on
- void writeSparseDynamicFileAndBinFile(const string &dynamic_basename, const string &bin_basename, ExprNodeOutputType output_type, const int mode) const;
+ void writeSparseDynamicMFile(const string &dynamic_basename, const string &basename, const int mode) const;
//! Computes jacobian and prepares for equation normalization
/*! Using values from initval/endval blocks and parameter initializations:
- computes the jacobian for the model w.r. to contemporaneous variables
@@ -145,8 +136,6 @@ public:
ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
//! Mode in which the ModelTree is supposed to work (Matlab, DLL or SparseDLL)
ModelTreeMode mode;
- //! Type of compiler used in matlab for SPARSE_DLL option: 0 = LCC or 1 = GCC or 2 = NO
- int compiler;
//! Absolute value under which a number is considered to be zero
double cutoff;
//! The weight of the Markowitz criteria to determine the pivot in the linear solver (simul_NG1 from simulate.cc)
@@ -181,7 +170,7 @@ public:
//! Adds informations for simulation in a binary file
void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename,
const int &num, int &u_count_int, bool &file_open) const;
-
+ //! Returns the number of equations in the model
int equation_number() const;
};
diff --git a/include/ParsingDriver.hh b/include/ParsingDriver.hh
index 3ce2065f..c9f60996 100644
--- a/include/ParsingDriver.hh
+++ b/include/ParsingDriver.hh
@@ -139,11 +139,6 @@ private:
ModFile *mod_file;
public:
- //! Constructor
- ParsingDriver();
- //! Destructor
- virtual ~ParsingDriver();
-
//! Starts parsing, and constructs the MOD file representation
/*! The returned pointer should be deleted after use */
ModFile *parse(istream &in, bool debug);
@@ -172,8 +167,6 @@ public:
void sparse_dll();
//! Sets mode of ModelTree class to block decompose the model and triggers the creation of the incidence matrix in Matlab context
void sparse();
- //! Sets the compiler type used in conjunction with SPARCE_DLL
- void init_compiler(int compiler_type);
//! Sets the FILENAME for the initial value in initval
void initval_file(string *filename);
//! Declares an endogenous variable
--
GitLab