Commit 8f364376 authored by Ferhat Mihoubi's avatar Ferhat Mihoubi
Browse files

- Extends the algorithms compatible with bytecode to compute the steady-state...

- Extends the algorithms compatible with bytecode to compute the  steady-state (ticket #11). The new values for solve_algo are:
   * 0: fsolve
   * 1: solve1
   * 2, 4: solve1 + block decomposition
   * 3: csolve
   * 5: LU decomposition with UMFPack (method handling sparse matrix in Matlab)
   * 6: GMRES
   * 7: BiCGStab
   * 8: bytecode own solver (use Gaussian elimination + sparse matrix)

- Bytecode can now evaluate a specific block instead of the overall blocks (new bytecode's option 'Block = block_number')
parent 4b824ad8
......@@ -30,8 +30,8 @@ function steady_()
global M_ oo_ it_ options_
if options_.bytecode && ...
(options_.solve_algo ~= 1 && options_.solve_algo ~= 2 && options_.solve_algo ~= 3 && options_.solve_algo ~= 4 && options_.solve_algo ~= 5)
error('STEADY: for the moment, you must use solve_algo=5 with bytecode option')
(options_.solve_algo < 0 || options_.solve_algo > 8)
error('STEADY: for the moment, you must use solve_algo=1, 2, 3, 4, 5, 6, 7, 8 with bytecode option')
end
if ~options_.bytecode && options_.solve_algo == 5
error('STEADY: you can''t yet use solve_algo=5 without bytecode option')
......@@ -123,8 +123,28 @@ elseif options_.block && ~options_.bytecode
oo_.exo_det_steady_state], M_.params);
end
elseif options_.bytecode
[check, oo_.steady_state] = bytecode('static');
mexErrCheck('bytecode', check);
if options_.solve_algo > 4
[check, oo_.steady_state] = bytecode('static');
mexErrCheck('bytecode', check);
else
for b = 1:size(M_.blocksMFS,1)
n = size(M_.blocksMFS{b}, 1);
ss = oo_.steady_state;
if n ~= 0
[y, check] = dynare_solve('block_bytecode_mfs_steadystate', ...
ss(M_.blocksMFS{b}), ...
options_.jacobian_flag, b);
if check ~= 0
error(['STEADY: convergence problems in block ' int2str(b)])
end
oo_.steady_state(M_.blocksMFS{b}) = y;
else
[check, r, g1, oo_.steady_state] = feval('bytecode', oo_.steady_state, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params, 1, options_.solve_algo, 'static', 'evaluate', ['block = ' int2str(b)]);
end;
end
end
else
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
oo_.steady_state,...
......
This diff is collapsed.
......@@ -64,11 +64,11 @@ private:
/*string remove_white(string str);*/
string add_underscore_to_fpe(const string &str);
string get_variable(const SymbolType variable_type, const unsigned int variable_num);
string error_location(bool evaluate, bool steady_state);
string error_location(bool evaluate, bool steady_state, int size, int block_num);
void compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state);
string print_expression(it_code_type it_code, bool evaluate);
string print_expression(it_code_type it_code, bool evaluate, int size, int block_num, bool steady_state);
void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0, int block = -1);
int simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num,
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
double *T;
......@@ -77,7 +77,7 @@ private:
it_code_type it_code;
int Block_Count, Per_u_, Per_y_;
int it_, nb_row_x, nb_row_xd, maxit_, size_of_direction;
double *g1, *r;
double *g2, *g1, *r;
double solve_tolf;
bool GaussSeidel;
double *x, *params;
......@@ -93,11 +93,12 @@ public:
double *direction_arg, int y_size_arg, int nb_row_x_arg,
int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg);
bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks);
bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, int block, int &nb_blocks);
inline mxArray* get_jacob(int block_num) {return jacobian_block[block_num];};
inline mxArray* get_jacob_exo(int block_num) {return jacobian_exo_block[block_num];};
inline mxArray* get_jacob_exo_det(int block_num) {return jacobian_det_exo_block[block_num];};
inline mxArray* get_jacob_other_endo(int block_num) {return jacobian_other_endo_block[block_num];};
inline vector<double> get_residual() {return residual;};
};
#endif
......@@ -313,7 +313,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
for (j = 0; j < Size; j++)
IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
}
else if (stack_solve_algo >= 1 || stack_solve_algo <= 4)
else if (stack_solve_algo >= 0 || stack_solve_algo <= 4)
{
for (i = 0; i < u_count_init-Size; i++)
{
......@@ -330,7 +330,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
}
else
{
if ((stack_solve_algo == 5 && !steady_state) || (solve_algo == 5 && steady_state))
if ((stack_solve_algo == 5 && !steady_state) || (solve_algo == 8 && steady_state))
{
for (i = 0; i < u_count_init; i++)
{
......@@ -341,7 +341,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
IM_i[make_pair(make_pair(eq, var), lag)] = j;
}
}
else if ( ((stack_solve_algo >= 1 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 1 || solve_algo <= 4) && steady_state) )
else if ( ((stack_solve_algo >= 0 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 5 || solve_algo <= 7) && steady_state) )
{
for (i = 0; i < u_count_init; i++)
{
......@@ -2643,7 +2643,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
mexPrintf("-----------------------------------\n");
}
if (solve_algo == 5)
if (solve_algo == 8)
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
else
{
......@@ -2664,13 +2664,13 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
Init_Matlab_Sparse_Simple(Size, IM_i, A_m, b_m);
}
if (solve_algo == 1 || solve_algo == 4 || solve_algo == 0)
if (solve_algo == 5)
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
else if (solve_algo == 2)
else if (solve_algo == 6)
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_);
else if (solve_algo == 3)
else if (solve_algo == 7)
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_);
else if (solve_algo == 5)
else if (solve_algo == 8)
Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
return;
}
......
......@@ -134,6 +134,7 @@ private:
long int nop_all, nop1, nop2;
map<pair<pair<int, int>, int>, int> IM_i;
protected:
vector<double> residual;
int u_count_alloc, u_count_alloc_save;
double *u, *y, *ya;
vector<double*> jac;
......
......@@ -32,8 +32,6 @@ Get_Argument(const char *argv)
return f;
}
int
main(int nrhs, const char *prhs[])
#else
......@@ -51,34 +49,29 @@ Get_Argument(const mxArray *prhs)
mxFree(first_argument);
return f;
}
#endif
/* The gateway routine */
void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#endif
{
mxArray *M_, *oo_, *options_;
Get_Arguments_and_global_variables(int nrhs,
#ifndef DEBUG_EX
mxArray *block_structur = NULL;
const mxArray *prhs[],
#else
load_global((char*)prhs[1]);
const char *prhs[],
#endif
//ErrorHandlingException error_handling;
int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd;
int y_kmin = 0, y_kmax = 0, y_decal = 0, periods = 1;
double *direction;
bool steady_state = false;
bool evaluate = false;
bool block = false;
double *params = NULL;
double *yd = NULL, *xd = NULL;
int count_array_argument = 0;
int &count_array_argument,
double *yd[], unsigned int &row_y, unsigned int &col_y,
double *xd[], unsigned int &row_x, unsigned int &col_x,
double *params[], unsigned int &periods,
#ifndef DEBUG_EX
mxArray *block_structur[],
#endif
bool &steady_state, bool &evaluate, int &block,
mxArray *M_[], mxArray *oo_[], mxArray *options_[])
{
#ifdef DEBUG_EX
for (i = 2; i < nrhs; i++)
for (int i = 2; i < nrhs; i++)
#else
for (i = 0; i < nrhs; i++)
for (int i = 0; i < nrhs; i++)
#endif
{
#ifndef DEBUG_EX
......@@ -87,26 +80,27 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
switch (count_array_argument)
{
case 0:
yd = mxGetPr(prhs[i]);
*yd = mxGetPr(prhs[i]);
row_y = mxGetM(prhs[i]);
col_y = mxGetN(prhs[i]);
break;
case 1:
xd = mxGetPr(prhs[i]);
*xd = mxGetPr(prhs[i]);
row_x = mxGetM(prhs[i]);
col_x = mxGetN(prhs[i]);
break;
case 2:
params = mxGetPr(prhs[i]);
*params = mxGetPr(prhs[i]);
break;
case 3:
periods = mxGetScalar(prhs[i]);
break;
case 4:
block_structur = mxDuplicateArray(prhs[i]);
*block_structur = mxDuplicateArray(prhs[i]);
break;
default:
mexPrintf("Unknown argument\n");
//mexPrintf("Unknown argument count_array_argument=%d\n",count_array_argument);
break;
}
count_array_argument++;
}
......@@ -118,23 +112,38 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
steady_state = false;
else if (Get_Argument(prhs[i]) == "evaluate")
evaluate = true;
else if (Get_Argument(prhs[i]) == "block")
block = true;
else
{
ostringstream tmp;
tmp << " in main, unknown argument : " << Get_Argument(prhs[i]) << "\n";
throw FatalExceptionHandling(tmp.str());
int pos = Get_Argument(prhs[i]).find("block");
if (pos != (int)string::npos)
{
int pos1 = Get_Argument(prhs[i]).find("=", pos+5);
if (pos1 != (int)string::npos)
pos = pos1 + 1;
else
pos += 5;
block = atoi(Get_Argument(prhs[i]).substr(pos, string::npos).c_str())-1;
}
else
{
ostringstream tmp;
tmp << " in main, unknown argument : " << Get_Argument(prhs[i]) << "\n";
throw FatalExceptionHandling(tmp.str());
}
}
}
if (count_array_argument > 0 && count_array_argument < 4)
{
ostringstream tmp;
tmp << " in main, missing arguments. All the following arguments have to be indicated y, x, params, it_\n";
throw FatalExceptionHandling(tmp.str());
if (count_array_argument == 3 && steady_state)
periods = 1;
else
{
ostringstream tmp;
tmp << " in main, missing arguments. All the following arguments have to be indicated y, x, params, it_\n";
throw FatalExceptionHandling(tmp.str());
}
}
M_ = mexGetVariable("global", "M_");
*M_ = mexGetVariable("global", "M_");
if (M_ == NULL)
{
ostringstream tmp;
......@@ -142,20 +151,70 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
throw FatalExceptionHandling(tmp.str());
}
/* Gets variables and parameters from global workspace of Matlab */
oo_ = mexGetVariable("global", "oo_");
*oo_ = mexGetVariable("global", "oo_");
if (oo_ == NULL)
{
ostringstream tmp;
tmp << " in main, global variable not found: oo_\n";
throw FatalExceptionHandling(tmp.str());
}
options_ = mexGetVariable("global", "options_");
*options_ = mexGetVariable("global", "options_");
if (options_ == NULL)
{
ostringstream tmp;
tmp << " in main, global variable not found: options_\n";
throw FatalExceptionHandling(tmp.str());
}
}
#ifdef DEBUG_EX
int
main(int nrhs, const char *prhs[])
#else
/* The gateway routine */
void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#endif
{
mxArray *M_, *oo_, *options_;
#ifndef DEBUG_EX
mxArray *block_structur = NULL;
#else
int nlhs = 0;
char *plhs[1];
load_global((char*)prhs[1]);
#endif
//ErrorHandlingException error_handling;
unsigned int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd;
int y_kmin = 0, y_kmax = 0, y_decal = 0;
unsigned int periods = 1;
double *direction;
bool steady_state = false;
bool evaluate = false;
int block = -1;
double *params = NULL;
double *yd = NULL, *xd = NULL;
int count_array_argument = 0;
try
{
Get_Arguments_and_global_variables(nrhs, prhs, count_array_argument,
&yd, row_y, col_y,
&xd, row_x, col_x,
&params, periods,
#ifndef DEBUG_EX
&block_structur,
#endif
steady_state, evaluate, block,
&M_, &oo_, &options_);
}
catch (GeneralExceptionHandling &feh)
{
DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
}
if (!count_array_argument)
params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
......@@ -246,7 +305,9 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int y_size = row_y;
int nb_row_x = row_x;
clock_t t0 = clock();
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo);
string f(fname);
mxFree(fname);
int nb_blocks = 0;
......@@ -266,6 +327,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
if (!steady_state && !evaluate && no_error)
mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
#ifndef DEBUG_EX
bool dont_store_a_structure = false;
if (nlhs > 0)
{
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
......@@ -276,14 +338,35 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
pind[0] = 1;
if (nlhs > 1)
{
plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
pind = mxGetPr(plhs[1]);
if (evaluate)
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i]-ya[i];
if (block >= 0)
{
if (evaluate)
{
vector<double> residual = interprete.get_residual();
plhs[1] = mxCreateDoubleMatrix(residual.size()/col_y, col_y, mxREAL);
pind = mxGetPr(plhs[1]);
for (i = 0; i < residual.size(); i++)
pind[i] = residual[i];
}
else
{
plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
pind = mxGetPr(plhs[1]);
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i];
}
}
else
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i];
{
plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
pind = mxGetPr(plhs[1]);
if (evaluate)
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i]-ya[i];
else
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i];
}
if (nlhs > 2)
{
int jacob_field_number = 0, jacob_exo_field_number = 0, jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0;
......@@ -295,14 +378,20 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
jacob_exo_det_field_number=2;
jacob_other_endo_field_number=2;
mwSize dims[1] = {nb_blocks };
plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
mexPrintf("the structure has been created\n");
block_structur = plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
}
else if (!mxIsStruct(block_structur))
DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, the third output argument must be a structure\n");
if (block >=0 )
{
block_structur = plhs[2] = mxDuplicateArray(interprete.get_jacob(0));
//mexCallMATLAB(0,NULL, 1, &block_structur, "disp");
dont_store_a_structure = true;
}
else
DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, the third output argument must be a structure\n");
else
{
mexPrintf("Adding Fields\n");
jacob_field_number = mxAddField(block_structur, "jacob");
if (jacob_field_number == -1)
DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob to the structArray\n");
......@@ -316,14 +405,26 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
if (jacob_other_endo_field_number == -1)
DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob_other_endo to the structArray\n");
}
for (int i = 0; i < nb_blocks; i++)
if (!dont_store_a_structure)
{
for (int i = 0; i < nb_blocks; i++)
{
mxSetFieldByNumber(block_structur,i,jacob_field_number,interprete.get_jacob(i));
if (!evaluate)
{
mxSetFieldByNumber(block_structur,i,jacob_exo_field_number,interprete.get_jacob_exo(i));
mxSetFieldByNumber(block_structur,i,jacob_exo_det_field_number,interprete.get_jacob_exo_det(i));
mxSetFieldByNumber(block_structur,i,jacob_other_endo_field_number,interprete.get_jacob_other_endo(i));
}
}
}
if (nlhs > 3)
{
mxSetFieldByNumber(block_structur,i,jacob_field_number,interprete.get_jacob(i));
mxSetFieldByNumber(block_structur,i,jacob_exo_field_number,interprete.get_jacob_exo(i));
mxSetFieldByNumber(block_structur,i,jacob_exo_det_field_number,interprete.get_jacob_exo_det(i));
mxSetFieldByNumber(block_structur,i,jacob_other_endo_field_number,interprete.get_jacob_other_endo(i));
plhs[3] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
pind = mxGetPr(plhs[3]);
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i];
}
plhs[2] = block_structur;
}
}
}
......
/*
* Copyright (C) 2007-2010 Dynare Team
*
......@@ -1182,6 +1183,7 @@ class CodeLoad
private:
uint8_t *code;
unsigned int nb_blocks;
vector<unsigned int> begin_block;
public:
inline unsigned int
......@@ -1189,6 +1191,12 @@ public:
{
return nb_blocks;
};
unsigned int inline
get_begin_block(int block)
{
return begin_block[block];
}
inline void *
get_current_code()
{
......@@ -1445,6 +1453,7 @@ public:
code = fbegin_block->load(code);
begin_block.push_back(tags_liste.size());
tags_liste.push_back(make_pair(FBEGINBLOCK, fbegin_block));
nb_blocks++;
}
......
This diff is collapsed.
......@@ -39,9 +39,12 @@ private:
//! Temporary terms for the file containing parameters dervicatives
temporary_terms_t params_derivs_temporary_terms;
//! Temporary terms for block decomposed models
//! global temporary terms for block decomposed models
vector<vector<temporary_terms_t> > v_temporary_terms;
//! local temporary terms for block decomposed models
vector<vector<temporary_terms_t> > v_temporary_terms_local;
vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
typedef map< pair< int, pair< int, int> >, expr_t> first_chain_rule_derivatives_t;
......@@ -61,7 +64,7 @@ private:
void writeModelEquationsOrdered_M(const string &dynamic_basename) const;
//! Writes the code of the Block reordred structure of the model in virtual machine bytecode
void writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_t map_idx) const;
void writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_t map_idx, vector<map_idx_t> map_idx2) const;
//! Writes the code of the model in virtual machine bytecode
void writeModelEquationsCode(const string file_name, const string bin_basename, map_idx_t map_idx) const;
......@@ -76,15 +79,17 @@ private:
map_idx_t map_idx;
vector<map_idx_t> map_idx2;
//! sorts the temporary terms in the blocks order
void computeTemporaryTermsOrdered();
//! creates a mapping from the index of temporary terms to a natural index
void computeTemporaryTermsMapping();
void computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, map_idx_t &map_idx);
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx) const;
void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, map_idx_t &map_idx) const;
void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
//! Get the type corresponding to a derivation ID
virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment