From 7d9b2a557b4e0a019e326267421efbe791090409 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Thu, 27 Jun 2019 17:00:12 +0200 Subject: [PATCH] perfect_foresight_problem MEX: number of threads is now configurable --- matlab/default_option_values.m | 1 + .../perfect_foresight_solver.m | 2 +- .../perfect_foresight_solver_core.m | 2 +- matlab/perfect-foresight-models/sim1.m | 2 +- .../perfect-foresight-models/solve_stacked_problem.m | 2 +- matlab/set_dynare_threads.m | 2 ++ .../perfect_foresight_problem.cc | 11 ++++++++--- 7 files changed, 15 insertions(+), 7 deletions(-) diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m index 68f7e31df..8cc2428f2 100644 --- a/matlab/default_option_values.m +++ b/matlab/default_option_values.m @@ -71,6 +71,7 @@ options_.huge_number = 1e7; % Default number of threads for parallelized mex files. options_.threads.kronecker.sparse_hessian_times_B_kronecker_C = num_procs; options_.threads.local_state_space_iteration_2 = 1; +options_.threads.perfect_foresight_problem = num_procs; % steady state options_.jacobian_flag = true; diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m index a67ad34a2..19c2ea4ae 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m @@ -191,7 +191,7 @@ if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_fores model_dynamic_g1_nz = str2func([M_.fname,'.dynamic_g1_nz']); [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz(); - residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll); + residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll, options_.threads.perfect_foresight_problem); if max(abs(residuals))< options_.dynatol.f oo_.deterministic_simulation.status = 1; diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m index 76a6e2d2f..e842e818f 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m @@ -140,7 +140,7 @@ if nargout>1 model_dynamic_g1_nz = str2func([M_.fname,'.dynamic_g1_nz']); [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz(); - residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll); + residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll, options_.threads.perfect_foresight_problem); end maxerror = max(max(abs(residuals))); end diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m index 1379ade92..58efa011b 100644 --- a/matlab/perfect-foresight-models/sim1.m +++ b/matlab/perfect-foresight-models/sim1.m @@ -67,7 +67,7 @@ h1 = clock; for iter = 1:options.simul.maxit h2 = clock; - [res, A] = perfect_foresight_problem(y, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll); + [res, A] = perfect_foresight_problem(y, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll, options.threads.perfect_foresight_problem); if options.endogenous_terminal_period && iter > 1 for it = 1:periods diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m index 797b23f9d..07d1804b3 100644 --- a/matlab/perfect-foresight-models/solve_stacked_problem.m +++ b/matlab/perfect-foresight-models/solve_stacked_problem.m @@ -56,7 +56,7 @@ else model_dynamic_g1_nz = str2func([M.fname,'.dynamic_g1_nz']); [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz(); - [y, check] = dynare_solve(@perfect_foresight_problem,z(:), options, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, options.periods, M.endo_nbr, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll); + [y, check] = dynare_solve(@perfect_foresight_problem,z(:), options, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, options.periods, M.endo_nbr, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll, options.threads.perfect_foresight_problem); end if all(imag(y)<.1*options.dynatol.x) diff --git a/matlab/set_dynare_threads.m b/matlab/set_dynare_threads.m index ca7618650..ca7827cb9 100644 --- a/matlab/set_dynare_threads.m +++ b/matlab/set_dynare_threads.m @@ -40,6 +40,8 @@ switch mexname options_.threads.kronecker.sparse_hessian_times_B_kronecker_C = n; case 'local_state_space_iteration_2' options_.threads.local_state_space_iteration_2 = n; + case 'perfect_foresight_problem' + options_.threads.perfect_foresight_problem = n; otherwise message = [ mexname ' is not a known parallel mex file.' ]; message_id = 'Dynare:Threads:UnknownParallelMex'; diff --git a/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc index f743c92b0..2fc9bfd80 100644 --- a/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc +++ b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc @@ -28,8 +28,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nlhs < 1 || nlhs > 2 || nrhs != 18) - mexErrMsgTxt("Must have 18 input arguments and 1 or 2 output arguments"); + if (nlhs < 1 || nlhs > 2 || nrhs != 19) + mexErrMsgTxt("Must have 19 input arguments and 1 or 2 output arguments"); bool compute_jacobian = nlhs == 2; // Give explicit names to input arguments @@ -51,6 +51,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) const mxArray *nzij_fwrd_mx = prhs[15]; const mxArray *has_external_function_mx = prhs[16]; const mxArray *use_dll_mx = prhs[17]; + const mxArray *num_threads_mx = prhs[18]; // Check input and map it to local variables if (!(mxIsChar(basename_mx) && mxGetM(basename_mx) == 1)) @@ -131,6 +132,10 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mexErrMsgTxt("use_dll should be a logical scalar"); bool use_dll = static_cast(mxGetScalar(use_dll_mx)); + if (!(mxIsScalar(num_threads_mx) && mxIsNumeric(num_threads_mx))) + mexErrMsgTxt("num_threads should be a numeric scalar"); + int num_threads = static_cast(mxGetScalar(num_threads_mx)); + // Allocate output matrices plhs[0] = mxCreateDoubleMatrix(periods*ny, 1, mxREAL); double *stacked_residual = mxGetPr(plhs[0]); @@ -178,7 +183,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* Parallelize the main loop, if use_dll and no external function (to avoid parallel calls to MATLAB) */ -#pragma omp parallel if (use_dll && !has_external_function) +#pragma omp parallel num_threads(num_threads) if (use_dll && !has_external_function) { // Allocate (thread-private) model evaluator (which allocates space for temporaries) std::unique_ptr m;