stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Adding Markowitz criteria in the linear solver (new option in simul: "Markowitz=val" - with val a strictly positive real) 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1298 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2007-06-03 22:35:30 +00:00
parent 500727c0cb
commit 77d6a24cf7
22 changed files with 2595 additions and 2232 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
matlab/dynare_m.exe
View File

Binary file not shown.

49
matlab/global_initialization.m
View File

@ -1,8 +1,8 @@
function global_initialization() function global_initialization()
% initializes global variables and options for DYNARE % initializes global variables and options for DYNARE
global oo_ M_ options_ ct_ endval_ rplottype_ global oo_ M_ options_ ct_ endval_ rplottype_
ct_=0; ct_=0;
endval_=0; endval_=0;
@ -21,8 +21,8 @@ function global_initialization()
options_.solve_tolf = eps^(2/3); options_.solve_tolf = eps^(2/3);
options_.solve_tolx = 3.7e-11; options_.solve_tolx = 3.7e-11;
options_.solve_maxit = 500; options_.solve_maxit = 500;
% steady state file % steady state file
if exist([M_.fname '_steadystate']) if exist([M_.fname '_steadystate'])
options_.steadystate_flag = 1; options_.steadystate_flag = 1;
@ -30,17 +30,17 @@ function global_initialization()
options_.steadystate_flag = 0; options_.steadystate_flag = 0;
end end
options_.steadystate_partial = []; options_.steadystate_partial = [];
% subset of the estimated deep parameters % subset of the estimated deep parameters
options_.ParamSubSet = 'None'; options_.ParamSubSet = 'None';
% bvar-dsge % bvar-dsge
options_.varlag = 4; options_.varlag = 4;
% Optimization algorithm [6] gmhmaxlik % Optimization algorithm [6] gmhmaxlik
options_.Opt6Iter = 3; options_.Opt6Iter = 3;
options_.Opt6Numb = 100000; options_.Opt6Numb = 100000;
% Graphics % Graphics
options_.graphics.nrows = 3; options_.graphics.nrows = 3;
options_.graphics.ncols = 3; options_.graphics.ncols = 3;
@ -48,8 +48,8 @@ function global_initialization()
options_.graphics.line_width = 1; options_.graphics.line_width = 1;
options_.nograph = 0; options_.nograph = 0;
options_.XTick = []; options_.XTick = [];
options_.XTickLabel = []; options_.XTickLabel = [];
% IRFs & other stoch_simul output % IRFs & other stoch_simul output
options_.irf = 40; options_.irf = 40;
options_.relative_irf = 0; options_.relative_irf = 0;
@ -63,21 +63,21 @@ function global_initialization()
options_.noprint = 0; options_.noprint = 0;
options_.simul = 0; options_.simul = 0;
options_.SpectralDensity = 0; options_.SpectralDensity = 0;
% TeX output % TeX output
options_.TeX = 0; options_.TeX = 0;
% Exel % Exel
options_.xls_sheet = ''; options_.xls_sheet = '';
options_.xls_range = ''; options_.xls_range = '';
% Prior draws % Prior draws
options_.forecast = 0; options_.forecast = 0;
options_.replic = 1; options_.replic = 1;
% Model % Model
options_.linear = 0; options_.linear = 0;
% Solution % Solution
options_.order = 2; options_.order = 2;
options_.dr_algo = 0; options_.dr_algo = 0;
@ -86,11 +86,11 @@ function global_initialization()
options_.replic = 50; options_.replic = 50;
options_.drop = 100; options_.drop = 100;
options_.simul_algo = 0; options_.simul_algo = 0;
% Ramsey policy % Ramsey policy
options_.planner_discount = 1.0; options_.planner_discount = 1.0;
options_.ramsey_policy = 0; options_.ramsey_policy = 0;
% estimation % estimation
options_.load_mh_file = 0; options_.load_mh_file = 0;
options_.first_obs = 1; options_.first_obs = 1;
@ -107,14 +107,14 @@ function global_initialization()
options_.mode_check = 0; options_.mode_check = 0;
options_.prior_trunc = 1e-10; options_.prior_trunc = 1e-10;
options_.mh_conf_sig = 0.90; options_.mh_conf_sig = 0.90;
options_.mh_mode = 1; options_.mh_mode = 1;
options_.mh_nblck = 2; options_.mh_nblck = 2;
options_.load_mh_file = 0; options_.load_mh_file = 0;
options_.nodiagnostic = 0; options_.nodiagnostic = 0;
options_.loglinear = 0; options_.loglinear = 0;
options_.unit_root_vars = []; options_.unit_root_vars = [];
options_.bayesian_irf = 0; options_.bayesian_irf = 0;
options_.bayesian_th_moments = 0; options_.bayesian_th_moments = 0;
options_.smoother = 0; options_.smoother = 0;
options_.moments_varendo = 0; options_.moments_varendo = 0;
options_.filtered_vars = 0; options_.filtered_vars = 0;
@ -126,8 +126,9 @@ function global_initialization()
options_.diffuse_d = []; options_.diffuse_d = [];
options_.logdata = 0; options_.logdata = 0;
options_.use_mh_covariance_matrix = 0; options_.use_mh_covariance_matrix = 0;
options_.noconstant = 0; options_.noconstant = 0;
options_.markowitz = 0.5;
% Misc % Misc
options_.conf_sig = 0.9; options_.conf_sig = 0.9;
oo_.exo_simul = []; oo_.exo_simul = [];
@ -136,9 +137,9 @@ function global_initialization()
oo_.exo_steady_state = []; oo_.exo_steady_state = [];
oo_.exo_det_steady_state = []; oo_.exo_det_steady_state = [];
oo_.exo_det_simul = []; oo_.exo_det_simul = [];
% Variance matrix for measurement errors % Variance matrix for measurement errors
M_.H = 0; M_.H = 0;
% BVAR % BVAR
M_.bvar = []; M_.bvar = [];

43
matlab/writedata_text.m
View File

@ -1,8 +1,8 @@
function writedata(fname) function writedata(fname)
% function writedata(fname) % function writedata(fname)
% store endogenous and exogenous variables in a XLS spreadsheet file % store endogenous and exogenous variables in a text file
% INPUT % INPUT
% fname: name of the XLS file % fname: name of the text file
% OUTPUT % OUTPUT
% none % none
% ALGORITHM % ALGORITHM
@ -13,20 +13,27 @@ function writedata(fname)
% part of DYNARE, copyright (2007) % part of DYNARE, copyright (2007)
% Gnu Public License. % Gnu Public License.
global M_ oo_ global M_ oo_
S=[fname '_endo.xls']; S=[fname '_endo.dat'];
n=size(oo_.endo_simul,2); fid = fopen(S,'w');
delete(S); for i = 1:size(M_.endo_names,1)
S=upper(cellstr(M_.endo_names)); fprintf(fid,'%s ',M_.endo_names(i,:)');
S1=cellstr([num2str((1:n)') char(65*ones(1,n))']); end;
xlswrite([fname '_endo'], S', 'endogenous', 'B1'); fprintf(fid,'\n');
xlswrite([fname '_endo'], S1, 'endogenous', 'A2'); for i = 1:size(oo_.endo_simul,2)
xlswrite([fname '_endo'], oo_.endo_simul', 'endogenous', 'B2'); fprintf(fid,'%15.7f ',oo_.endo_simul(:,i));
S=[fname '_exo.xls']; fprintf(fid,'\n');
n=size(oo_.exo_simul,1); end
delete(S); fclose(fid);
S=upper(cellstr(M_.exo_names));
S1=cellstr([num2str((1:n)') char(65*ones(1,n))']); S=[fname '_exo.dat'];
xlswrite([fname '_exo'], S','exogenous', 'B1'); fid = fopen(S,'w');
xlswrite([fname '_exo'], S1, 'exogenous', 'A2'); for i = 1:size(M_.exo_names,1)
xlswrite([fname '_exo'], oo_.exo_simul,'exogenous', 'B2'); fprintf(fid,'%s ',M_.exo_names(i,:));
end;
fprintf(fid,'\n');
for i = 1:size(oo_.exo_simul,1)
fprintf(fid,'%15.7f ',oo_.exo_simul(i,:));
fprintf(fid,'\n');
end
fclose(fid);
return; return;

4
parser.src/BlockTriangular.cc
View File

@ -777,7 +777,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
SIM0[iter->first.first*n+iter->first.second]=1; SIM0[iter->first.first*n+iter->first.second]=1;
if(!IM_0[iter->first.first*n+iter->first.second]) if(!IM_0[iter->first.first*n+iter->first.second])
{ {
cout << "error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << "\n"; cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << "\n";
} }
} }
else else
@ -791,7 +791,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM); OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM);
suppressed=suppress; suppressed=suppress;
if(!OK) if(!OK)
bi/=1.05; bi/=1.07;
if(bi>1e-14) if(bi>1e-14)
free(SIM00); free(SIM00);
} }

11
parser.src/ComputingTasks.cc
View File

@ -252,6 +252,17 @@ CutoffStatement::writeOutput(ostream &output, const string &basename) const
output << "options_.cutoff = " << cutoff << ";" << endl; output << "options_.cutoff = " << cutoff << ";" << endl;
} }
MarkowitzStatement::MarkowitzStatement(double markowitz_arg) : markowitz(markowitz_arg)
{
}
void
MarkowitzStatement::writeOutput(ostream &output, const string &basename) const
{
output << "options_.markowitz = " << markowitz << ";" << endl;
}
DsampleStatement::DsampleStatement(int val1_arg) : val1(val1_arg), val2(-1) DsampleStatement::DsampleStatement(int val1_arg) : val1(val1_arg), val2(-1)
{ {
} }

3553
parser.src/DynareBison.cc
View File

File diff suppressed because it is too large Load Diff

17
parser.src/DynareBison.yy
View File

