diff --git a/matlab/perfect-foresight-models/perfect_foresight_problem.m b/matlab/perfect-foresight-models/perfect_foresight_problem.m
deleted file mode 100644
index 07220cdcef0609fbe995bab06946c7ab4a2f5349..0000000000000000000000000000000000000000
--- a/matlab/perfect-foresight-models/perfect_foresight_problem.m
+++ /dev/null
@@ -1,102 +0,0 @@
-function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, Y0, YT, ...
- exo_simul, params, steady_state, ...
- maximum_lag, T, ny, i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, ...
- i_cols_j, i_cols_0, i_cols_J0)
-
-% Computes the residuals and the Jacobian matrix for a perfect foresight problem over T periods.
-%
-% INPUTS
-% - y [double] N*1 array, terminal conditions for the endogenous variables
-% - dynamic_function [handle] function handle to _dynamic-file
-% - Y0 [double] N*1 array, initial conditions for the endogenous variables
-% - YT [double] N*1 array, terminal conditions for the endogenous variables
-% - exo_simul [double] nperiods*M_.exo_nbr matrix of exogenous variables (in declaration order)
-% for all simulation periods
-% - params [double] nparams*1 array, parameter values
-% - steady_state [double] endo_nbr*1 vector of steady state values
-% - maximum_lag [scalar] maximum lag present in the model
-% - T [scalar] number of simulation periods
-% - ny [scalar] number of endogenous variables
-% - i_cols [double] indices of variables appearing in M.lead_lag_incidence
-% and that need to be passed to _dynamic-file
-% - i_cols_J1 [double] indices of contemporaneous and forward looking variables
-% appearing in M.lead_lag_incidence
-% - i_cols_1 [double] indices of contemporaneous and forward looking variables in
-% M.lead_lag_incidence in dynamic Jacobian (relevant in first period)
-% - i_cols_T [double] columns of dynamic Jacobian related to contemporaneous and backward-looking
-% variables (relevant in last period)
-% - i_cols_j [double] indices of contemporaneous variables in M.lead_lag_incidence
-% in dynamic Jacobian (relevant in intermediate periods)
-% - i_cols_0 [double] indices of contemporaneous variables in M.lead_lag_incidence in dynamic
-% Jacobian (relevant in problems with periods=1)
-% - i_cols_J0 [double] indices of contemporaneous variables appearing in M.lead_lag_incidence (relevant in problems with periods=1)
-%
-% OUTPUTS
-% - residuals [double] (N*T)*1 array, residuals of the stacked problem
-% - JJacobian [double] (N*T)*(N*T) array, Jacobian of the stacked problem
-%
-% ALGORITHM
-% None
-%
-% SPECIAL REQUIREMENTS
-% None.
-
-% Copyright (C) 1996-2019 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-
-YY = [Y0; y; YT];
-
-residuals = zeros(T*ny,1);
-if nargout == 2
- iJacobian = cell(T,1);
-end
-
-i_rows = 1:ny;
-i_cols_J = i_cols;
-offset = 0;
-
-for it = maximum_lag+(1:T)
- if nargout == 1
- residuals(i_rows) = dynamic_function(YY(i_cols), exo_simul, params, steady_state, it);
- elseif nargout == 2
- [residuals(i_rows),jacobian] = dynamic_function(YY(i_cols), exo_simul, params, steady_state, it);
- if T==1 && it==maximum_lag+1
- [rows, cols, vals] = find(jacobian(:,i_cols_0));
- iJacobian{1} = [rows, i_cols_J0(cols), vals];
- elseif it == maximum_lag+1
- [rows,cols,vals] = find(jacobian(:,i_cols_1));
- iJacobian{1} = [offset+rows, i_cols_J1(cols), vals];
- elseif it == maximum_lag+T
- [rows,cols,vals] = find(jacobian(:,i_cols_T));
- iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals];
- else
- [rows,cols,vals] = find(jacobian(:,i_cols_j));
- iJacobian{it-maximum_lag} = [offset+rows, i_cols_J(cols), vals];
- i_cols_J = i_cols_J + ny;
- end
- offset = offset + ny;
- end
- i_rows = i_rows + ny;
- i_cols = i_cols + ny;
-end
-
-if nargout == 2
- iJacobian = cat(1,iJacobian{:});
- JJacobian = sparse(iJacobian(:,1), iJacobian(:,2), iJacobian(:,3), T*ny, T*ny);
-end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index 9816d65967663de18797d6efc1749f55aa6ce995..a67ad34a2bee5d14e73efe3ded6efdf6b28f38a5 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -176,26 +176,23 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
end
-if ~isreal(oo_.endo_simul(:)) % can only happen without bytecode
- y0 = real(oo_.endo_simul(:,1));
- yT = real(oo_.endo_simul(:,periods+2));
- yy = real(oo_.endo_simul(:,2:periods+1));
- illi = M_.lead_lag_incidence';
- [i_cols,~,i_cols_j] = find(illi(:));
- illi = illi(:,2:3);
- [i_cols_J1,~,i_cols_1] = find(illi(:));
- i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
- if periods==1
- i_cols_0 = nonzeros(M_.lead_lag_incidence(2,:)');
- i_cols_J0 = find(M_.lead_lag_incidence(2,:)');
+if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_foresight_problem DLL
+ if M_.maximum_lag > 0
+ y0 = real(oo_.endo_simul(:, M_.maximum_lag));
else
- i_cols_0 = [];
- i_cols_J0 = [];
+ y0 = NaN(ny, 1);
end
- residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ...
- oo_.exo_simul,M_.params,oo_.steady_state, ...
- M_.maximum_lag, periods, M_.endo_nbr, i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
+ if M_.maximum_lead > 0
+ yT = real(oo_.endo_simul(:, M_.maximum_lag+periods+1));
+ else
+ yT = NaN(ny, 1);
+ end
+ yy = real(oo_.endo_simul(:,M_.maximum_lag+(1:periods)));
+ 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);
+
if max(abs(residuals))< options_.dynatol.f
oo_.deterministic_simulation.status = 1;
oo_.endo_simul=real(oo_.endo_simul);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 24464e8f1996eb1affc9099eb983c41b745d5143..76a6e2d2f55352fd9441c2277b2b636158b1a073 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -120,31 +120,27 @@ else
end
if nargout>1
- y0 = oo_.endo_simul(:,1);
- yT = oo_.endo_simul(:,periods+2);
- yy = oo_.endo_simul(:,2:periods+1);
- illi = M_.lead_lag_incidence';
- [i_cols,~,i_cols_j] = find(illi(:));
- illi = illi(:,2:3);
- [i_cols_J1,~,i_cols_1] = find(illi(:));
- i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
- if periods==1
- i_cols_0 = nonzeros(M_.lead_lag_incidence(2,:)');
- i_cols_J0 = find(M_.lead_lag_incidence(2,:)');
- else
- i_cols_0 = [];
- i_cols_J0 = [];
- end
if options_.block && ~options_.bytecode
maxerror = oo_.deterministic_simulation.error;
else
if options_.bytecode
[chck, residuals, ~]= bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
else
- residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ...
- oo_.exo_simul,M_.params,oo_.steady_state, ...
- M_.maximum_lag, periods,M_.endo_nbr,i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
+ if M_.maximum_lag > 0
+ y0 = oo_.endo_simul(:, M_.maximum_lag);
+ else
+ y0 = NaN(ny, 1);
+ end
+ if M_.maximum_lead > 0
+ yT = oo_.endo_simul(:, M_.maximum_lag+periods+1);
+ else
+ yT = NaN(ny, 1);
+ end
+ yy = oo_.endo_simul(:,M_.maximum_lag+(1:periods));
+ 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);
end
maxerror = max(max(abs(residuals)));
end
diff --git a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
index 55290fc9f97c2e3a8e13af44c867725caa3eefd3..5677add41d9678b82ea0c2d240cf737f3c839b6c 100644
--- a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
+++ b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
@@ -70,8 +70,16 @@ elseif (options.solve_algo == 11)
end
end
-y0 = endogenousvariables(:,M.maximum_lag);
-yT = endogenousvariables(:,M.maximum_lag+periods+1);
+if M.maximum_lag > 0
+ y0 = endogenousvariables(:, M.maximum_lag);
+else
+ y0 = NaN(ny, 1);
+end
+if M.maximum_lead > 0
+ yT = endogenousvariables(:, M.maximum_lag+periods+1);
+else
+ yT = NaN(ny, 1);
+end
z = endogenousvariables(:,M.maximum_lag+(1:periods));
illi = M.lead_lag_incidence';
[i_cols,~,i_cols_j] = find(illi(:));
@@ -85,4 +93,4 @@ else
i_cols_0 = [];
i_cols_J0 = [];
end
-dynamicmodel = str2func([M.fname,'.dynamic']);
\ No newline at end of file
+dynamicmodel = str2func([M.fname,'.dynamic']);
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 1453875a7875427353961a512ea6032ed0bd3c2a..1379ade92c0b1c57333de4dea992d84e02d1c6ad 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -1,21 +1,17 @@
function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options)
-
-% Performs deterministic simulations with lead or lag on one period. Uses sparse matrices.
+% Performs deterministic simulations with lead or lag of one period, using
+% a basic Newton solver on sparse matrices.
+% Uses perfect_foresight_problem DLL to construct the stacked problem.
%
% INPUTS
-% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial guess).
+% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (initial condition + initial guess + terminal condition).
% - exogenousvariables [double] T*M array, paths for the exogenous variables.
% - steadystate [double] N*1 array, steady state for the endogenous variables.
% - M [struct] contains a description of the model.
% - options [struct] contains various options.
% OUTPUTS
-% - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
+% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model).
% - info [struct] contains informations about the results.
-% ALGORITHM
-% ...
-%
-% SPECIAL REQUIREMENTS
-% None.
% Copyright (C) 1996-2019 Dynare Team
%
@@ -36,38 +32,25 @@ function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariab
verbose = options.verbosity;
-endogenous_terminal_period = options.endogenous_terminal_period;
-vperiods = options.periods*ones(1,options.simul.maxit);
-azero = options.dynatol.f/1e7;
-
-lead_lag_incidence = M.lead_lag_incidence;
-
ny = M.endo_nbr;
+periods = options.periods;
+vperiods = periods*ones(1,options.simul.maxit);
-maximum_lag = M.maximum_lag;
-max_lag = M.maximum_endo_lag;
-
-nyp = nnz(lead_lag_incidence(1,:)) ;
-ny0 = nnz(lead_lag_incidence(2,:)) ;
-nyf = nnz(lead_lag_incidence(3,:)) ;
-
-nd = nyp+ny0+nyf;
-stop = 0 ;
+if M.maximum_lag > 0
+ y0 = endogenousvariables(:, M.maximum_lag);
+else
+ y0 = NaN(ny, 1);
+end
-periods = options.periods;
-params = M.params;
+if M.maximum_lead > 0
+ yT = endogenousvariables(:, M.maximum_lag+periods+1);
+else
+ yT = NaN(ny, 1);
+end
-i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)');
-i_cols_A1 = find(lead_lag_incidence(2:3,:)');
-i_cols_A1 = i_cols_A1(:);
-i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
-i_cols_0 = nonzeros(lead_lag_incidence(2,:)');
-i_cols_A0 = find(lead_lag_incidence(2,:)');
-i_cols_A0 = i_cols_A0(:);
-i_cols_j = (1:nd)';
-i_upd = maximum_lag*ny+(1:periods*ny);
+y = reshape(endogenousvariables(:, M.maximum_lag+(1:periods)), ny*periods, 1);
-Y = endogenousvariables(:);
+stop = false;
if verbose
skipline()
@@ -76,121 +59,78 @@ if verbose
skipline()
end
-model_dynamic = str2func([M.fname,'.dynamic']);
-
-res = zeros(periods*ny,1);
-
-o_periods = periods;
-
-if endogenous_terminal_period
- ZERO = zeros(length(i_upd),1);
-end
+model_dynamic_g1_nz = str2func([M.fname,'.dynamic_g1_nz']);
+[nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
-h1 = clock ;
-iA = zeros(periods*M.NNZDerivatives(1),3);
+h1 = clock;
for iter = 1:options.simul.maxit
- h2 = clock ;
- i_rows = (1:ny)';
- i_cols_A = find(lead_lag_incidence');
- i_cols_A = i_cols_A(:);
- i_cols = i_cols_A+(maximum_lag-1)*ny;
- m = 0;
- for it = (maximum_lag+1):(maximum_lag+periods)
- [d1,jacobian] = model_dynamic(Y(i_cols), exogenousvariables, params, steadystate,it);
- if periods==1 && it==maximum_lag+1
- [r,c,v] = find(jacobian(:,i_cols_0));
- iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A0(c(:)),v(:)];
- elseif it == maximum_lag+periods
- [r,c,v] = find(jacobian(:,i_cols_T));
- iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A(i_cols_T(c(:))),v(:)];
- elseif it == maximum_lag+1
- [r,c,v] = find(jacobian(:,i_cols_1));
- iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A1(c(:)),v(:)];
- else
- [r,c,v] = find(jacobian(:,i_cols_j));
- iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A(c(:)),v(:)];
- end
- m = m + length(v);
- res(i_rows) = d1;
- if endogenous_terminal_period && iter>1
- dr = max(abs(d1));
- if dr<azero
- vperiods(iter) = it;
- periods = it-maximum_lag+1;
+ 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);
+
+ if options.endogenous_terminal_period && iter > 1
+ for it = 1:periods
+ if max(abs(res((it-1)*ny+(1:ny)))) < options.dynatol.f/1e7
+ if it < periods
+ res = res(1:(it*ny));
+ A = A(1:(it*ny), 1:(it*ny));
+ yT = y(it*ny+(1:ny));
+ endogenousvariables(:, M.maximum_lag+((it+1):periods)) = reshape(y(it*ny+(1:(ny*(periods-it)))), ny, periods-it);
+ y = y(1:(it*ny));
+ periods = it;
+ end
break
end
end
- i_rows = i_rows + ny;
- i_cols = i_cols + ny;
- if it > maximum_lag+1
- i_cols_A = i_cols_A + ny;
- end
+ vperiods(iter) = periods;
end
+
err = max(abs(res));
if options.debug
fprintf('\nLargest absolute residual at iteration %d: %10.3f\n',iter,err);
- if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y))
+ if any(isnan(res)) || any(isinf(res)) || any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables)))
fprintf('\nWARNING: NaN or Inf detected in the residuals or endogenous variables.\n');
end
- if ~isreal(res) || ~isreal(Y)
- fprintf('\nWARNING: Imaginary parts detected in the residuals or endogenous variables.\n');
- end
skipline()
end
if verbose
- str = sprintf('Iter: %s,\t err. = %s, \t time = %s',num2str(iter),num2str(err), num2str(etime(clock,h2)));
- disp(str);
+ fprintf('Iter: %d,\t err. = %g,\t time = %g\n', iter, err, etime(clock,h2));
end
if err < options.dynatol.f
- stop = 1 ;
+ stop = true;
break
end
- iA = iA(1:m,:);
- A = sparse(iA(:,1),iA(:,2),iA(:,3),periods*ny,periods*ny);
- if endogenous_terminal_period && iter>1
- dy = ZERO;
- if options.simul.robust_lin_solve
- dy(1:i_rows(end)) = -lin_solve_robust( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)),verbose );
- else
- dy(1:i_rows(end)) = -lin_solve( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)), verbose );
- end
+ if options.simul.robust_lin_solve
+ dy = -lin_solve_robust(A, res, verbose);
else
- if options.simul.robust_lin_solve
- dy = -lin_solve_robust( A, res, verbose );
- else
- dy = -lin_solve( A, res, verbose );
- end
+ dy = -lin_solve(A, res, verbose);
end
- if any(~isreal(dy)) || any(isnan(dy)) || any(isinf(dy))
+ if any(isnan(dy)) || any(isinf(dy))
if verbose
display_critical_variables(reshape(dy,[ny periods])', M);
end
end
- Y(i_upd) = Y(i_upd) + dy;
+ y = y + dy;
end
-if endogenous_terminal_period
- err = evaluate_max_dynamic_residual(model_dynamic, Y, exogenousvariables, params, steadystate, o_periods, ny, max_lag, lead_lag_incidence);
- periods = o_periods;
-end
+endogenousvariables(:, M.maximum_lag+(1:periods)) = reshape(y, ny, periods);
+if options.endogenous_terminal_period
+ periods = options.periods;
+ err = evaluate_max_dynamic_residual(str2func([M.fname,'.dynamic']), endogenousvariables, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_endo_lag, M.lead_lag_incidence);
+end
if stop
- if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) || ~isreal(res) || ~isreal(Y)
+ if any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables)))
info.status = false;% NaN or Inf occurred
info.error = err;
info.iterations = iter;
info.periods = vperiods(1:iter);
- endogenousvariables = reshape(Y,ny,periods+maximum_lag+M.maximum_lead);
if verbose
skipline()
- disp(sprintf('Total time of simulation: %s.', num2str(etime(clock,h1))))
- if ~isreal(res) || ~isreal(Y)
- disp('Simulation terminated with imaginary parts in the residuals or endogenous variables.')
- else
- disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.')
- end
+ fprintf('Total time of simulation: %g.\n', etime(clock,h1))
+ disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.')
display_critical_variables(reshape(dy,[ny periods])', M);
disp('There is most likely something wrong with your model. Try model_diagnostics or another simulation method.')
printline(105)
@@ -198,19 +138,18 @@ if stop
else
if verbose
skipline();
- disp(sprintf('Total time of simulation: %s', num2str(etime(clock,h1))))
+ fprintf('Total time of simulation: %g.\n', etime(clock,h1))
printline(56)
end
info.status = true;% Convergency obtained.
info.error = err;
info.iterations = iter;
info.periods = vperiods(1:iter);
- endogenousvariables = reshape(Y,ny,periods+maximum_lag+M.maximum_lead);
end
elseif ~stop
if verbose
skipline();
- disp(sprintf('Total time of simulation: %s.', num2str(etime(clock,h1))))
+ fprintf('Total time of simulation: %g.\n', etime(clock,h1))
disp('Maximum number of iterations is reached (modify option maxit).')
printline(62)
end
@@ -224,21 +163,21 @@ if verbose
skipline();
end
-function x = lin_solve( A, b,verbose)
-if norm( b ) < sqrt( eps ) % then x = 0 is a solution
+function x = lin_solve(A, b, verbose)
+if norm(b) < sqrt(eps) % then x = 0 is a solution
x = 0;
return
end
x = A\b;
-x( ~isfinite( x ) ) = 0;
-relres = norm( b - A * x ) / norm( b );
+x(~isfinite(x)) = 0;
+relres = norm(b - A*x) / norm(b);
if relres > 1e-6 && verbose
- fprintf( 'WARNING : Failed to find a solution to the linear system.\n' );
+ fprintf('WARNING : Failed to find a solution to the linear system.\n');
end
-function [ x, flag, relres ] = lin_solve_robust( A, b , verbose)
-if norm( b ) < sqrt( eps ) % then x = 0 is a solution
+function [ x, flag, relres ] = lin_solve_robust(A, b ,verbose)
+if norm(b) < sqrt(eps) % then x = 0 is a solution
x = 0;
flag = 0;
relres = 0;
@@ -246,8 +185,8 @@ if norm( b ) < sqrt( eps ) % then x = 0 is a solution
end
x = A\b;
-x( ~isfinite( x ) ) = 0;
-[ x, flag, relres ] = bicgstab( A, b, [], [], [], [], x ); % returns immediately if x is a solution
+x(~isfinite(x)) = 0;
+[ x, flag, relres ] = bicgstab(A, b, [], [], [], [], x); % returns immediately if x is a solution
if flag == 0
return
end
@@ -255,38 +194,38 @@ end
disp( relres );
if verbose
- fprintf( 'Initial bicgstab failed, trying alternative start point.\n' );
+ fprintf('Initial bicgstab failed, trying alternative start point.\n');
end
old_x = x;
old_relres = relres;
-[ x, flag, relres ] = bicgstab( A, b );
+[ x, flag, relres ] = bicgstab(A, b);
if flag == 0
return
end
if verbose
- fprintf( 'Alternative start point also failed with bicgstab, trying gmres.\n' );
+ fprintf('Alternative start point also failed with bicgstab, trying gmres.\n');
end
if old_relres < relres
x = old_x;
end
-[ x, flag, relres ] = gmres( A, b, [], [], [], [], [], x );
+[ x, flag, relres ] = gmres(A, b, [], [], [], [], [], x);
if flag == 0
return
end
if verbose
- fprintf( 'Initial gmres failed, trying alternative start point.\n' );
+ fprintf('Initial gmres failed, trying alternative start point.\n');
end
old_x = x;
old_relres = relres;
-[ x, flag, relres ] = gmres( A, b );
+[ x, flag, relres ] = gmres(A, b);
if flag == 0
return
end
if verbose
- fprintf( 'Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n' );
+ fprintf('Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n');
end
if old_relres < relres
x = old_x;
@@ -294,15 +233,15 @@ if old_relres < relres
end
old_x = x;
old_relres = relres;
-x = pinv( full( A ) ) * b;
-relres = norm( b - A * x ) / norm( b );
+x = pinv(full(A)) * b;
+relres = norm(b - A*x) / norm(b);
if old_relres < relres
x = old_x;
relres = old_relres;
end
flag = relres > 1e-6;
if flag ~= 0 && verbose
- fprintf( 'WARNING : Failed to find a solution to the linear system\n' );
+ fprintf('WARNING : Failed to find a solution to the linear system\n');
end
function display_critical_variables(dyy, M)
@@ -321,10 +260,3 @@ if any(isinf(dyy))
fprintf('%s, ', endo_names{:}),
fprintf('\n'),
end
-if any(~isreal(dyy))
- indx = find(any(~isreal(dyy)));
- endo_names = M.endo_names(indx);
- disp('Last iteration provided complex number for the following variables:')
- fprintf('%s, ', endo_names{:}),
- fprintf('\n'),
-end
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
index 6290d41765d93bfa34b04655f36a3997a769ad30..797b23f9dadb35a2d9d29df9a28ac1bd09b79450 100644
--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -53,11 +53,10 @@ if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
eq_index);
else
- [y, check] = dynare_solve(@perfect_foresight_problem,z(:),options, ...
- dynamicmodel, y0, yT, ...
- exogenousvariables, M.params, steadystate, ...
- M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
+ 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);
end
if all(imag(y)<.1*options.dynatol.x)
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index 0f722508266b7d86cb0ff9bb037202b6aa9bfb65..8d555f1adcd7edb408258c8670464bf720635681 100644
--- a/mex/build/matlab/Makefile.am
+++ b/mex/build/matlab/Makefile.am
@@ -1,7 +1,7 @@
ACLOCAL_AMFLAGS = -I ../../../m4
if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations
+SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations perfect_foresight_problem
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if HAVE_MATIO
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index eb8f8b7b0332a3f268f15c618d6a090d2099cf3d..906aeee2e7a0ddab38ec70cd2f6436a43f761546 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -172,6 +172,7 @@ AC_CONFIG_FILES([Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile
sobol/Makefile
- local_state_space_iterations/Makefile])
+ local_state_space_iterations/Makefile
+ perfect_foresight_problem/Makefile])
AC_OUTPUT
diff --git a/mex/build/matlab/perfect_foresight_problem/Makefile.am b/mex/build/matlab/perfect_foresight_problem/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..d2ad7ff72c2f88a0e4c0a12cafcb7371769cc07e
--- /dev/null
+++ b/mex/build/matlab/perfect_foresight_problem/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../perfect_foresight_problem.am
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index 704a86bf0850ce10c54d4b04f739a754f2be5a12..4703f18c61225408b483b4fa76cbe4910dae4028 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/Makefile.am
@@ -1,7 +1,7 @@
ACLOCAL_AMFLAGS = -I ../../../m4
if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations
+SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations perfect_foresight_problem
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if HAVE_MATIO
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 478b9b22f739337250d347cb288b6bff1e6b11ba..7bb01f7047edbbb676450a54eee64dd7de1cacf6 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -1,6 +1,6 @@
dnl Process this file with autoconf to produce a configure script.
-dnl Copyright © 2009-2017 Dynare Team
+dnl Copyright © 2009-2019 Dynare Team
dnl
dnl This file is part of Dynare.
dnl
@@ -135,6 +135,7 @@ AC_CONFIG_FILES([Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile
sobol/Makefile
- local_state_space_iterations/Makefile])
+ local_state_space_iterations/Makefile
+ perfect_foresight_problem/Makefile])
AC_OUTPUT
diff --git a/mex/build/octave/perfect_foresight_problem/Makefile.am b/mex/build/octave/perfect_foresight_problem/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..65ceb979fa5bde26c66cefb3a0c29d8988100b15
--- /dev/null
+++ b/mex/build/octave/perfect_foresight_problem/Makefile.am
@@ -0,0 +1,3 @@
+EXEEXT = .mex
+include ../mex.am
+include ../../perfect_foresight_problem.am
diff --git a/mex/build/perfect_foresight_problem.am b/mex/build/perfect_foresight_problem.am
new file mode 100644
index 0000000000000000000000000000000000000000..e295c7e25257115040eca56a1c8541512a57d522
--- /dev/null
+++ b/mex/build/perfect_foresight_problem.am
@@ -0,0 +1,15 @@
+mex_PROGRAMS = perfect_foresight_problem
+
+TOPDIR = $(top_srcdir)/../../sources/perfect_foresight_problem
+
+nodist_perfect_foresight_problem_SOURCES = perfect_foresight_problem.cc DynamicModelCaller.cc
+
+perfect_foresight_problem_CPPFLAGS = $(AM_CPPFLAGS) -I$(TOPDIR)
+perfect_foresight_problem_CXXFLAGS = $(AM_CXXFLAGS) -fopenmp
+perfect_foresight_problem_LDADD = $(LIBADD_DLOPEN)
+
+BUILT_SOURCES = $(nodist_perfect_foresight_problem_SOURCES)
+CLEANFILES = $(nodist_perfect_foresight_problem_SOURCES)
+
+%.cc: $(TOPDIR)/%.cc
+ $(LN_S) -f $< $@
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 654100d3ab78f1c431c2bd4e53dbe04627610ad4..6b469e1cfd4c208fddc143b1d52d47e9491f8853 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -13,7 +13,8 @@ EXTRA_DIST = \
sobol \
local_state_space_iterations \
gensylv \
- dynare_simul_
+ dynare_simul_ \
+ perfect_foresight_problem
clean-local:
rm -rf `find mex/sources -name *.o`
diff --git a/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc b/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2a59e0937e06f75e3ded44126996960b87d49c32
--- /dev/null
+++ b/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
@@ -0,0 +1,215 @@
+/*
+ * Copyright © 2019 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#if defined(_WIN32) || defined(__CYGWIN32__)
+# ifndef NOMINMAX
+# define NOMINMAX // Do not define "min" and "max" macros
+# endif
+# include <windows.h>
+#else
+# include <dlfcn.h> // unix/linux DLL (.so) handling routines
+#endif
+
+#include <dynmex.h>
+
+#include "DynamicModelCaller.hh"
+
+#include <algorithm>
+
+std::string DynamicModelCaller::error_msg;
+
+void *DynamicModelDllCaller::dynamic_mex{nullptr};
+DynamicModelDllCaller::dynamic_tt_fct DynamicModelDllCaller::residual_tt_fct{nullptr}, DynamicModelDllCaller::g1_tt_fct{nullptr};
+DynamicModelDllCaller::dynamic_deriv_fct DynamicModelDllCaller::residual_fct{nullptr}, DynamicModelDllCaller::g1_fct{nullptr};
+
+void
+DynamicModelDllCaller::load_dll(const std::string &basename)
+{
+ // Load symbols from dynamic MEX
+ std::string mex_name;
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ mex_name = "./";
+#endif
+ mex_name += "+" + basename + "/dynamic" + MEXEXT;
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ dynamic_mex = dlopen(mex_name.c_str(), RTLD_NOW);
+#else
+ dynamic_mex = LoadLibrary(mex_name.c_str());
+#endif
+ if (!dynamic_mex)
+ mexErrMsgTxt("Can't load dynamic MEX file");
+
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ residual_tt_fct = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamic_mex, "dynamic_resid_tt"));
+ residual_fct = reinterpret_cast<dynamic_deriv_fct>(dlsym(dynamic_mex, "dynamic_resid"));
+ g1_tt_fct = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamic_mex, "dynamic_g1_tt"));
+ g1_fct = reinterpret_cast<dynamic_deriv_fct>(dlsym(dynamic_mex, "dynamic_g1"));
+#else
+ residual_tt_fct = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamic_mex, "dynamic_resid_tt"));
+ residual_fct = reinterpret_cast<dynamic_deriv_fct>(GetProcAddress(dynamic_mex, "dynamic_resid"));
+ g1_tt_fct = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamic_mex, "dynamic_g1_tt"));
+ g1_fct = reinterpret_cast<dynamic_deriv_fct>(GetProcAddress(dynamic_mex, "dynamic_g1"));
+#endif
+ if (!residual_tt_fct || !residual_fct || !g1_tt_fct || !g1_fct)
+ mexErrMsgTxt("Can't load functions in dynamic MEX file");
+}
+
+void
+DynamicModelDllCaller::unload_dll()
+{
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ dlclose(dynamic_mex);
+#else
+ FreeLibrary(dynamic_mex);
+#endif
+}
+
+DynamicModelDllCaller::DynamicModelDllCaller(size_t ntt, mwIndex nx, mwIndex ny, size_t ndynvars, const double *x_arg, size_t nb_row_x_arg, const double *params_arg, const double *steady_state_arg, bool compute_jacobian_arg) :
+ DynamicModelCaller{compute_jacobian_arg},
+ nb_row_x{nb_row_x_arg}, x{x_arg}, params{params_arg}, steady_state{steady_state_arg}
+{
+ tt = std::make_unique<double[]>(ntt);
+ y_p = std::make_unique<double[]>(ndynvars);
+ if (compute_jacobian)
+ jacobian_p = std::make_unique<double[]>((ndynvars+nx)*ny);
+}
+
+void
+DynamicModelDllCaller::eval(int it, double *resid)
+{
+ residual_tt_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get());
+ residual_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get(), resid);
+ if (compute_jacobian)
+ {
+ g1_tt_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get());
+ g1_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get(), jacobian_p.get());
+ }
+}
+
+DynamicModelMatlabCaller::DynamicModelMatlabCaller(std::string basename_arg, size_t ntt, size_t ndynvars, const mxArray *x_mx_arg, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, bool compute_jacobian_arg) :
+ DynamicModelCaller{compute_jacobian_arg},
+ basename{std::move(basename_arg)},
+ T_mx{mxCreateDoubleMatrix(ntt, 1, mxREAL)},
+ y_mx{mxCreateDoubleMatrix(ndynvars, 1, mxREAL)},
+ it_mx{mxCreateDoubleScalar(0)},
+ T_flag_mx{mxCreateLogicalScalar(false)},
+ jacobian_mx{nullptr},
+ x_mx{mxDuplicateArray(x_mx_arg)},
+ params_mx{mxDuplicateArray(params_mx_arg)},
+ steady_state_mx{mxDuplicateArray(steady_state_mx_arg)}
+{
+}
+
+DynamicModelMatlabCaller::~DynamicModelMatlabCaller()
+{
+ mxDestroyArray(T_mx);
+ mxDestroyArray(y_mx);
+ mxDestroyArray(it_mx);
+ mxDestroyArray(T_flag_mx);
+ if (jacobian_mx)
+ mxDestroyArray(jacobian_mx);
+ mxDestroyArray(x_mx);
+ mxDestroyArray(params_mx);
+ mxDestroyArray(steady_state_mx);
+}
+
+void
+DynamicModelMatlabCaller::eval(int it, double *resid)
+{
+ *mxGetPr(it_mx) = it + 1;
+
+ {
+ // Compute temporary terms
+ std::string funcname = basename + (compute_jacobian ? ".dynamic_g1_tt" : ".dynamic_resid_tt");
+ mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx };
+
+ mxArray *exception = mexCallMATLABWithTrap(1, plhs, 6, prhs, funcname.c_str());
+ if (exception)
+ error_msg = std::string{"An error occurred when calling "} + funcname;
+
+ mxDestroyArray(T_mx);
+ T_mx = plhs[0];
+ }
+
+ {
+ // Compute residuals
+ std::string funcname = basename + ".dynamic_resid";
+ mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx, T_flag_mx };
+
+ mxArray *exception = mexCallMATLABWithTrap(1, plhs, 7, prhs, funcname.c_str());
+ if (exception)
+ error_msg = std::string{"An error occurred when calling "} + funcname;
+
+ if (mxIsComplex(plhs[0]))
+ plhs[0] = cmplxToReal(plhs[0]);
+
+ std::copy_n(mxGetPr(plhs[0]), mxGetNumberOfElements(plhs[0]), resid);
+ mxDestroyArray(plhs[0]);
+ }
+
+ if (compute_jacobian)
+ {
+ // Compute Jacobian
+ std::string funcname = basename + ".dynamic_g1";
+ mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx, T_flag_mx };
+
+ mxArray *exception = mexCallMATLABWithTrap(1, plhs, 7, prhs, funcname.c_str());
+ if (exception)
+ error_msg = std::string{"An error occurred when calling "} + funcname;
+
+ if (jacobian_mx)
+ {
+ mxDestroyArray(jacobian_mx);
+ jacobian_mx = nullptr;
+ }
+
+ if (mxIsComplex(plhs[0]))
+ jacobian_mx = cmplxToReal(plhs[0]);
+ else
+ jacobian_mx = plhs[0];
+ }
+}
+
+mxArray *
+DynamicModelMatlabCaller::cmplxToReal(mxArray *cmplx_mx)
+{
+ mxArray *real_mx = mxCreateDoubleMatrix(mxGetM(cmplx_mx), mxGetN(cmplx_mx), mxREAL);
+
+#if MX_HAS_INTERLEAVED_COMPLEX
+ mxComplexDouble *cmplx = mxGetComplexDoubles(cmplx_mx);
+#else
+ double *cmplx_real = mxGetPr(cmplx_mx);
+ double *cmplx_imag = mxGetPi(cmplx_mx);
+#endif
+ double *real = mxGetPr(real_mx);
+
+ for (size_t i = 0; i < mxGetNumberOfElements(cmplx_mx); i++)
+#if MX_HAS_INTERLEAVED_COMPLEX
+ if (cmplx[i].imag == 0.0)
+ real[i] = cmplx[i].real;
+#else
+ if (cmplx_imag[i] == 0.0)
+ real[i] = cmplx_real[i];
+#endif
+ else
+ real[i] = std::numeric_limits<double>::quiet_NaN();
+
+ mxDestroyArray(cmplx_mx);
+ return real_mx;
+}
diff --git a/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh b/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8c64938fdb6faf8a420cb62dfa9b5d5981a237c0
--- /dev/null
+++ b/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh
@@ -0,0 +1,83 @@
+/*
+ * Copyright © 2019 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <string>
+#include <memory>
+#include <limits>
+
+class DynamicModelCaller
+{
+public:
+ const bool compute_jacobian;
+
+ // Used to store error messages (as exceptions cannot cross the OpenMP boundary)
+ static std::string error_msg;
+
+ DynamicModelCaller(bool compute_jacobian_arg) : compute_jacobian{compute_jacobian_arg} {};
+ virtual ~DynamicModelCaller() = default;
+ virtual double &y(size_t i) const = 0;
+ virtual double jacobian(size_t i) const = 0;
+ virtual void eval(int it, double *resid) = 0;
+};
+
+class DynamicModelDllCaller : public DynamicModelCaller
+{
+private:
+ static void *dynamic_mex;
+ using dynamic_tt_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T);
+ using dynamic_deriv_fct = void (*) (const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *deriv);
+ static dynamic_tt_fct residual_tt_fct, g1_tt_fct;
+ static dynamic_deriv_fct residual_fct, g1_fct;
+ size_t nb_row_x;
+ const double *x, *params, *steady_state;
+ std::unique_ptr<double[]> tt, y_p, jacobian_p;
+
+public:
+ DynamicModelDllCaller(size_t ntt, mwIndex nx, mwIndex ny, size_t ndynvars, const double *x_arg, size_t nb_row_x_arg, const double *params_arg, const double *steady_state_arg, bool compute_jacobian_arg);
+ virtual ~DynamicModelDllCaller() = default;
+ double &y(size_t i) const override { return y_p[i]; };
+ double jacobian(size_t i) const override { return jacobian_p[i]; };
+ void eval(int it, double *resid) override;
+ static void load_dll(const std::string &basename);
+ static void unload_dll();
+};
+
+class DynamicModelMatlabCaller : public DynamicModelCaller
+{
+private:
+ std::string basename;
+ mxArray *T_mx, *y_mx, *it_mx, *T_flag_mx, *jacobian_mx, *x_mx, *params_mx, *steady_state_mx;
+ /* Given a complex matrix, returns a real matrix of same size.
+ Real elements of the original matrix are copied as-is to the new one.
+ Complex elements are replaced by NaNs.
+ Destroys the original matrix. */
+ static mxArray *cmplxToReal(mxArray *m);
+public:
+ DynamicModelMatlabCaller(std::string basename_arg, size_t ntt, size_t ndynvars, const mxArray *x_mx_arg, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, bool compute_jacobian_arg);
+ ~DynamicModelMatlabCaller() override;
+ double &y(size_t i) const override { return mxGetPr(y_mx)[i]; };
+ double jacobian(size_t i) const override { return jacobian_mx ? mxGetPr(jacobian_mx)[i] : std::numeric_limits<double>::quiet_NaN(); };
+ void eval(int it, double *resid) override;
+ class Exception {
+ public:
+ const std::string msg;
+ Exception(std::string msg_arg) : msg{std::move(msg_arg)} {};
+ };
+};
+
diff --git a/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f743c92b023fe6b90883e1d7c79fb7bd442d9cf7
--- /dev/null
+++ b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
@@ -0,0 +1,277 @@
+/*
+ * Copyright © 2019 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <string>
+#include <memory>
+#include <algorithm>
+
+#include <dynmex.h>
+
+#include "DynamicModelCaller.hh"
+
+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");
+ bool compute_jacobian = nlhs == 2;
+
+ // Give explicit names to input arguments
+ const mxArray *y_mx = prhs[0];
+ const mxArray *basename_mx = prhs[1];
+ const mxArray *ntt_mx = prhs[2];
+ const mxArray *y0_mx = prhs[3];
+ const mxArray *yT_mx = prhs[4];
+ const mxArray *exo_path_mx = prhs[5];
+ const mxArray *params_mx = prhs[6];
+ const mxArray *steady_state_mx = prhs[7];
+ const mxArray *periods_mx = prhs[8];
+ const mxArray *ny_mx = prhs[9];
+ const mxArray *maximum_lag_mx = prhs[10];
+ const mxArray *maximum_endo_lag_mx = prhs[11];
+ const mxArray *lead_lag_incidence_mx = prhs[12];
+ const mxArray *nzij_pred_mx = prhs[13];
+ const mxArray *nzij_current_mx = prhs[14];
+ const mxArray *nzij_fwrd_mx = prhs[15];
+ const mxArray *has_external_function_mx = prhs[16];
+ const mxArray *use_dll_mx = prhs[17];
+
+ // Check input and map it to local variables
+ if (!(mxIsChar(basename_mx) && mxGetM(basename_mx) == 1))
+ mexErrMsgTxt("basename should be a character string");
+ char *basename = mxArrayToString(basename_mx);
+
+ if (!(mxIsScalar(ntt_mx) && mxIsNumeric(ntt_mx)))
+ mexErrMsgTxt("ntt should be a numeric scalar");
+ size_t ntt = mxGetScalar(ntt_mx);
+
+ if (!(mxIsScalar(periods_mx) && mxIsNumeric(periods_mx)))
+ mexErrMsgTxt("periods should be a numeric scalar");
+ mwIndex periods = static_cast<mwIndex>(mxGetScalar(periods_mx));
+
+ if (!(mxIsScalar(ny_mx) && mxIsNumeric(ny_mx)))
+ mexErrMsgTxt("ny should be a numeric scalar");
+ mwIndex ny = static_cast<mwIndex>(mxGetScalar(ny_mx));
+
+ if (!(mxIsScalar(maximum_lag_mx) && mxIsNumeric(maximum_lag_mx)))
+ mexErrMsgTxt("maximum_lag should be a numeric scalar");
+ mwIndex maximum_lag = static_cast<mwIndex>(mxGetScalar(maximum_lag_mx));
+
+ if (!(mxIsScalar(maximum_endo_lag_mx) && mxIsNumeric(maximum_endo_lag_mx)))
+ mexErrMsgTxt("maximum_endo_lag should be a numeric scalar");
+ mwIndex maximum_endo_lag = static_cast<mwIndex>(mxGetScalar(maximum_endo_lag_mx));
+
+ if (!(mxIsDouble(y_mx) && mxGetM(y_mx) == static_cast<size_t>(ny*periods) && mxGetN(y_mx) == 1))
+ mexErrMsgTxt("y should be a double precision column-vector of ny*periods elements");
+ const double *y = mxGetPr(y_mx);
+
+ if (!(mxIsDouble(y0_mx) && mxGetM(y0_mx) == static_cast<size_t>(ny) && mxGetN(y0_mx) == 1))
+ mexErrMsgTxt("y0 should be a double precision column-vector of ny elements");
+ const double *y0 = mxGetPr(y0_mx);
+
+ if (!(mxIsDouble(yT_mx) && mxGetM(yT_mx) == static_cast<size_t>(ny) && mxGetN(yT_mx) == 1))
+ mexErrMsgTxt("yT should be a double precision column-vector of ny elements");
+ const double *yT = mxGetPr(yT_mx);
+
+ if (!(mxIsDouble(exo_path_mx) && mxGetM(exo_path_mx) >= static_cast<size_t>(periods+maximum_lag)))
+ mexErrMsgTxt("exo_path should be a double precision matrix with at least periods+maximum_lag rows");
+ mwIndex nx = static_cast<mwIndex>(mxGetN(exo_path_mx));
+ size_t nb_row_x = mxGetM(exo_path_mx);
+ const double *exo_path = mxGetPr(exo_path_mx);
+
+ if (!(mxIsDouble(params_mx) && mxGetN(params_mx) == 1))
+ mexErrMsgTxt("params should be a double precision column-vector");
+ const double *params = mxGetPr(params_mx);
+
+ if (!(mxIsDouble(steady_state_mx) && mxGetN(steady_state_mx) == 1))
+ mexErrMsgTxt("steady_state should be a double precision column-vector");
+ const double *steady_state = mxGetPr(steady_state_mx);
+
+ if (!(mxIsDouble(lead_lag_incidence_mx) && mxGetM(lead_lag_incidence_mx) == static_cast<size_t>(2+maximum_endo_lag)
+ && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(ny)))
+ mexErrMsgTxt("lead_lag_incidence should be a double precision matrix with 2+maximum_endo_lag rows and endo_nbr columns");
+ const double *lead_lag_incidence = mxGetPr(lead_lag_incidence_mx);
+
+ if (!(mxIsInt32(nzij_pred_mx) && mxGetN(nzij_pred_mx) == 2))
+ mexErrMsgTxt("nzij_pred should be an int32 matrix with 2 columns");
+ size_t nnz_pred = mxGetM(nzij_pred_mx);
+ const int32_T *nzij_pred = static_cast<const int32_T *>(mxGetData(nzij_pred_mx));
+
+ if (!(mxIsInt32(nzij_current_mx) && mxGetN(nzij_current_mx) == 2))
+ mexErrMsgTxt("nzij_current should be an int32 matrix with 2 columns");
+ size_t nnz_current = mxGetM(nzij_current_mx);
+ const int32_T *nzij_current = static_cast<const int32_T *>(mxGetData(nzij_current_mx));
+
+ if (!(mxIsInt32(nzij_fwrd_mx) && mxGetN(nzij_fwrd_mx) == 2))
+ mexErrMsgTxt("nzij_fwrd should be an int32 matrix with 2 columns");
+ size_t nnz_fwrd = mxGetM(nzij_fwrd_mx);
+ const int32_T *nzij_fwrd = static_cast<const int32_T *>(mxGetData(nzij_fwrd_mx));
+
+ if (!(mxIsLogicalScalar(has_external_function_mx)))
+ mexErrMsgTxt("has_external_function should be a logical scalar");
+ bool has_external_function = static_cast<bool>(mxGetScalar(has_external_function_mx));
+
+ if (!(mxIsLogicalScalar(use_dll_mx)))
+ mexErrMsgTxt("use_dll should be a logical scalar");
+ bool use_dll = static_cast<bool>(mxGetScalar(use_dll_mx));
+
+ // Allocate output matrices
+ plhs[0] = mxCreateDoubleMatrix(periods*ny, 1, mxREAL);
+ double *stacked_residual = mxGetPr(plhs[0]);
+
+ mwIndex nzmax = periods*nnz_current+(periods-1)*(nnz_pred+nnz_fwrd);
+
+ double *stacked_jacobian = nullptr;
+ mwIndex *ir = nullptr, *jc = nullptr;
+ if (compute_jacobian)
+ {
+ plhs[1] = mxCreateSparse(periods*ny, periods*ny, nzmax, mxREAL);
+ stacked_jacobian = mxGetPr(plhs[1]);
+ ir = mxGetIr(plhs[1]);
+ jc = mxGetJc(plhs[1]);
+
+ /* Create the index vectors (IR, JC) of the sparse stacked jacobian. This
+ makes parallelization across periods possible when evaluating the model,
+ since all indices are known ex ante. */
+ mwIndex k = 0;
+ jc[0] = 0;
+ for (mwIndex T = 0; T < periods; T++)
+ {
+ size_t row_pred = 0, row_current = 0, row_fwrd = 0;
+ for (int32_T j = 0; j < static_cast<int32_T>(ny); j++)
+ {
+ if (T != 0)
+ while (row_fwrd < nnz_fwrd && nzij_fwrd[row_fwrd+nnz_fwrd]-1 == j)
+ ir[k++] = (T-1)*ny + nzij_fwrd[row_fwrd++]-1;
+ while (row_current < nnz_current && nzij_current[row_current+nnz_current]-1 == j)
+ ir[k++] = T*ny + nzij_current[row_current++]-1;
+ if (T != periods-1)
+ while (row_pred < nnz_pred && nzij_pred[row_pred+nnz_pred]-1 == j)
+ ir[k++] = (T+1)*ny + nzij_pred[row_pred++]-1;
+ jc[T*ny+j+1] = k;
+ }
+ }
+ }
+
+ size_t ndynvars = static_cast<size_t>(*std::max_element(lead_lag_incidence, lead_lag_incidence+(maximum_endo_lag+2)*ny));
+
+ if (use_dll)
+ DynamicModelDllCaller::load_dll(basename);
+
+ DynamicModelCaller::error_msg.clear();
+
+ /* 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)
+ {
+ // Allocate (thread-private) model evaluator (which allocates space for temporaries)
+ std::unique_ptr<DynamicModelCaller> m;
+ if (use_dll)
+ m = std::make_unique<DynamicModelDllCaller>(ntt, nx, ny, ndynvars, exo_path, nb_row_x, params, steady_state, compute_jacobian);
+ else
+ m = std::make_unique<DynamicModelMatlabCaller>(basename, ntt, ndynvars, exo_path_mx, params_mx, steady_state_mx, compute_jacobian);
+
+ // Main computing loop
+#pragma omp for
+ for (mwIndex T = 0; T < periods; T++)
+ {
+ // Fill vector of dynamic variables
+ for (mwIndex j = 0; j < maximum_endo_lag+2; j++)
+ for (mwIndex i = 0; i < ny; i++)
+ {
+ int idx = static_cast<int>(lead_lag_incidence[j+i*(2+maximum_endo_lag)])-1;
+ if (idx != -1)
+ {
+ if (T+j == maximum_endo_lag-1)
+ m->y(idx) = y0[i];
+ else if (T+j == maximum_endo_lag+periods)
+ m->y(idx) = yT[i];
+ else
+ m->y(idx) = y[i+(T+j-maximum_endo_lag)*ny];
+ }
+ }
+
+ // Compute the residual and Jacobian, and fill the stacked residual
+ m->eval(T+maximum_lag, stacked_residual+T*ny);
+
+ if (compute_jacobian)
+ {
+ // Fill the stacked jacobian
+ for (mwIndex col = T > maximum_endo_lag ? (T-maximum_endo_lag)*ny : 0; // We can't use std::max() here, because mwIndex is unsigned under MATLAB
+ col < std::min(periods*ny, (T+2)*ny); col++)
+ {
+ mwIndex k = jc[col];
+ while (k < jc[col+1])
+ {
+ if (ir[k] < T*ny)
+ {
+ k++;
+ continue;
+ }
+ if (ir[k] >= (T+1)*ny)
+ break;
+
+ mwIndex eq = ir[k]-T*ny;
+ mwIndex lli_row = col/ny-(T-maximum_endo_lag); // 0, 1 or 2
+ mwIndex lli_col = col%ny;
+ mwIndex dynvar = static_cast<mwIndex>(lead_lag_incidence[lli_row+lli_col*(2+maximum_endo_lag)])-1;
+ stacked_jacobian[k] = m->jacobian(eq+dynvar*ny);
+ k++;
+ }
+ }
+ }
+ }
+ }
+
+ /* Mimic a try/catch using a global string, since exceptions are not allowed
+ to cross OpenMP boundary */
+ if (!DynamicModelCaller::error_msg.empty())
+ mexErrMsgTxt(DynamicModelCaller::error_msg.c_str());
+
+ if (compute_jacobian)
+ {
+ /* The stacked jacobian so far constructed can still reference some zero
+ elements, because some expressions that are symbolically different from
+ zero may still evaluate to zero for some endogenous/parameter values. The
+ following code further compresses the Jacobian by removing the references
+ to those extra zeros. This makes a significant speed difference when
+ inversing the Jacobian for some large models. */
+ mwIndex k_orig = 0, k_new = 0;
+ for (mwIndex col = 0; col < periods*ny; col++)
+ {
+ while (k_orig < jc[col+1])
+ {
+ if (stacked_jacobian[k_orig] != 0.0)
+ {
+ if (k_new != k_orig)
+ {
+ stacked_jacobian[k_new] = stacked_jacobian[k_orig];
+ ir[k_new] = ir[k_orig];
+ }
+ k_new++;
+ }
+ k_orig++;
+ }
+ jc[col+1] = k_new;
+ }
+ }
+
+ if (use_dll)
+ DynamicModelDllCaller::unload_dll();
+}