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

- The Temporary terms management with bytecode when the model is block decomposed and solve_algo<5 is now compatible with octave 

- Gets rid of warning message during the compilation of bytecode
time-shift
Ferhat Mihoubi 2011-01-14 19:22:29 +01:00
parent 96b6f1bf05
commit ce07223628
9 changed files with 51 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
matlab/block_bytecode_mfs_steadystate.m
View File

@ -1,4 +1,4 @@
function [r, g1, temporary_terms] = block_bytecode_mfs_steadystate(y, b, y_all, temporary_terms)
function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all)
% Wrapper around the *_static.m file, for use with dynare_solve,
% when block_mfs option is given to steady.
@ -20,7 +20,8 @@ function [r, g1, temporary_terms] = block_bytecode_mfs_steadystate(y, b, y_all,
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ oo_
global temporary_terms;
indx = M_.blocksMFS{b};
y_all(indx) = y;
x = [oo_.exo_steady_state; oo_.exo_det_steady_state];
eval(['[chk, r, g1] = bytecode( y_all, x, M_.params, 1, y_all, ''evaluate'', ''static'', ''block = ' int2str(b) ''', ''global_temporary_terms'');']);
eval(['[chk, r, g1, nulldev, temporary_terms] = bytecode( y_all, x, M_.params, 1, y_all, temporary_terms, ''evaluate'', ''static'', ''block = ' int2str(b) ''', ''global_temporary_terms'');']);

2
matlab/dynare_solve.m
View File

@ -57,7 +57,7 @@ if options_.solve_algo == 0
% Under Octave, use a wrapper, since fsolve() does not have a 4th arg
func2 = str2func(func);
func = @(x) func2(x, varargin{:});
% Octave do not converge when it starts from the solution
% The Octave version of fsolve do not converge when it starts from the solution
[fvec,fjac] = feval(func,x,varargin{:});
if max(abs(fvec)) >= options_.solve_tolf
[x,fval,exitval,output] = fsolve(func,x,options);

3
matlab/dynare_solve_block_or_bytecode.m
View File

@ -48,13 +48,12 @@ elseif options_.bytecode
mexErrCheck('bytecode', check);
info = check;
elseif options_.block
temporary_terms = [];
for b = 1:size(M_.blocksMFS,1)
n = size(M_.blocksMFS{b}, 1);
if n ~= 0
[y, check] = dynare_solve('block_bytecode_mfs_steadystate', ...
x(M_.blocksMFS{b}), ...
options_.jacobian_flag, b, x, temporary_terms);
options_.jacobian_flag, b, x);
if check ~= 0
error(['STEADY: convergence problems in block ' int2str(b)])
end

17
matlab/solve_one_boundary.m
View File

@ -323,11 +323,20 @@ for it_=start:incr:finish
end
while(flag1>0)
[L1, U1]=luinc(g1,luinc_tol);
phat = ya + U1 \ (L1 \ r);
res = g1 * phat;
if max(abs(res)) >= options_.solve_tolf
phat = ya - U1 \ (L1 \ r);
if(is_dynamic)
y(it_,y_index_eq) = phat;
else
y(y_index_eq) = phat;
end;
if(is_dynamic)
[r, y, g1, g2, g3] = feval(fname, y, x, params, it_, 0);
else
[r, y, g1] = feval(fname, y, x, params);
end;
if max(abs(r)) >= options_.solve_tolf
[dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
else
else
flag1 = 0;
dx = phat - ya;
end;

4
mex/sources/bytecode/ErrorHandling.hh
View File

@ -361,13 +361,13 @@ print_expression(it_code_type it_code, bool evaluate, int size, int block_num, b
double v1f, v2f, v3f =0.0;
bool go_on = true;
double ll;
ExpressionType equation_type ;
ExpressionType equation_type = TemporaryTerm;
unsigned int equation_num;
unsigned int dvar1, dvar2, dvar3;
int lag1, lag2, lag3;
size_t found;
double *jacob = NULL, *jacob_other_endo = NULL, *jacob_exo = NULL, *jacob_exo_det = NULL;
external_function_type function_type ;
external_function_type function_type = ExternalFunctionWithoutDerivative;
if (evaluate)
{

11
mex/sources/bytecode/Interpreter.cc
View File

@ -33,7 +33,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
int nb_row_x_arg, int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
int maxit_arg_, double solve_tolf_arg, int size_of_direction_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg,
string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
bool global_temporary_terms_arg, bool print_arg)
bool global_temporary_terms_arg, bool print_arg, mxArray* GlobalTemporaryTerms_arg)
{
params = params_arg;
y = y_arg;
@ -63,6 +63,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
solve_algo = solve_algo_arg;
global_temporary_terms = global_temporary_terms_arg;
print = print_arg;
GlobalTemporaryTerms = GlobalTemporaryTerms_arg;
}
double
@ -132,7 +133,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
double *jacob = NULL, *jacob_other_endo = NULL, *jacob_exo = NULL, *jacob_exo_det = NULL;
EQN_block = block_num;
stack<double> Stack;
external_function_type function_type;
external_function_type function_type = ExternalFunctionWithoutDerivative;
#ifdef DEBUG
mexPrintf("compute_block_time\n");
@ -2574,7 +2575,6 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
mxFree(T);
if (global_temporary_terms)
{
GlobalTemporaryTerms = mexGetVariable("caller", "temporary_terms");
if (GlobalTemporaryTerms == NULL)
{
mexPrintf("GlobalTemporaryTerms is NULL\n");mexEvalString("drawnow;");
@ -2585,6 +2585,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
}
else
T = (double *) mxMalloc(var*sizeof(double));
if (block >= 0)
it_code = code_liste.begin() + code.get_begin_block(block);
else
@ -2597,12 +2598,10 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
throw FatalExceptionHandling(tmp.str());
}
}
if (global_temporary_terms)
mexPutVariable("caller", "temporary_terms", GlobalTemporaryTerms);
mxFree(Init_Code->second);
nb_blocks = Block_Count+1;
if (T)
if (T and !global_temporary_terms)
mxFree(T);
return result;
}