@ -47,7 +47,7 @@ class ParsingDriver;
%token <string_val> INT_NUMBER %token <string_val> INT_NUMBER
%token INV_GAMMA_PDF IRF %token INV_GAMMA_PDF IRF
%token KALMAN_ALGO KALMAN_TOL %token KALMAN_ALGO KALMAN_TOL
%token LAPLACE LCC_COMPILER LIK_ALGO LIK_INIT LINEAR LOAD_MH_FILE LOGLINEAR %token LAPLACE LCC_COMPILER LIK_ALGO LIK_INIT LINEAR LOAD_MH_FILE LOGLINEAR MARKOWITZ
%token MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER %token MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER
%token MODE_CHECK MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MSHOCKS %token MODE_CHECK MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MSHOCKS
%token MODEL_COMPARISON_APPROXIMATION MODIFIEDHARMONICMEAN MOMENTS_VARENDO %token MODEL_COMPARISON_APPROXIMATION MODIFIEDHARMONICMEAN MOMENTS_VARENDO
@ -89,6 +89,7 @@ class ParsingDriver;
: declaration : declaration
| periods | periods
| cutoff | cutoff
| markowitz
| model | model
| initval | initval
| endval | endval
@ -239,6 +240,18 @@ cutoff
} }
; ;
markowitz
: MARKOWITZ FLOAT_NUMBER ';'
{
driver.markowitz($2);
}
| MARKOWITZ EQUAL FLOAT_NUMBER ';'
{
driver.markowitz($3);
}
;
init_param init_param
: NAME EQUAL expression ';' : NAME EQUAL expression ';'
{driver.init_param($1, $3);} {driver.init_param($1, $3);}
@ -346,6 +359,7 @@ cutoff
LCC_COMPILER { driver.init_compiler(0); } LCC_COMPILER { driver.init_compiler(0); }
| GCC_COMPILER { driver.init_compiler(1); } | GCC_COMPILER { driver.init_compiler(1); }
| o_cutoff | o_cutoff
| o_markowitz
; ;
model model
@ -1118,6 +1132,7 @@ cutoff
o_hp_ngrid: HP_NGRID EQUAL INT_NUMBER {driver.option_num("hp_ngrid", $3);}; o_hp_ngrid: HP_NGRID EQUAL INT_NUMBER {driver.option_num("hp_ngrid", $3);};
o_periods: PERIODS EQUAL INT_NUMBER {driver.option_num("periods", $3); driver.option_num("simul", "1");}; o_periods: PERIODS EQUAL INT_NUMBER {driver.option_num("periods", $3); driver.option_num("simul", "1");};
o_cutoff: CUTOFF EQUAL FLOAT_NUMBER {driver.option_num("cutoff", $3);} o_cutoff: CUTOFF EQUAL FLOAT_NUMBER {driver.option_num("cutoff", $3);}
o_markowitz: MARKOWITZ EQUAL FLOAT_NUMBER {driver.option_num("markowitz", $3);}
o_simul: SIMUL {driver.option_num("simul", "1");}; o_simul: SIMUL {driver.option_num("simul", "1");};
o_simul_seed: SIMUL_SEED EQUAL INT_NUMBER { driver.option_num("simul_seed", $3)}; o_simul_seed: SIMUL_SEED EQUAL INT_NUMBER { driver.option_num("simul_seed", $3)};
o_qz_criterium: QZ_CRITERIUM EQUAL INT_NUMBER { driver.option_num("qz_criterium", $3)} o_qz_criterium: QZ_CRITERIUM EQUAL INT_NUMBER { driver.option_num("qz_criterium", $3)}

3
parser.src/DynareFlex.ll
View File

@ -58,6 +58,7 @@ int sigma_e = 0;
<INITIAL>parameters {BEGIN DYNARE_STATEMENT; return token::PARAMETERS;} <INITIAL>parameters {BEGIN DYNARE_STATEMENT; return token::PARAMETERS;}
<INITIAL>periods {BEGIN DYNARE_STATEMENT; return token::PERIODS;} <INITIAL>periods {BEGIN DYNARE_STATEMENT; return token::PERIODS;}
<INITIAL>cutoff {BEGIN DYNARE_STATEMENT; return token::CUTOFF;} <INITIAL>cutoff {BEGIN DYNARE_STATEMENT; return token::CUTOFF;}
<INITIAL>markowitz {BEGIN DYNARE_STATEMENT; return token::MARKOWITZ;}
<INITIAL>estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;} <INITIAL>estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;}
<INITIAL>prior_analysis {BEGIN DYNARE_STATEMENT; return token::PRIOR_ANALYSIS;} <INITIAL>prior_analysis {BEGIN DYNARE_STATEMENT; return token::PRIOR_ANALYSIS;}
<INITIAL>posterior_analysis {BEGIN DYNARE_STATEMENT; return token::POSTERIOR_ANALYSIS;} <INITIAL>posterior_analysis {BEGIN DYNARE_STATEMENT; return token::POSTERIOR_ANALYSIS;}
@ -150,6 +151,7 @@ int sigma_e = 0;
<DYNARE_STATEMENT>optim {return token::OPTIM;} <DYNARE_STATEMENT>optim {return token::OPTIM;}
<DYNARE_STATEMENT>periods {return token::PERIODS;} <DYNARE_STATEMENT>periods {return token::PERIODS;}
<DYNARE_STATEMENT>cutoff {return token::CUTOFF;} <DYNARE_STATEMENT>cutoff {return token::CUTOFF;}
<DYNARE_STATEMENT>markowitz {return token::MARKOWITZ;}
<DYNARE_STATEMENT>model_comparison_approximation {return token::MODEL_COMPARISON;} <DYNARE_STATEMENT>model_comparison_approximation {return token::MODEL_COMPARISON;}
<DYNARE_STATEMENT>laplace {return token::LAPLACE;} <DYNARE_STATEMENT>laplace {return token::LAPLACE;}
<DYNARE_STATEMENT>modifiedharmonicmean {return token::MODIFIEDHARMONICMEAN;} <DYNARE_STATEMENT>modifiedharmonicmean {return token::MODIFIEDHARMONICMEAN;}
@ -170,6 +172,7 @@ int sigma_e = 0;
<DYNARE_BLOCK>corr {return token::CORR;} <DYNARE_BLOCK>corr {return token::CORR;}
<DYNARE_BLOCK>periods {return token::PERIODS;} <DYNARE_BLOCK>periods {return token::PERIODS;}
<DYNARE_BLOCK>cutoff {return token::CUTOFF;} <DYNARE_BLOCK>cutoff {return token::CUTOFF;}
<DYNARE_BLOCK>markowitz {return token::MARKOWITZ;}
<DYNARE_BLOCK>filename {return token::FILENAME;} <DYNARE_BLOCK>filename {return token::FILENAME;}
<DYNARE_BLOCK>gamma_pdf {return token::GAMMA_PDF;} <DYNARE_BLOCK>gamma_pdf {return token::GAMMA_PDF;}
<DYNARE_BLOCK>beta_pdf {return token::BETA_PDF;} <DYNARE_BLOCK>beta_pdf {return token::BETA_PDF;}

9
parser.src/ExprNode.cc
View File

