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

Bug fix and efficiency change. 

- Even in models where there is only one endogenous variable in the
   LHS and where all the LHS are unique, it may be that because of the
   preprocessor transformations an auxiliary variable appears in more
   than one LHS. If diff(X) is on the LHS of an equation in the original
   model, the preprocessor will create an auxiliary variable AUX_DIFF
   which will appear in the the original equation, replacing diff(X),
   and in a new equation defining the auxiliary variable. In this case
   the, the Dulmage-Mendelsohn decomposition will associate AUX_DIFF
   with the original equation and X with the equation. This was
   problematic in the previous version of the algorithm, since it was
   assumed that each equation determines the LHS variable (here AUX_DIFF
   = X - X(-1) determines a RHS variable (X).

 - Changed the expression for evaluating an LHS variable under a log.

 - Improved efficiency by not evaluating the residuals of the model if
   not required for solving the current univariate block.
time-shift
Stéphane Adjemian (Charybdis) 2020-02-27 22:08:21 +01:00
parent 84e7d4bd3d
commit 67b18710e8
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
1 changed files with 48 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
matlab/dynare_solve.m
View File

@ -72,7 +72,10 @@ end
% Get first element of varargin if solve_algo ∈ {12,14} and rename varargin.
if ismember(options.solve_algo, [12, 14])
isloggedlhs = varargin{1};
arguments = varargin(2:end);
isauxdiffloggedrhs = varargin{2};
endo_names = varargin{3};
lhs = varargin{4};
arguments = varargin(5:end);
else
arguments = varargin;
end
@ -229,7 +232,7 @@ elseif ismember(options.solve_algo, [2, 12, 4, 14])
else
solver = @trust_region;
end
specializedunivariateblocks = ismember(options.solve_algo, [12, 14]);
specializedunivariateblocks = ismember(options.solve_algo, [12, 14]);
if ~jacobian_flag
fjac = zeros(nn,nn) ;
dh = max(abs(x), options.gstep(1)*ones(nn,1))*eps^(1/3);
@ -239,10 +242,13 @@ elseif ismember(options.solve_algo, [2, 12, 4, 14])
fjac(:,j) = (feval(f, xdh, arguments{:})-fvec)./dh(j) ;
end
end
[j1,j2,r] = dmperm(fjac);
[j1,j2,r,s] = dmperm(fjac);
JAC = abs(fjac(j1,j2))>tolf;
if options.debug
disp(['DYNARE_SOLVE (solve_algo=2|4|12|14): number of blocks = ' num2str(length(r))]);
end
l = 0;
fre = false;
for i=length(r)-1:-1:1
blocklength = r(i+1)-r(i);
if options.debug
@ -254,11 +260,39 @@ elseif ismember(options.solve_algo, [2, 12, 4, 14])
dprintf('DYNARE_SOLVE (solve_algo=2|4|12|14): solving block %u by evaluating RHS.', i);
end
if isequal(blocklength, 1)
z = feval(f, x, arguments{:});
if isloggedlhs(j1(j))
x(j2(j)) = x(j2(j))*exp(-z(j1(j)));
if i<length(r)-1
if fre || any(JAC(r(i), s(i)+(1:l)))
% Reevaluation of the residuals is required because the current RHS depends on
% variables that potentially have been updated previously.
z = feval(f, x, arguments{:});
l = 0;
fre = false;
end
else
x(j2(j)) = x(j2(j))-z(j1(j));
% First iteration requires the evaluation of the residuals.
z = feval(f, x, arguments{:});
end
l = l+1;
if isequal(lhs{j1(j)}, endo_names{j2(j)}) || isequal(lhs{j1(j)}, sprintf('log(%s)', endo_names{j2(j)}))
if isloggedlhs(j1(j))
x(j2(j)) = exp(log(x(j2(j)))-z(j1(j)));
else
x(j2(j)) = x(j2(j))-z(j1(j));
end
else
if options.debug
dprintf('LHS variable is not determined by RHS expression (%u).', j1(j))
dprintf('%s -> %s', lhs{j1(j)}, endo_names{j2(j)})
end
if ~isempty(regexp(lhs{j1(j)}, '\<AUX_DIFF_(\d*)\>', 'once'))
if isauxdiffloggedrhs(j1(j))
x(j2(j)) = exp(log(x(j2(j)))+z(j1(j)));
else
x(j2(j)) = x(j2(j))+z(j1(j));
end
else
error('Algorithm solve_algo=%u cannot be used with this nonlinear problem.', options.solve_algo)
end
end
continue
end
@ -268,18 +302,20 @@ elseif ismember(options.solve_algo, [2, 12, 4, 14])
end
end
[x, errorflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
options.gstep, ...
tolf, options.solve_tolx, ...
maxit, options.debug, arguments{:});
options.gstep, ...
tolf, options.solve_tolx, ...
maxit, options.debug, arguments{:});
fre = true;
if errorflag
return
end
end
fvec = feval(f, x, arguments{:});
if max(abs(fvec))>tolf
disp('Problem.')
[x, errorflag] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
options.gstep, tolf, options.solve_tolx, ...
maxit, options.debug, arguments{:});
options.gstep, tolf, options.solve_tolx, ...
maxit, options.debug, arguments{:});
end
elseif options.solve_algo==3
if jacobian_flag