From ec05451d1a3d4cea27da1f6cd91c736e0d171c07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Mon, 17 Jun 2019 16:17:49 +0200 Subject: [PATCH] Remove symmetric elements in 3rd model derivatives --- .../get_first_order_solution_params_deriv.m | 1 + matlab/unfold_g3.m | 45 +++++++++++++++++++ .../k_order_perturbation/k_ord_dynare.cc | 6 +-- preprocessor | 2 +- 4 files changed, 50 insertions(+), 4 deletions(-) create mode 100644 matlab/unfold_g3.m diff --git a/matlab/get_first_order_solution_params_deriv.m b/matlab/get_first_order_solution_params_deriv.m index 9c8d75df8..6cf4b0cc5 100644 --- a/matlab/get_first_order_solution_params_deriv.m +++ b/matlab/get_first_order_solution_params_deriv.m @@ -384,6 +384,7 @@ elseif (kronflag == 0 || kronflag == 1) d2g1_full = sparse(endo_nbr*yy0ex0_nbr, param_nbr*param_nbr); %initialize dyy0ex0 = sparse([dyy0; zeros(yy0ex0_nbr-yy0_nbr,param_nbr)]); %Jacobian (wrt model parameters) of steady state of dynamic (endogenous and auxiliary) and exogenous variables + g3 = unfold_g3(g3, yy0ex0_nbr); g3_tmp = reshape(g3,[endo_nbr*yy0ex0_nbr*yy0ex0_nbr yy0ex0_nbr]); d2g1_part4_left = sparse(endo_nbr*yy0ex0_nbr*yy0ex0_nbr,param_nbr); for j = 1:param_nbr diff --git a/matlab/unfold_g3.m b/matlab/unfold_g3.m new file mode 100644 index 000000000..e25e0eb11 --- /dev/null +++ b/matlab/unfold_g3.m @@ -0,0 +1,45 @@ +function g3_unfolded = unfold_g3(g3, ny) +% Given the 3rd order derivatives stored in a sparse matrix and without +% symmetric elements (as returned by the static/dynamic files) and the number +% of (static or dynamic )variables in the jacobian, returns +% an unfolded version of the same matrix (i.e. with symmetric elements). + +% Copyright (C) 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 . + +[i, j, v] = find(g3); + +i_unfolded = []; +j_unfolded = []; +v_unfolded = []; + +for k = 1:length(v) + l1 = rem(j(k)-1, ny); + j2 = floor((j(k)-1)/ny); + l2 = rem(j2, ny); + l3 = floor(j2/ny); + + p = unique(perms([l1 l2 l3]), 'rows'); + np = rows(p); + + i_unfolded = [i_unfolded; repmat(i(k), np, 1)]; + j_unfolded = [j_unfolded; 1 + p(:,1) + ny*(p(:,2) + ny*p(:,3))]; + v_unfolded = [v_unfolded; repmat(v(k), np, 1)]; +end + +g3_unfolded = sparse(i_unfolded, j_unfolded, v_unfolded, size(g3, 1), size(g3, 2)); + diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc index ac86a9ddb..e64ae03c9 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.cc +++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc @@ -130,9 +130,9 @@ KordpDynare::populateDerivativesContainer(int ord) i1 /= nJcols; } - if ((ord == 2 || ord == 3) && !s.isSorted()) - continue; // Skip symmetric elements (only needed at order 2 and 3) - else if (ord > 3) + if (ord == 2 && !s.isSorted()) + continue; // Skip symmetric elements (only needed at order 2) + else if (ord > 2) s.sort(); // For higher order, canonicalize the multi-index double x = g.get(i, 2); diff --git a/preprocessor b/preprocessor index 33c2f9b88..271a57980 160000 --- a/preprocessor +++ b/preprocessor @@ -1 +1 @@ -Subproject commit 33c2f9b88bb7c4477479fa6f553f001b299911ee +Subproject commit 271a57980847fb0b7c4d5c0f78e339f21ec847ad