@ -337,13 +337,14 @@ VariableNode::eval(const eval_context_type &eval_context) const throw (EvalExcep
{ {
/*if (lag != 0) /*if (lag != 0)
throw EvalException();*/ throw EvalException();*/
if(&eval_context==NULL)
throw EvalException();
eval_context_type::const_iterator it = eval_context.find(make_pair(symb_id, type)); eval_context_type::const_iterator it = eval_context.find(make_pair(symb_id, type));
if (it == eval_context.end()) if (it == eval_context.end())
{ {
cout << "Error: the variable or parameter (" << datatree.symbol_table.getNameByID( type, symb_id) << ") has not been initialized (in derivatives evaluation)\n"; if (eval_context.size()>0)
cout.flush(); {
cout << "Error: the variable or parameter (" << datatree.symbol_table.getNameByID( type, symb_id) << ") has not been initialized (in derivatives evaluation)\n";
cout.flush();
}
throw EvalException(); throw EvalException();
} }
return it->second; return it->second;

14
parser.src/ModelTree.cc
View File

@ -15,6 +15,7 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
mode(eStandardMode), mode(eStandardMode),
compiler(LCC_COMPILE), compiler(LCC_COMPILE),
cutoff(1e-12), cutoff(1e-12),
markowitz(0.7),
new_SGE(true), new_SGE(true),
computeJacobian(false), computeJacobian(false),
computeJacobianExo(false), computeJacobianExo(false),
@ -1473,7 +1474,7 @@ ModelTree::writeSparseDLLDynamicCFileAndBinFile(const string &dynamic_basename,
mDynamicModelFile << " Read_SparseMatrix(\"" << reform(bin_basename) << "\"," mDynamicModelFile << " Read_SparseMatrix(\"" << reform(bin_basename) << "\","
<< block_triangular.ModelBlock->Block_List[i].Size << ", periods, y_kmin, y_kmax" << block_triangular.ModelBlock->Block_List[i].Size << ", periods, y_kmin, y_kmax"
<< ");\n"; << ");\n";
mDynamicModelFile << " u_count=" << u_count_int << "*(periods+y_kmax);\n"; mDynamicModelFile << " u_count=" << u_count_int << "*(periods+y_kmax+y_kmin);\n";
} }
else else
{ {
@ -1509,11 +1510,10 @@ ModelTree::writeSparseDLLDynamicCFileAndBinFile(const string &dynamic_basename,
mDynamicModelFile << " }\n"; mDynamicModelFile << " }\n";
mDynamicModelFile << " }\n"; mDynamicModelFile << " }\n";
mDynamicModelFile << " cvg=(max_res<solve_tolf);\n"; mDynamicModelFile << " cvg=(max_res<solve_tolf);\n";
mDynamicModelFile << " if(!cvg)\n";
if(new_SGE) 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);\n"; 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 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 << " 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 << " iter++;\n";
mDynamicModelFile << " }\n"; mDynamicModelFile << " }\n";
mDynamicModelFile << " if (!cvg)\n"; mDynamicModelFile << " if (!cvg)\n";
@ -1538,7 +1538,7 @@ ModelTree::writeSparseDLLDynamicCFileAndBinFile(const string &dynamic_basename,
mDynamicModelFile << "#endif\n"; mDynamicModelFile << "#endif\n";
mDynamicModelFile << " }\n"; mDynamicModelFile << " }\n";
if(new_SGE) 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);\n"; 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 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 << " 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";
} }
@ -1617,6 +1617,7 @@ ModelTree::writeSparseDLLDynamicCFileAndBinFile(const string &dynamic_basename,
mDynamicModelFile << " periods=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"periods\"))))));\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 << " maxit_=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"maxit_\"))))));\n";
mDynamicModelFile << " slowc=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"slowc\")))));\n"; mDynamicModelFile << " slowc=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"slowc\")))));\n";
mDynamicModelFile << " markowitz_c=double(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"markowitz\")))));\n";
} }
else else
{ {
@ -1624,7 +1625,8 @@ ModelTree::writeSparseDLLDynamicCFileAndBinFile(const string &dynamic_basename,
mDynamicModelFile << " y_kmax=(int)floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"maximum_lead\")))));\n"; mDynamicModelFile << " y_kmax=(int)floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"maximum_lead\")))));\n";
mDynamicModelFile << " periods=(int)floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"periods\")))));\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 << " maxit_=(int)floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"maxit_\")))));\n";
mDynamicModelFile << " slowc=(int)floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,\"slowc\")))));\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 << " col_y=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,\"endo_simul\")));;\n"; mDynamicModelFile << " col_y=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,\"endo_simul\")));;\n";
mDynamicModelFile << " if (col_y<row_x)\n"; mDynamicModelFile << " if (col_y<row_x)\n";

2
parser.src/NumericalConstants.cc
View File

@ -12,7 +12,7 @@ int
NumericalConstants::AddConstant(const string &iConst) NumericalConstants::AddConstant(const string &iConst)
{ {
map<string, int>::iterator iter = numConstantsIndex.find(iConst); map<string, int>::iterator iter = numConstantsIndex.find(iConst);
//cout << "iConst=" << iConst << "\n" ;
if (iter != numConstantsIndex.end()) if (iter != numConstantsIndex.end())
return iter->second; return iter->second;

9
parser.src/ParsingDriver.cc
View File

@ -192,6 +192,15 @@ ParsingDriver::cutoff(string *cutoff)
delete cutoff; delete cutoff;
} }
void
ParsingDriver::markowitz(string *markowitz)
{
double markowitz_val = atof(markowitz->c_str());
mod_file->addStatement(new MarkowitzStatement(markowitz_val));
delete markowitz;
}
void void
ParsingDriver::dsample(string *arg1) ParsingDriver::dsample(string *arg1)
{ {

10
parser.src/include/ComputingTasks.hh
View File

@ -110,6 +110,16 @@ public:
virtual void writeOutput(ostream &output, const string &basename) const; virtual void writeOutput(ostream &output, const string &basename) const;
}; };
class MarkowitzStatement : public Statement
{
private:
const double markowitz;
public:
MarkowitzStatement(double markowitz_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class DsampleStatement : public Statement class DsampleStatement : public Statement
{ {
private: private:

8
parser.src/include/DataTree.hh
View File

@ -152,8 +152,7 @@ DataTree::AddUnaryOp(UnaryOpcode op_code, NodeID arg)
{ {
try try
{ {
eval_context_type *evc=NULL; double argval = arg->eval(eval_context_type());
double argval = arg->eval(/*eval_context_type()*/*evc);
double val = UnaryOpNode::eval_opcode(op_code, argval); double val = UnaryOpNode::eval_opcode(op_code, argval);
return AddPossiblyNegativeConstant(val); return AddPossiblyNegativeConstant(val);
} }
@ -174,9 +173,8 @@ DataTree::AddBinaryOp(NodeID arg1, BinaryOpcode op_code, NodeID arg2)
// Try to reduce to a constant // Try to reduce to a constant
try try
{ {
eval_context_type *evc=NULL; double argval1 = arg1->eval(eval_context_type());
double argval1 = arg1->eval(/*eval_context_type()*/*evc); double argval2 = arg2->eval(eval_context_type());
double argval2 = arg2->eval(/*eval_context_type()*/*evc);
double val = BinaryOpNode::eval_opcode(argval1, op_code, argval2); double val = BinaryOpNode::eval_opcode(argval1, op_code, argval2);
return AddPossiblyNegativeConstant(val); return AddPossiblyNegativeConstant(val);
} }

195
parser.src/include/DynareBison.hh
View File

@ -177,102 +177,103 @@ namespace yy
LINEAR = 305, LINEAR = 305,
LOAD_MH_FILE = 306, LOAD_MH_FILE = 306,
LOGLINEAR = 307, LOGLINEAR = 307,
MH_DROP = 308, MARKOWITZ = 308,
MH_INIT_SCALE = 309, MH_DROP = 309,
MH_JSCALE = 310, MH_INIT_SCALE = 310,
MH_MODE = 311, MH_JSCALE = 311,
MH_NBLOCKS = 312, MH_MODE = 312,
MH_REPLIC = 313, MH_NBLOCKS = 313,
MH_RECOVER = 314, MH_REPLIC = 314,
MODE_CHECK = 315, MH_RECOVER = 315,
MODE_COMPUTE = 316, MODE_CHECK = 316,
MODE_FILE = 317, MODE_COMPUTE = 317,
MODEL = 318, MODE_FILE = 318,
MODEL_COMPARISON = 319, MODEL = 319,
MSHOCKS = 320, MODEL_COMPARISON = 320,
MODEL_COMPARISON_APPROXIMATION = 321, MSHOCKS = 321,
MODIFIEDHARMONICMEAN = 322, MODEL_COMPARISON_APPROXIMATION = 322,
MOMENTS_VARENDO = 323, MODIFIEDHARMONICMEAN = 323,
NAME = 324, MOMENTS_VARENDO = 324,
NOBS = 325, NAME = 325,
NOCONSTANT = 326, NOBS = 326,
NOCORR = 327, NOCONSTANT = 327,
NODIAGNOSTIC = 328, NOCORR = 328,
NOFUNCTIONS = 329, NODIAGNOSTIC = 329,
NOGRAPH = 330, NOFUNCTIONS = 330,
NOMOMENTS = 331, NOGRAPH = 331,
NOPRINT = 332, NOMOMENTS = 332,
NORMAL_PDF = 333, NOPRINT = 333,
OBSERVATION_TRENDS = 334, NORMAL_PDF = 334,
OLR = 335, OBSERVATION_TRENDS = 335,
OLR_INST = 336, OLR = 336,
OLR_BETA = 337, OLR_INST = 337,
OPTIM = 338, OLR_BETA = 338,
OPTIM_WEIGHTS = 339, OPTIM = 339,
ORDER = 340, OPTIM_WEIGHTS = 340,
OSR = 341, ORDER = 341,
OSR_PARAMS = 342, OSR = 342,
PARAMETERS = 343, OSR_PARAMS = 343,
PERIODS = 344, PARAMETERS = 344,
PLANNER_OBJECTIVE = 345, PERIODS = 345,
PREFILTER = 346, PLANNER_OBJECTIVE = 346,
PRESAMPLE = 347, PREFILTER = 347,
PRINT = 348, PRESAMPLE = 348,
PRIOR_TRUNC = 349, PRINT = 349,
PRIOR_ANALYSIS = 350, PRIOR_TRUNC = 350,
POSTERIOR_ANALYSIS = 351, PRIOR_ANALYSIS = 351,
QZ_CRITERIUM = 352, POSTERIOR_ANALYSIS = 352,
RELATIVE_IRF = 353, QZ_CRITERIUM = 353,
REPLIC = 354, RELATIVE_IRF = 354,
RPLOT = 355, REPLIC = 355,
SHOCKS = 356, RPLOT = 356,
SIGMA_E = 357, SHOCKS = 357,
SIMUL = 358, SIGMA_E = 358,
SIMUL_ALGO = 359, SIMUL = 359,
SIMUL_SEED = 360, SIMUL_ALGO = 360,
SMOOTHER = 361, SIMUL_SEED = 361,
SOLVE_ALGO = 362, SMOOTHER = 362,
SPARSE_DLL = 363, SOLVE_ALGO = 363,
STDERR = 364, SPARSE_DLL = 364,
STEADY = 365, STDERR = 365,
STOCH_SIMUL = 366, STEADY = 366,
TEX = 367, STOCH_SIMUL = 367,
RAMSEY_POLICY = 368, TEX = 368,
PLANNER_DISCOUNT = 369, RAMSEY_POLICY = 369,
TEX_NAME = 370, PLANNER_DISCOUNT = 370,
UNIFORM_PDF = 371, TEX_NAME = 371,
UNIT_ROOT_VARS = 372, UNIFORM_PDF = 372,
USE_DLL = 373, UNIT_ROOT_VARS = 373,
VALUES = 374, USE_DLL = 374,
VAR = 375, VALUES = 375,
VAREXO = 376, VAR = 376,
VAREXO_DET = 377, VAREXO = 377,
VAROBS = 378, VAREXO_DET = 378,
XLS_SHEET = 379, VAROBS = 379,
XLS_RANGE = 380, XLS_SHEET = 380,
COMMA = 381, XLS_RANGE = 381,
MINUS = 382, COMMA = 382,
PLUS = 383, MINUS = 383,
DIVIDE = 384, PLUS = 384,
TIMES = 385, DIVIDE = 385,
UMINUS = 386, TIMES = 386,
POWER = 387, UMINUS = 387,
EXP = 388, POWER = 388,
LOG = 389, EXP = 389,
LOG10 = 390, LOG = 390,
SIN = 391, LOG10 = 391,
COS = 392, SIN = 392,
TAN = 393, COS = 393,
ASIN = 394, TAN = 394,
ACOS = 395, ASIN = 395,
ATAN = 396, ACOS = 396,
SINH = 397, ATAN = 397,
COSH = 398, SINH = 398,
TANH = 399, COSH = 399,
ASINH = 400, TANH = 400,
ACOSH = 401, ASINH = 401,
ATANH = 402, ACOSH = 402,
SQRT = 403 ATANH = 403,
SQRT = 404
}; };
}; };
@ -367,7 +368,7 @@ namespace yy
static const short int yytable_[]; static const short int yytable_[];
static const signed char yytable_ninf_; static const signed char yytable_ninf_;
static const unsigned short int yycheck_[]; static const short int yycheck_[];
/// For a state, its accessing symbol. /// For a state, its accessing symbol.
static const unsigned short int yystos_[]; static const unsigned short int yystos_[];

2
parser.src/include/ModelTree.hh
View File

@ -108,6 +108,8 @@ public:
int compiler; int compiler;
//! Absolute value under which a number is considered to be zero //! Absolute value under which a number is considered to be zero
double cutoff; double cutoff;
//! The weight of the Markowitz criteria to determine the pivot in the linear solver (simul_NG1 from simulate.cc)
double markowitz;
//! Use a graphical and symbolic version of the symbolic gaussian elimination new_SGE = false //! Use a graphical and symbolic version of the symbolic gaussian elimination new_SGE = false
//! or use direct gaussian elimination new_SGE = true //! or use direct gaussian elimination new_SGE = true
bool new_SGE; bool new_SGE;

6
parser.src/include/ParsingDriver.hh
View File

@ -179,6 +179,8 @@ public:
void periods(string *periods); void periods(string *periods);
//! Adds a "cutoff" statement //! Adds a "cutoff" statement
void cutoff(string *cutoff); void cutoff(string *cutoff);
//! Adds a weight of the "markowitz" criteria statement
void markowitz(string *markowitz);
//! Adds a "dsample" statement //! Adds a "dsample" statement
void dsample(string *arg1); void dsample(string *arg1);
//! Adds a "dsample" statement //! Adds a "dsample" statement
@ -211,9 +213,9 @@ public:
void add_covar_shock(string *var1, string *var2, NodeID value); void add_covar_shock(string *var1, string *var2, NodeID value);
//! Adds a correlated chock //! Adds a correlated chock
void add_correl_shock(string *var1, string *var2, NodeID value); void add_correl_shock(string *var1, string *var2, NodeID value);
//! Adds a shock period range //! Adds a shock period range
void add_period(string *p1, string *p2); void add_period(string *p1, string *p2);
//! Adds a shock period //! Adds a shock period
void add_period(string *p1); void add_period(string *p1);
//! Adds a shock value (when only a numerical constant) //! Adds a shock value (when only a numerical constant)
void add_value_const(string *value); void add_value_const(string *value);

225
tests/ferhat/MARK3_endo.dat
View File

File diff suppressed because one or more lines are too long

35
tests/ferhat/MARK3_exo.dat
View File

File diff suppressed because one or more lines are too long

2
tests/ferhat/ireland.mod
View File

@ -82,7 +82,7 @@ options_.dynatol=4e-8;
shocks; shocks;
var e; var e;
periods 1; periods 1;
values -0.0002; values 0.0002;
end; end;
simul(periods=60); simul(periods=60);
rplot y; rplot y;

11
tests/ferhat/multimod.mod
View File

@ -2655,7 +2655,7 @@ W0906=0.0800069594276;
W0907=0.147854375051; W0907=0.147854375051;
W0908=0.206834342322; W0908=0.206834342322;
W0909=-1; W0909=-1;
model(SPARSE_DLL,gcc_compiler); model(SPARSE_DLL,gcc_compiler,markowitz=2.0);
( log(US_CPI)-(log(US_CPI(-1)))) = US_CPI1*( log(US_PIM)-(log(US_PIM(-1))))+US_CPI2*( log(US_PGNP)-(log(US_PGNP(-1))))+(1-US_CPI1-US_CPI2)*log(US_CPI(-1)/US_CPI(-2))+RES_US_CPI ; ( log(US_CPI)-(log(US_CPI(-1)))) = US_CPI1*( log(US_PIM)-(log(US_PIM(-1))))+US_CPI2*( log(US_PGNP)-(log(US_PGNP(-1))))+(1-US_CPI1-US_CPI2)*log(US_CPI(-1)/US_CPI(-2))+RES_US_CPI ;
US_UNR_A = US_UNR_FE+US_UNR_1*100*log(US_GDP/US_GDP_FE)+US_UNR_2*(US_UNR(-1)-US_UNR_FE(-1))+RES_US_UNR_A ; US_UNR_A = US_UNR_FE+US_UNR_1*100*log(US_GDP/US_GDP_FE)+US_UNR_2*(US_UNR(-1)-US_UNR_FE(-1))+RES_US_UNR_A ;
US_UNR = /*MAX(US_UNR_A;0.1)*/US_UNR_A ; US_UNR = /*MAX(US_UNR_A;0.1)*/US_UNR_A ;
@ -4551,8 +4551,15 @@ options_.slowc = 1.0;
options_.dynatol = 1e-4; options_.dynatol = 1e-4;
simul(periods=117,datafile=mark3); simul(periods=50,datafile=mark3);
shocks;
var US_G;
periods 1;
values 4330.714737;
end;
simul(periods=50);
rplot WTRADER; rplot WTRADER;
rplot US_GDP; rplot US_GDP;

619
tests/ferhat/simulate.cc
View File

@ -1,4 +1,4 @@
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// simulate.c // // simulate.c //
// simulate file designed for GNU GCC C++ compiler // // simulate file designed for GNU GCC C++ compiler //
// use GCC_COMPILER option in MODEL command // // use GCC_COMPILER option in MODEL command //
@ -7,6 +7,7 @@
//#define PRINT_OUT //#define PRINT_OUT
#define INDIRECT_SIMULATE #define INDIRECT_SIMULATE
#define PRINT_OUT_p #define PRINT_OUT_p
#define MARKOVITZ
//#define PRINT_OUT_y1 //#define PRINT_OUT_y1
//#define PRINT_u //#define PRINT_u
//#define PRINT_OUT_b //#define PRINT_OUT_b
@ -14,6 +15,7 @@
//#define EXTENDED //#define EXTENDED
//#define FLOAT //#define FLOAT
//#define WRITE_u //#define WRITE_u
//#define MEMORY_LEAKS
//#define N_MX_ALLOC //#define N_MX_ALLOC
//#define MEM_ALLOC_CHK //#define MEM_ALLOC_CHK
#define NEW_ALLOC #define NEW_ALLOC
@ -56,6 +58,13 @@ typedef struct NonZeroElem
NonZeroElem *NZE_R_N, *NZE_C_N; NonZeroElem *NZE_R_N, *NZE_C_N;
}; };
typedef struct t_save_op_s
{
short int lag, operat;
int first, second;
};
/*#include "simulate.hh"*/ /*#include "simulate.hh"*/
@ -70,11 +79,11 @@ int *pivot, *pivotk;
longd *pivotv, *pivotva=NULL; longd *pivotv, *pivotva=NULL;
bool *line_done; bool *line_done;
bool symbolic=true, alt_symbolic=false, record_all=false; bool symbolic=true, alt_symbolic=false, record_all=false;
long int nop_all, nopa_all; long int nop_all, nopa_all, nop1, nop2;
const longd very_big=1e24; const longd very_big=1e24;
int Per_y_, Per_u_, it_, nb_row_x, u_size, y_size, x_size, y_kmin, y_kmax, y_decal; int Per_y_, Per_u_, it_, nb_row_x, u_size, y_size, x_size, y_kmin, y_kmax, y_decal;
int periods, maxit_, max_u=0, min_u=0x7FFFFFFF, g_nop_all=0; int periods, maxit_, max_u=0, min_u=0x7FFFFFFF, g_nop_all=0;
double *params; longd *params, markowitz_c;
longd *u, *y, *x, *r, *g1, *g2, *ya, *direction; longd *u, *y, *x, *r, *g1, *g2, *ya, *direction;
longd slowc, slowc_save, solve_tolf, max_res, res1, res2, res1a=9.0e60; longd slowc, slowc_save, solve_tolf, max_res, res1, res2, res1a=9.0e60;
bool cvg, print_err, swp_f; bool cvg, print_err, swp_f;
@ -82,6 +91,8 @@ int swp_f_b=0;
pctimer_t t0, t1; pctimer_t t0, t1;
int u_count_alloc, u_count_alloc_save, size_of_direction; int u_count_alloc, u_count_alloc_save, size_of_direction;
int i, j, k, nb_endo, u_count, u_count_init, iter; int i, j, k, nb_endo, u_count, u_count_init, iter;
const int alt_symbolic_count_max=1;
int alt_symbolic_count=0;
longd err; longd err;
std::string filename; std::string filename;
NonZeroElem **FNZE_R, **FNZE_C; NonZeroElem **FNZE_R, **FNZE_C;
@ -105,17 +116,6 @@ typedef struct t_table_u
std::fstream SaveCode, SaveCode_swp; std::fstream SaveCode, SaveCode_swp;
/*longd pow1(longd a, longd b)
{
if(a<0)
{
mexPrintf("Attempt to compute (-X)^Y at time %d\n",it_);
return(-pow(-a,b));
}
else
return(pow(a,b));
}
*/
longd pow1(longd a, longd b) longd pow1(longd a, longd b)
{ {
@ -1085,6 +1085,7 @@ class SparseMatrix
int At_Col(int c, int lag, NonZeroElem **first); int At_Col(int c, int lag, NonZeroElem **first);
int NRow(int r); int NRow(int r);
int NCol(int c); int NCol(int c);
int Union_Row(int row1, int row2);
void Print(int Size,int *b); void Print(int Size,int *b);
int Get_u(); int Get_u();
void Delete_u(int pos); void Delete_u(int pos);
@ -1212,6 +1213,39 @@ int SparseMatrix::At_Row(int r, NonZeroElem **first)
return NbNZRow[r]; return NbNZRow[r];
} }
int
SparseMatrix::Union_Row(int row1, int row2)
{
NonZeroElem *first1, *first2;
int n1=At_Row(row1, &first1);
int n2=At_Row(row2, &first2);
int i1=0, i2=0, nb_elem=0;
while(i1<n1 && i2<n2)
{
if(first1->c_index==first2->c_index)
{
nb_elem++;
i1++;
i2++;
first1=first1->NZE_R_N;
first2=first2->NZE_R_N;
}
else if(first1->c_index<first2->c_index)
{
nb_elem++;
i1++;
first1=first1->NZE_R_N;
}
else
{
nb_elem++;
i2++;
first2=first2->NZE_R_N;
}
}
return nb_elem;
}
int int
SparseMatrix::At_Pos(int r, int c, NonZeroElem **first) SparseMatrix::At_Pos(int r, int c, NonZeroElem **first)
{ {
@ -1940,6 +1974,7 @@ void SparseMatrix::End(int Size)
{ {
mxFree(NZE_Mem_add[i*CHUNK_BLCK_SIZE]); mxFree(NZE_Mem_add[i*CHUNK_BLCK_SIZE]);
} }
mxFree(NZE_Mem_add);
init_Mem(); init_Mem();
#else #else
for (int i=0;i<Size*periods;i++) for (int i=0;i<Size*periods;i++)
@ -1959,16 +1994,20 @@ void SparseMatrix::End(int Size)
mxFree(NbNZCol); mxFree(NbNZCol);
mxFree(b); mxFree(b);
mxFree(line_done); mxFree(line_done);
mxFree(pivot);
mxFree(pivotk);
mxFree(pivotv);
mxFree(pivotva);
} }
bool bool
compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size, long int *ndiv, long int *nsub) compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size, long int *ndiv, long int *nsub)
{ {
//mexPrintf("=>in compare beg_t=%d\n",beg_t); //mexPrintf("=>in compare beg_t=%d\n",beg_t);
long int i,j,nop=nop4/4, t, index_d, k; long int i,j,/*nop=nop4/4*/nop=nop4/2, t, index_d, k;
longd r=0.0; longd r=0.0;
bool OK=true; bool OK=true;
t_save_op_s *save_op_s, *save_opa_s, *save_opaa_s;
int *diff1, *diff2; int *diff1, *diff2;
//mexPrintf("nop=%d\n",nop); //mexPrintf("nop=%d\n",nop);
//g_save_op=(int*)mxMalloc(nop*5*sizeof(int)); //g_save_op=(int*)mxMalloc(nop*5*sizeof(int));
@ -1983,70 +2022,99 @@ compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, lo
#ifdef MEM_ALLOC_CHK #ifdef MEM_ALLOC_CHK
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
j=k=0; j=k=i=0;
for (i=0;i<nop4 && OK;i+=4) while(i<nop4 && OK)
{
save_op_s=(t_save_op_s*)&(save_op[i]);
save_opa_s=(t_save_op_s*)&(save_opa[i]);
save_opaa_s=(t_save_op_s*)&(save_opaa[i]);
//mexPrintf("i=%d nop4=%d save_op_s->operat=%d\n",i,nop4,save_op_s->operat);
diff1[j]=save_op_s->first-save_opa_s->first;
switch(save_op_s->operat)
{
case FLD:
case FDIV:
OK=(save_op_s->operat==save_opa_s->operat && save_opa_s->operat==save_opaa_s->operat
&& diff1[j]==(save_opa_s->first-save_opaa_s->first));
i+=2;
break;
case FLESS:
case FSUB:
diff2[j]=save_op_s->second-save_opa_s->second;
OK=(save_op_s->operat==save_opa_s->operat && save_opa_s->operat==save_opaa_s->operat
&& diff1[j]==(save_opa_s->first-save_opaa_s->first)
&& diff2[j]==(save_opa_s->second-save_opaa_s->second));
i+=3;
break;
default:
mexPrintf("unknown operator = %d ",save_op_s->operat);
filename+=" stopped";
mexErrMsgTxt(filename.c_str());
break;
}
j++;
}
/*j=k=i=0;
save_op_s=(t_save_op_s*)&(*save_op);
save_opa_s=(t_save_op_s*)&(*save_opa);
save_opaa_s=(t_save_op_s*)&(*save_opaa);
while(i<nop4 && OK)
{
mexPrintf("i=%d nop4=%d save_op_s->operat=%d\n",i,nop4,save_op_s->operat);
diff1[j]=save_op_s->first-save_opa_s->first;
switch(save_op_s->operat)
{
case FLD:
case FDIV:
OK=(save_op_s->operat==save_opa_s->operat && save_opa_s->operat==save_opaa_s->operat
&& diff1[j]==(save_opa_s->first-save_opaa_s->first));
i+=2;
if(i<nop4)
{
save_op_s+=2;
save_opa_s+=2;
save_opaa_s+=2;
}
break;
case FLESS:
case FSUB:
diff2[j]=save_op_s->second-save_opa_s->second;
OK=(save_op_s->operat==save_opa_s->operat && save_opa_s->operat==save_opaa_s->operat
&& diff1[j]==(save_opa_s->first-save_opaa_s->first)
&& diff2[j]==(save_opa_s->second-save_opaa_s->second));
i+=3;
if(i<nop4)
{
save_op_s+=3;
save_opa_s+=3;
save_opaa_s+=3;
}
break;
default:
mexPrintf("unknown operator = %d ",save_op_s->operat);
filename+=" stopped";
mexErrMsgTxt(filename.c_str());
break;
}
j++;
}
*/
/*for (i=0;i<nop4 && OK;i+=4)
{ {
diff1[j]=save_op[i+1]-save_opa[i+1]; diff1[j]=save_op[i+1]-save_opa[i+1];
diff2[j]=save_op[i+2]-save_opa[i+2]; diff2[j]=save_op[i+2]-save_opa[i+2];
/*g_save_op[k]=save_op[i];
g_save_op[k+1]=save_op[i+1]-beg_t*diff1[j];
g_save_op[k+2]=diff1[j];
if(save_op[i]==FSUB || save_op[i]==FLESS)
{
g_save_op[k+3]=save_op[i+2]-beg_t*diff2[j];
g_save_op[k+4]=diff2[j];
}
k+=5;*/
/*if(save_op[i+2]==48616 )
mexPrintf("save_op[%d]=%d save_opa[%d]=%d save_opa[%d]=%d diff2[%d]=%d\n",i+2,save_op[i+2],i+2,save_opa[i+2],i+2,save_opaa[i+2],j,diff2[j]);*/
OK=(save_op[i]==save_opa[i] && save_opa[i]==save_opaa[i] OK=(save_op[i]==save_opa[i] && save_opa[i]==save_opaa[i]
&& diff1[j]==(save_opa[i+1]-save_opaa[i+1]) && diff1[j]==(save_opa[i+1]-save_opaa[i+1])
&& diff2[j]==(save_opa[i+2]-save_opaa[i+2])); && diff2[j]==(save_opa[i+2]-save_opaa[i+2]));
//mexPrintf("diff1[%d]=%d, diff2[%d]=%d\n",j,diff1[j],j,diff2[j]);
/*mexPrintf("save_op[%d]=%d, save_opa[%d]=%d, save_opaa[%d]=%d\n",i,save_op[i],i,save_opa[i],i,save_opaa[i]);
mexPrintf("save_op[%d+1]=%d, save_opa[%d+1]=%d, save_opaa[%d+1]=%d\n",i,save_op[i+1],i,save_opa[i+1],i,save_opaa[i+1]);
mexPrintf("save_op[%d+2]=%d, save_opa[%d+2]=%d, save_opaa[%d+2]=%d\n",i,save_op[i+2],i,save_opa[i+2],i,save_opaa[i+2]);*/
j++; j++;
}
//mexPrintf("j=%d nop=%d\n",j,nop);
//OK=false;
/*if(!OK)
{
mexPrintf("save_op[%d](%d)==save_opa[%d](%d)==save_opaa[%d](%d)\n diff1[%d](%d)==(save_opa[%d]-save_opaa[%d])(%d)\n diff2[%d](%d)==(save_opa[%d]-save_opaa[%d])(%d)\n",
i,save_op[i],i,save_opa[i],i,save_opaa[i],j,diff1[j],i+1,i+1,(save_opa[i+1]-save_opaa[i+1]),j,diff2[j],i+2,i+2,(save_opa[i+2]-save_opaa[i+2]));
}*/ }*/
//g_nop_all=k;
//mexPrintf("1 - OK=%d t=%d\n",OK,beg_t);
if (OK) if (OK)
{ for (i=beg_t;i<periods;i++)
/*for(i=beg_t-2;i<beg_t;i++) for (j=0;j<Size;j++)
{ pivot[i*Size+j]=pivot[(i-1)*Size+j]+Size;
for(j=0;j<Size;j++) /*mexPrintf("2 - OK=%d t=%d\n",OK,beg_t);
{ mexPrintf("from beg_t=%d to periods=%d\n",beg_t, periods);*/
mexPrintf("%d==%d ",pivot[i*Size+j]+Size,pivot[(i+1)*Size+j]);
OK = OK && ((pivot[i*Size+j]+Size)==pivot[(i+1)*Size+j]);
}
mexPrintf("\n");
}
if(OK)
{*/
for (i=beg_t;i<periods;i++)
{
/*if(i+1==periods)
mexPrintf("t=%d ",i);*/
for (j=0;j<Size;j++)
{
pivot[i*Size+j]=pivot[(i-1)*Size+j]+Size;
/*if(i+1==periods)
mexPrintf("%d==%d ",pivot[(i-1)*Size+j]+Size,pivot[i*Size+j]);*/
}
/*if(i+1==periods)
mexPrintf("\n");*/
}
/*}*/
}
//mexPrintf("2 - OK=%d t=%d\n",OK,beg_t);
//mexPrintf("from beg_t=%d to periods=%d\n",beg_t, periods);
if (OK) if (OK)
{ {
@ -2057,10 +2125,140 @@ compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, lo
#endif #endif
t=1; t=1;
//mexPrintf("nop4=%d\n",nop4);
for (;t<periods-beg_t-max(y_kmax,y_kmin);t++) for (;t<periods-beg_t-max(y_kmax,y_kmin);t++)
{ {
//mexPrintf("t=%d\n",t); i=j=0;
while(i<nop4)
{
save_op_s=(t_save_op_s*)(&(save_op[i]));
index_d=save_op_s->first+t*diff1[j];
if (index_d>=u_count_alloc)
{
u_count_alloc+=5*u_count_alloc_save;
#ifdef MEM_ALLOC_CHK
mexPrintf("u=(longd*)mxRealloc(u,u_count_alloc*sizeof(longd))\n",u_count_alloc);
#endif
u=(longd*)mxRealloc(u,u_count_alloc*sizeof(longd));
#ifdef MEM_ALLOC_CHK
mexPrintf("ok\n");
#endif
if (!u)
{
mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n",u_count_alloc*sizeof(longd));
mexErrMsgTxt("Exit from Dynare");
}
}
switch (save_op_s->operat)
{
case FLD :
r=u[index_d];
#ifdef PRINT_u
mexPrintf("FLD u[%d] (%f)\n",index_d,u[index_d]);
#endif
i+=2;
break;
case FDIV :
u[index_d]/=r;
#ifdef PRINT_u
mexPrintf("FDIV u[%d](%f)/=r(%f)=(%f)\n",index_d,u[index_d],r,u[index_d]);
#endif
i+=2;
break;
case FSUB :
u[index_d]-=u[save_op_s->second+t*diff2[j]]*r;
#ifdef PRINT_u
mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,save_op_s->second+t*diff2[j],u[save_op_s->second+t*diff2[j]],r,u[index_d] );
#endif
i+=3;
break;
case FLESS:
u[index_d]=-u[save_op_s->second+t*diff2[j]]*r;
#ifdef PRINT_u
mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,save_op_s->second+t*diff2[j],u[save_op_s->second+t*diff2[j]],r,u[index_d] );
#endif
i+=3;
break;
}
j++;
}
}
int t1=t;
for (t=t1;t<periods-beg_t;t++)
{
i=j=0;
while(i<nop4)
{
save_op_s=(t_save_op_s*)(&(save_op[i]));
if(save_op_s->lag<((periods-beg_t)-t))
{
index_d=save_op_s->first+t*diff1[j];
if (index_d>=u_count_alloc)
{
u_count_alloc+=5*u_count_alloc_save;
#ifdef MEM_ALLOC_CHK
mexPrintf("u=(longd*)mxRealloc(u,u_count_alloc*sizeof(longd))\n",u_count_alloc);
#endif
u=(longd*)mxRealloc(u,u_count_alloc*sizeof(longd));
#ifdef MEM_ALLOC_CHK
mexPrintf("ok\n");
#endif
if (!u)
{
mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n",u_count_alloc*sizeof(longd));
mexErrMsgTxt("Exit from Dynare");
}
}
switch (save_op_s->operat)
{
case FLD :
r=u[index_d];
#ifdef PRINT_u
mexPrintf("FLD u[%d] (%f)\n",index_d,u[index_d]);
#endif
i+=2;
break;
case FDIV :
u[index_d]/=r;
#ifdef PRINT_u
mexPrintf("FDIV u[%d](%f)/=r(%f)=(%f)\n",index_d,u[index_d],r,u[index_d]);
#endif
i+=2;
break;
case FSUB :
u[index_d]-=u[save_op_s->second+t*diff2[j]]*r;
#ifdef PRINT_u
mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,save_op_s->second+t*diff2[j],u[save_op_s->second+t*diff2[j]],r,u[index_d] );
#endif
i+=3;
break;
case FLESS:
u[index_d]=-u[save_op_s->second+t*diff2[j]]*r;
#ifdef PRINT_u
mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,save_op_s->second+t*diff2[j],u[save_op_s->second+t*diff2[j]],r,u[index_d] );
#endif
i+=3;
break;
}
}
else
{
switch (save_op_s->operat)
{
case FLD :
case FDIV :
i+=2;
break;
case FSUB :
case FLESS :
i+=3;
break;
}
}
j++;
}
}
/* for (;t<periods-beg_t-max(y_kmax,y_kmin);t++)
{
#ifdef WRITE_u #ifdef WRITE_u
toto << "t=" << t << endl; toto << "t=" << t << endl;
#endif #endif
@ -2098,36 +2296,18 @@ compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, lo
#ifdef PRINT_u #ifdef PRINT_u
mexPrintf("FDIV u[%d](%f)/=r(%f)=(%f)\n",index_d,u[index_d],r,u[index_d]); mexPrintf("FDIV u[%d](%f)/=r(%f)=(%f)\n",index_d,u[index_d],r,u[index_d]);
#endif #endif
/*if((t>=periods-beg_t-y_kmax))
mexPrintf("u[%d]/=%f=%f\n",index_d,double(r),double(u[index_d]));*/
/*toto << i_toto << " u[" << index_d << "]/=" << r << "=" << u[index_d] << endl;
i_toto++;*/
//(*ndiv)++;
break; break;
case FSUB : case FSUB :
u[index_d]-=u[save_op[i+2]+t*diff2[j]]*r; u[index_d]-=u[save_op[i+2]+t*diff2[j]]*r;
#ifdef PRINT_u #ifdef PRINT_u
mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],r,u[index_d] ); mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],r,u[index_d] );
#endif #endif
/*if((t==periods-beg_t-y_kmax))
mexPrintf("u[%d]-=u[%d]*%f=%f\n",index_d,save_op[i+2]+t*diff2[j],double(r),double(u[index_d]));*/
/*toto << i_toto << " u[" << index_d << "]-=u[" << save_op[i+2]+t*diff2[j] << "]*" << r << "=" << u[index_d] << endl;
i_toto++;*/
//(*nsub)++;
break; break;
case FLESS: case FLESS:
u[index_d]=-u[save_op[i+2]+t*diff2[j]]*r; u[index_d]=-u[save_op[i+2]+t*diff2[j]]*r;
#ifdef PRINT_u #ifdef PRINT_u
mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],r,u[index_d] ); mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],r,u[index_d] );
#endif #endif
/*if((t==periods-beg_t-y_kmax))
mexPrintf("u[%d]=-u[%d]*%f=%f\n",index_d,save_op[i+2]+t*diff2[j],double(r),double(u[index_d]));*/
/*toto << i_toto << "u[" << index_d << "]=-u[" << save_op[i+2]+t*diff2[j] << "]*" << r << "=" << u[index_d] << endl;
i_toto++;*/
//(*nsub)++;
break; break;
} }
j++; j++;
@ -2138,7 +2318,7 @@ compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, lo
{ {
j=0; j=0;
#ifdef WRITE_u #ifdef WRITE_u
toto << "t=" << t << " t+beg_t=" << t+beg_t/*<< "----------------------------------------------------------------------------------------------------------------"*/ <<endl; toto << "t=" << t << " t+beg_t=" << t+beg_t <<endl;
#endif #endif
//mexPrintf("t=%d----------------------------------------------------------------------------------------------------------------\n",t); //mexPrintf("t=%d----------------------------------------------------------------------------------------------------------------\n",t);
for (i=0;i<nop4;i+=4) for (i=0;i<nop4;i+=4)
@ -2177,10 +2357,9 @@ compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, lo
mexPrintf("FDIV u[%d](%f)/=r(%f)=(%f)\n",index_d,u[index_d],r,u[index_d]); mexPrintf("FDIV u[%d](%f)/=r(%f)=(%f)\n",index_d,u[index_d],r,u[index_d]);
#endif #endif
/*if(index_d==81726)
mexPrintf("(0) u[%d]/=%f=%f save_op[%d]=%d diff1[%d]=%d\n",index_d,double(r),double(u[index_d]),i+1,save_op[i+1],j,diff1[j]);*/
#ifdef WRITE_u #ifdef WRITE_u
toto << i_toto << " u[" /*<< index_d*/ << "]/=" << r << "=" << u[index_d] << /*" lag_index=" << save_op[i+3] << " periods-beg_t-t=" << periods-beg_t-t <<*/ endl; toto << i_toto << " u[" << "]/=" << r << "=" << u[index_d] << endl;
i_toto++; i_toto++;
#endif #endif
//(*ndiv)++; //(*ndiv)++;
@ -2191,10 +2370,8 @@ compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, lo
#ifdef PRINT_u #ifdef PRINT_u
mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],r,u[index_d] ); mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],r,u[index_d] );
#endif #endif
/*if(index_d==81726)
mexPrintf("(1) u[%d]-=u[%d](%f)*%f=%f save_op[%d]=%d diff1[%d]=%d\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],double(r),double(u[index_d]),i+1,save_op[i+1],j,diff1[j]);*/
#ifdef WRITE_u #ifdef WRITE_u
toto << i_toto << " u[" /*<< index_d*/ << "]-=u[" /*<< save_op[i+2]+t*diff2[j]*/ << "]*" << r << "=" << u[index_d] << /*" lag_index=" << save_op[i+3] << " periods-beg_t-t=" << periods-beg_t-t <<*/ endl; toto << i_toto << " u[" << "]-=u[" << "]*" << r << "=" << u[index_d] << endl;
i_toto++; i_toto++;
#endif #endif
//(*nsub)++; //(*nsub)++;
@ -2205,10 +2382,8 @@ compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, lo
#ifdef PRINT_u #ifdef PRINT_u
mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],r,u[index_d] ); mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],r,u[index_d] );
#endif #endif
/*if(index_d==81726)
mexPrintf("(3) u[%d]=-u[%d](%f)*r(%f)=%f save_op[%d]=%d diff1[%d]=%d\n",index_d,save_op[i+2]+t*diff2[j],u[save_op[i+2]+t*diff2[j]],double(r),double(u[index_d]),i+1,save_op[i+1],j,diff1[j]);*/
#ifdef WRITE_u #ifdef WRITE_u
toto << i_toto << " u[" /*<< index_d*/ << "]=-u[" /*<< save_op[i+2]+t*diff2[j]*/ << "]*" << r << "=" << u[index_d] << /*" lag_index=" << save_op[i+3] << " periods-beg_t-t=" << periods-beg_t-t <<*/ endl; toto << i_toto << " u[" << "]=-u[" << "]*" << r << "=" << u[index_d] << endl;
i_toto++; i_toto++;
#endif #endif
//(*nsub)++; //(*nsub)++;
@ -2217,7 +2392,7 @@ compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, lo
} }
j++; j++;
} }
} }*/
#ifdef WRITE_u #ifdef WRITE_u
toto.close(); toto.close();
filename+=" stopped"; filename+=" stopped";
@ -2459,9 +2634,6 @@ SparseMatrix::complete(int beg_t, int Size, int periods, int *b)
int *diff; int *diff;
longd yy=0.0, err, ferr; longd yy=0.0, err, ferr;
/*mexPrintf("in complete t=%d\n",beg_t);
mexPrintf("u_count=%d\n",u_count);*/
int size_of_save_code=(1+y_kmax)*Size*(Size+1+4)/2*4; int size_of_save_code=(1+y_kmax)*Size*(Size+1+4)/2*4;
#ifdef MEM_ALLOC_CHK #ifdef MEM_ALLOC_CHK
mexPrintf("save_code=(int*)mxMalloc(%d*sizeof(int))\n",size_of_save_code); mexPrintf("save_code=(int*)mxMalloc(%d*sizeof(int))\n",size_of_save_code);
@ -2475,7 +2647,6 @@ SparseMatrix::complete(int beg_t, int Size, int periods, int *b)
#ifdef MEM_ALLOC_CHK #ifdef MEM_ALLOC_CHK
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
//mexPrintf("Size=%d\n",Size);
cal_y=y_size*y_kmin; cal_y=y_size*y_kmin;
i=(beg_t+1)*Size-1; i=(beg_t+1)*Size-1;
@ -2818,9 +2989,12 @@ void chk_avail_mem(int **save_op_all,long int *nop_all,long int *nopa_all,int ad
} }
int int
simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it) simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg)
{ {
/*Triangularisation at each period of a block using a simple gaussian Elimination*/ /*Triangularisation at each period of a block using a simple gaussian Elimination*/
const int maskhi=0xFFFF0000;
const int masklo=0x0000FFFF;
t_save_op_s *save_op_s;
bool record=false; bool record=false;
int *save_op=NULL, *save_opa=NULL, *save_opaa=NULL; int *save_op=NULL, *save_opa=NULL, *save_opaa=NULL;
long int nop=0, nopa=0, nopaa=0; long int nop=0, nopa=0, nopaa=0;
@ -2840,6 +3014,39 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
pctimer_t t01; pctimer_t t01;
pctimer_t t1 = pctimer(); pctimer_t t1 = pctimer();
tdelete1=0; tdelete2=0; tdelete21=0; tdelete22=0; tdelete221=0; tdelete222=0; tdelete1=0; tdelete2=0; tdelete21=0; tdelete22=0; tdelete221=0; tdelete222=0;
/*int Z[3];
int *p;
p=Z;
mexPrintf("p=%x\n",p);
p=p+1;
mexPrintf("p=%x\n",p);*/
/*
int Z[3];
Z[0]=0x56CE000A;
Z[1]=2;
Z[2]=3;
save_op_s=(t_save_op_s*)Z;
*/
/*for(i=0;i<1;i++)
mexPrintf("Z[i]=%d lo(Z)=%d hi(Z)=%d\n",Z[i],Z[i] & masklo, Z[i] & maskhi);*/
/* mexPrintf("save_op_s->lag=%d sizeof(save_op_s->lag)=%d\n",int(save_op_s->lag),sizeof(save_op_s->lag));
mexPrintf("save_op_s->operat=%d sizeof(save_op_s->operat)=%d\n",int(save_op_s->operat),sizeof(save_op_s->operat));
mexPrintf("save_op_s->first=%d sizeof(save_op_s->first)=%d\n",save_op_s->first,sizeof(save_op_s->first));
mexPrintf("save_op_s->second=%d sizeof(save_op_s->second)=%d\n",save_op_s->second,sizeof(save_op_s->second));
mexPrintf("save_op_s->lag=%d sizeof(save_op_s->lag)=%d\n",int(save_op_s->lag),sizeof(save_op_s->lag));
*/
/*int z=0x56CECDEF;*/
/*mexPrintf("lo(z)=%d hi(z)=%d\n",z & masklo, z & maskhi);*/
/*filename+=" stopped";
mexErrMsgTxt(filename.c_str());*/
#ifdef MEMORY_LEAKS
mexEvalString("feature('memstats');");
#endif
if (isnan(res1) || isinf(res1)) if (isnan(res1) || isinf(res1))
{ {
/*if(record_all and nop_all) /*if(record_all and nop_all)
@ -2870,11 +3077,12 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
/*}*/ /*}*/
} }
u_count+=u_count_init; u_count+=u_count_init;
if (alt_symbolic) if (alt_symbolic && alt_symbolic_count<alt_symbolic_count_max)
{ {
mexPrintf("Pivoting method will be only applied to the first periods.\n"); mexPrintf("Pivoting method will be applied only to the first periods.\n");
alt_symbolic=false; alt_symbolic=false;
symbolic=true; symbolic=true;
alt_symbolic_count++;
} }
if (((res1/res1a-1)>-0.1) && symbolic) if (((res1/res1a-1)>-0.1) && symbolic)
{ {
@ -2891,6 +3099,8 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
mexPrintf("-----------------------------------\n"); mexPrintf("-----------------------------------\n");
//mexPrintf("record_all=%d save_op_all=%x\n",record_all,save_op_all); //mexPrintf("record_all=%d save_op_all=%x\n",record_all,save_op_all);
mexEvalString("drawnow;"); mexEvalString("drawnow;");
if(cvg)
return(0);
#ifdef PRINT_OUT #ifdef PRINT_OUT
mexPrintf("Size=%d y_size=%d y_kmin=%d y_kmax=%d u_count=%d u_count_alloc=%d periods=%d\n",Size,y_size,y_kmin,y_kmax,u_count,u_count_alloc,periods); mexPrintf("Size=%d y_size=%d y_kmin=%d y_kmax=%d u_count=%d u_count_alloc=%d periods=%d\n",Size,y_size,y_kmin,y_kmax,u_count,u_count_alloc,periods);
#endif #endif
@ -2979,6 +3189,7 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
#else #else
save_op=(int*)mxMalloc(nop*sizeof(int)); save_op=(int*)mxMalloc(nop*sizeof(int));
#endif #endif
nopa=nop;
#ifdef MEM_ALLOC_CHK #ifdef MEM_ALLOC_CHK
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
@ -3010,8 +3221,14 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
#ifdef PRINT_OUT #ifdef PRINT_OUT
mexPrintf("nb_eq=%d\n",nb_eq); mexPrintf("nb_eq=%d\n",nb_eq);
#endif #endif
if /*(t<=y_kmin)*/((symbolic && t<=y_kmin) || !symbolic) if /*(t<=y_kmin) */((symbolic && t<=y_kmin) || !symbolic)
{ {
#ifdef MARKOVITZ
double piv_v[Size];
int pivj_v[Size], pivk_v[Size], NR[Size], l=0, N_max=0;
bool one=false;
piv_abs=0;
#endif
for (j=0;j<nb_eq;j++) for (j=0;j<nb_eq;j++)
{ {
#ifdef PRINT_OUT #ifdef PRINT_OUT
@ -3026,6 +3243,39 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
#endif #endif
int jj=first->r_index; int jj=first->r_index;
int NRow_jj=sparse_matrix.NRow(jj); int NRow_jj=sparse_matrix.NRow(jj);
//Union_Row(jj, int row2)
#ifdef MARKOVITZ
piv_v[l]=u[k];
longd piv_fabs=fabs(u[k]);
//mexPrintf("piv_v[%d]=%f\n",l,piv_v[l]);
pivj_v[l]=jj;
pivk_v[l]=k;
NR[l]=NRow_jj;
if(NRow_jj==1 && !one)
{
one=true;
piv_abs=piv_fabs;
N_max=NRow_jj;
}
if(!one)
{
if(piv_fabs>piv_abs)
piv_abs=piv_fabs;
if(NRow_jj>N_max)
N_max=NRow_jj;
}
else
{
if(NRow_jj==1)
{
if(piv_fabs>piv_abs)
piv_abs=piv_fabs;
if(NRow_jj>N_max)
N_max=NRow_jj;
}
}
l++;
#else
if (piv_abs<fabs(u[k])||NRow_jj==1) if (piv_abs<fabs(u[k])||NRow_jj==1)
{ {
piv=u[k]; piv=u[k];
@ -3035,10 +3285,48 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
if (NRow_jj==1) if (NRow_jj==1)
break; break;
} }
#endif
//ncomp++; //ncomp++;
} }
first=first->NZE_C_N; first=first->NZE_C_N;
} }
#ifdef MARKOVITZ
double markovitz=0, markovitz_max=-9e70;
//mexPrintf("N_max=%d piv_abs=%f l=%d nb_eq=%d\n",N_max,piv_abs,l,nb_eq);
if(!one)
{
for(j=0;j<l;j++)
{
markovitz=exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(double(NR[j])/double(N_max)));
//mexPrintf("none j=%d markovitz=%f piv_v=%f NR=%d term1=%f term2=%f\n",j,markovitz,piv_v[j],NR[j],log(fabs(piv_v[j])/piv_abs),log(double(NR[j])/double(N_max)));
if(markovitz>markovitz_max)
{
piv=piv_v[j];
//piv_abs=fabs(piv);
pivj=pivj_v[j]; //Line number
pivk=pivk_v[j]; //positi
markovitz_max=markovitz;
}
}
}
else
{
for(j=0;j<l;j++)
{
markovitz=exp(log(fabs(piv_v[j])/piv_abs)-markowitz_c*log(NR[j]/N_max));
//mexPrintf("one j=%d markovitz=%f piv_v=%f NR=%d\n",j,markovitz,piv_v[j],NR[j]);
if(markovitz>markovitz_max && NR[j]==1)
{
piv=piv_v[j];
//piv_abs=fabs(piv);
pivj=pivj_v[j]; //Line number
pivk=pivk_v[j]; //positi
markovitz_max=markovitz;
}
}
}
//smexPrintf("selected pivj=%d piv=%f piv_abs=%f markovitz_max=%f one=%d \n",pivj,piv,piv_abs,markovitz_max,one);
#endif
#ifdef PROFILER #ifdef PROFILER
tpivot+=pctimer()-td0; tpivot+=pctimer()-td0;
#endif #endif
@ -3071,7 +3359,7 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
{ {
if (record) if (record)
{ {
if (nop+2>=nopa) if (nop+1>=nopa)
{ {
//mexPrintf("Error: out of save_op[%d] nopa=%d\n",nop+2,nopa); //mexPrintf("Error: out of save_op[%d] nopa=%d\n",nop+2,nopa);
nopa=int(1.5*nopa); nopa=int(1.5*nopa);
@ -3087,12 +3375,17 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
} }
save_op[nop]=FLD; save_op_s=(t_save_op_s*)(&(save_op[nop]));
save_op_s->operat=FLD;
save_op_s->first=pivk;
save_op_s->lag=0;
/*save_op[nop]=FLD;
save_op[nop+1]=pivk; save_op[nop+1]=pivk;
save_op[nop+2]=0; save_op[nop+2]=0;
save_op[nop+3]=0; save_op[nop+3]=0;*/
} }
nop+=4; /*nop+=4;*/
nop+=2;
} }
else if (record_all) else if (record_all)
{ {
@ -3155,7 +3448,7 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
{ {
if (record) if (record)
{ {
if (nop+2>=nopa) if (nop+1>=nopa)
{ {
//mexPrintf("Error: out of save_op[%d] nopa=%d\n",nop+2,nopa); //mexPrintf("Error: out of save_op[%d] nopa=%d\n",nop+2,nopa);
nopa=int(1.5*nopa); nopa=int(1.5*nopa);
@ -3171,13 +3464,17 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
} }
save_op_s=(t_save_op_s*)(&(save_op[nop]));
save_op[nop]=FDIV; save_op_s->operat=FDIV;
save_op_s->first=first->u_index;
save_op_s->lag=first->lag_index;
/*save_op[nop]=FDIV;
save_op[nop+1]=first->u_index; save_op[nop+1]=first->u_index;
save_op[nop+2]=0; save_op[nop+2]=0;
save_op[nop+3]=first->lag_index; save_op[nop+3]=first->lag_index;*/
} }
nop+=4; /*nop+=4;*/
nop+=2;
} }
else if (record_all) else if (record_all)
{ {
@ -3225,7 +3522,7 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
{ {
if (record) if (record)
{ {
if (nop+2>=nopa) if (nop+1>=nopa)
{ {
//mexPrintf("Error: out of save_op[%d] nopa=%d\n",nop+2,nopa); //mexPrintf("Error: out of save_op[%d] nopa=%d\n",nop+2,nopa);
nopa=int(1.5*nopa); nopa=int(1.5*nopa);
@ -3241,12 +3538,17 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
} }
save_op[nop]=FDIV; save_op_s=(t_save_op_s*)(&(save_op[nop]));
save_op_s->operat=FDIV;
save_op_s->first=b[pivj];
save_op_s->lag=0;
/*save_op[nop]=FDIV;
save_op[nop+1]=b[pivj]; save_op[nop+1]=b[pivj];
save_op[nop+2]=0; save_op[nop+2]=0;
save_op[nop+3]=0; save_op[nop+3]=0;*/
} }
nop+=4; /*nop+=4;*/
nop+=2;
} }
else if (record_all) else if (record_all)
{ {
@ -3295,7 +3597,7 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
{ {
if (record) if (record)
{ {
if (nop+2>=nopa) if (nop+1>=nopa)
{ {
//mexPrintf("Error: out of save_op[%d] nopa=%d\n",nop+2,nopa); //mexPrintf("Error: out of save_op[%d] nopa=%d\n",nop+2,nopa);
nopa=int(1.5*nopa); nopa=int(1.5*nopa);
@ -3311,12 +3613,17 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
} }
save_op[nop]=FLD; save_op_s=(t_save_op_s*)(&(save_op[nop]));
save_op_s->operat=FLD;
save_op_s->first=first->u_index;
save_op_s->lag=abs(first->lag_index);
/*save_op[nop]=FLD;
save_op[nop+1]=first->u_index; save_op[nop+1]=first->u_index;
save_op[nop+2]=0; save_op[nop+2]=0;
save_op[nop+3]=abs(first->lag_index); save_op[nop+3]=abs(first->lag_index);*/
} }
nop+=4; /*nop+=4;*/
nop+=2;
} }
else if (record_all) else if (record_all)
{ {
@ -3419,12 +3726,18 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
} }
save_op[nop]=FLESS; save_op_s=(t_save_op_s*)(&(save_op[nop]));
save_op_s->operat=FLESS;
save_op_s->first=tmp_u_count;
save_op_s->second=first_piv->u_index;
save_op_s->lag=max(first_piv->lag_index,abs(tmp_lag));
/*save_op[nop]=FLESS;
save_op[nop+1]=tmp_u_count; save_op[nop+1]=tmp_u_count;
save_op[nop+2]=first_piv->u_index; save_op[nop+2]=first_piv->u_index;
save_op[nop+3]=max(first_piv->lag_index,abs(tmp_lag)/*abs(lag)*/); save_op[nop+3]=max(first_piv->lag_index,abs(tmp_lag));*/
} }
nop+=4; /*nop+=4;*/
nop+=3;
} }
else if (record_all) else if (record_all)
{ {
@ -3539,12 +3852,18 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
} }
save_op[nop]=FSUB; save_op_s=(t_save_op_s*)(&(save_op[nop]));
save_op_s->operat=FSUB;
save_op_s->first=first_sub->u_index;
save_op_s->second=first_piv->u_index;
save_op_s->lag=max(abs(tmp_lag),first_piv->lag_index);
/*save_op[nop]=FSUB;
save_op[nop+1]=first_sub->u_index; save_op[nop+1]=first_sub->u_index;
save_op[nop+2]=first_piv->u_index; save_op[nop+2]=first_piv->u_index;
save_op[nop+3]=max(abs(/*first_sub->lag_index*/tmp_lag),first_piv->lag_index); save_op[nop+3]=max(abs(tmp_lag),first_piv->lag_index);*/
} }
nop+=4; /*nop+=4;*/
nop+=3;
} }
else if (record_all) else if (record_all)
{ {
@ -3630,12 +3949,18 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
} }
save_op[nop]=FSUB; save_op_s=(t_save_op_s*)(&(save_op[nop]));
save_op_s->operat=FSUB;
save_op_s->first=b[row];
save_op_s->second=b[pivj];
save_op_s->lag=abs(tmp_lag);
/*save_op[nop]=FSUB;
save_op[nop+1]=b[row]; save_op[nop+1]=b[row];
save_op[nop+2]=b[pivj]; save_op[nop+2]=b[pivj];
save_op[nop+3]=abs(tmp_lag); save_op[nop+3]=abs(tmp_lag);*/
} }
nop+=4; /*nop+=4;*/
nop+=3;
} }
else if (record_all) else if (record_all)
{ {
@ -3680,10 +4005,10 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
if (symbolic) if (symbolic)
{ {
//mexPrintf("record=%d, nop=%d, nopa=%d\n",record,nop,nopa); //mexPrintf("record=%d, nop=%d, nopa=%d\n",record,nop,nopa);
if (record && (nop==nopa)) if (record && (nop==nop1))
{ {
#ifdef PRINT_OUT #ifdef PRINT_OUT
mexPrintf("nop=%d, nopa=%d, record=%d save_opa=%x save_opaa=%x \n",nop, nopa, record, save_opa, save_opaa); mexPrintf("nop=%d, nop1=%d, record=%d save_opa=%x save_opaa=%x \n",nop, nop1, record, save_opa, save_opaa);
#endif #endif
//mexPrintf("save_opa=%x, save_opaa=%x\n",save_opa, save_opaa); //mexPrintf("save_opa=%x, save_opaa=%x\n",save_opa, save_opaa);
if (save_opa && save_opaa) if (save_opa && save_opaa)
@ -3714,23 +4039,27 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
save_opaa=NULL; save_opaa=NULL;
} }
#ifdef MEM_ALLOC_CHK #ifdef MEM_ALLOC_CHK
mexPrintf("save_opaa=(int*)mxMalloc(%d*sizeof(int))\n",nopa); mexPrintf("save_opaa=(int*)mxMalloc(%d*sizeof(int))\n",nop1);
#endif #endif
#ifdef N_MX_ALLOC #ifdef N_MX_ALLOC
save_opaa=malloc_std(nopa); save_opaa=malloc_std(nop1);
#else #else
save_opaa=(int*)mxMalloc(nopa*sizeof(int)); save_opaa=(int*)mxMalloc(nop1*sizeof(int));
#endif #endif
#ifdef MEM_ALLOC_CHK #ifdef MEM_ALLOC_CHK
mexPrintf("ok\n"); mexPrintf("ok\n");
#endif #endif
memcpy(save_opaa, save_opa, nopa*sizeof(int)); memcpy(save_opaa, save_opa, nop1*sizeof(int));
//mexPrintf("1 - end memcpy\n"); //mexPrintf("1 - end memcpy\n");
} }
//mexPrintf("okb\n"); //mexPrintf("okb\n");
if (save_opa) if (save_opa)
{ {
#ifdef N_MX_ALLOC
free(save_opa);
#else
mxFree(save_opa); mxFree(save_opa);
#endif
save_opa=NULL; save_opa=NULL;
} }
#ifdef MEM_ALLOC_CHK #ifdef MEM_ALLOC_CHK
@ -3749,7 +4078,7 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
} }
else else
{ {
if (nop==nopa) if (nop==nop1)
record=true; record=true;
else else
{ {
@ -3775,8 +4104,8 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
} }
//mexPrintf("nop=%d, nopa=%d, record=%d time=%f\n",nop, nopa, record,1000*(pctimer()-by_time_t0)); //mexPrintf("nop=%d, nopa=%d, record=%d time=%f\n",nop, nopa, record,1000*(pctimer()-by_time_t0));
} }
nopaa=nopa; nop2=nop1;
nopa=nop; nop1=nop;
} }
record=true; record=true;
} }