3
mex/sources/bytecode/Interpreter.hh
View File

@ -81,13 +81,14 @@ 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 global_temporary_terms_arg, bool print_arg);
bool global_temporary_terms_arg, bool print_arg, mxArray* GlobalTemporaryTerms_arg);
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;};
inline mxArray* get_Temporary_Terms() {return GlobalTemporaryTerms;};
};
#endif

8
mex/sources/bytecode/SparseMatrix.cc
View File

@ -1806,7 +1806,7 @@ SparseMatrix::Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int S
B2_j[++B2_var] = B2_nze;
while(A3_var < Size)
A3_j[++A3_var] = A3_nze;
mxArray *d1;
mxArray *d1 = NULL;
vector<pair<mxArray*, mxArray*> > triangular_form;
int last_t = 0;
double sumc=0, C_sumc = 1000;
@ -2051,7 +2051,7 @@ SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double sl
void
SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray* x0_m, bool steady_state)
{
int n = mxGetM(A_m);
unsigned int n = mxGetM(A_m);
/*[L1, U1]=luinc(g1a,luinc_tol);*/
mxArray *lhs0[2];
mxArray *rhs0[2];
@ -2151,7 +2151,7 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double
{
double *res = mxGetPr(z);
if (is_two_boundaries)
for (int i = 0; i < n; i++)
for (unsigned int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
@ -2159,7 +2159,7 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double
y[eq] += slowc * yy;
}
else
for (int i = 0; i < n; i++)
for (unsigned int i = 0; i < n; i++)
{
int eq = index_vara[i];
double yy = - (res[i] + y[eq+it_*y_size]);

23
mex/sources/bytecode/bytecode.cc
View File

@ -67,7 +67,8 @@ Get_Arguments_and_global_variables(int nrhs,
#endif
bool &steady_state, bool &evaluate, int &block,
mxArray *M_[], mxArray *oo_[], mxArray *options_[], bool &global_temporary_terms,
bool &print)
bool &print,
mxArray* GlobalTemporaryTerms[])
{
#ifdef DEBUG_EX
for (int i = 2; i < nrhs; i++)
@ -99,6 +100,10 @@ Get_Arguments_and_global_variables(int nrhs,
case 4:
*block_structur = mxDuplicateArray(prhs[i]);
break;
case 5:
global_temporary_terms = true;
*GlobalTemporaryTerms = mxDuplicateArray(prhs[i]);
break;
default:
//mexPrintf("Unknown argument count_array_argument=%d\n",count_array_argument);
break;
@ -182,6 +187,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#endif
{
mxArray *M_, *oo_, *options_;
mxArray* GlobalTemporaryTerms;
#ifndef DEBUG_EX
mxArray *block_structur = NULL;
#else
@ -215,7 +221,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#endif
steady_state, evaluate, block,
&M_, &oo_, &options_, global_temporary_terms,
print);
print, &GlobalTemporaryTerms);
}
catch (GeneralExceptionHandling &feh)
{
@ -315,7 +321,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int nb_row_x = row_x;
clock_t t0 = clock();
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print);
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, global_temporary_terms, print, GlobalTemporaryTerms);
string f(fname);
mxFree(fname);
@ -427,7 +433,18 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
pind = mxGetPr(plhs[3]);
for (i = 0; i < row_y*col_y; i++)
pind[i] = y[i];
if (nlhs > 4)
{
mxArray* GlobalTemporaryTerms = interprete.get_Temporary_Terms();
unsigned int nb_temp_terms = mxGetM(GlobalTemporaryTerms);
plhs[4] = mxCreateDoubleMatrix(nb_temp_terms, 1, mxREAL);
pind = mxGetPr(plhs[4]);
double *tt = mxGetPr(GlobalTemporaryTerms);
for (i = 0; i < nb_temp_terms; i++)
pind[i] = tt[i];
}
}
}
}
}