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

Improve fallback code when mjdgges DLL is absent 

The new code relies on qz(..., 'real'), ordqz and ordeig, and returns a real
decomposition. The previous version was using Sims' qzdiv and returned a
complex decomposition.

As a consequence, we can drop options_.qzdiv, which was used to detect when
imaginary parts had to be dropped.

This code does not work on Octave for the time being, but this is acceptable
since it is only a fallback.
time-shift
Sébastien Villemot 2018-05-16 17:06:36 +02:00
parent c180777818
commit 59ad26a3c9
4 changed files with 31 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
matlab/dr_block.m
View File

@ -593,7 +593,7 @@ for i = 1:Size
else
[err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
end
if options_.aim_solver ~= 1 && options_.use_qzdiv
if options_.aim_solver ~= 1
% Necessary when using Sims' routines for QZ
ghx_other = real(ghx_other);
end
@ -650,7 +650,7 @@ for i = 1:Size
end
if options_.aim_solver ~= 1 && options_.use_qzdiv
if options_.aim_solver ~= 1
% Necessary when using Sims' routines for QZ
ghx = real(ghx);
ghu = real(ghu);

2
matlab/dyn_first_order_solver.m
View File

@ -297,7 +297,7 @@ end
dr.ghx = ghx;
dr.ghu = ghu;
if DynareOptions.aim_solver ~= 1 && DynareOptions.use_qzdiv
if DynareOptions.aim_solver ~= 1
% Necessary when using Sims' routines for QZ
dr.ghx = real(ghx);
dr.ghu = real(ghu);

9
matlab/global_initialization.m
View File

@ -334,15 +334,6 @@ options_.linear = 0;
options_.replic = 50;
options_.simul_replic = 1;
options_.drop = 100;
% if mjdgges.dll (or .mexw32 or ....) doesn't exist, matlab/qz is added to the path.
% There exists now qz/mjdgges.m that contains the calls to the old Sims code
% Hence, if mjdgges.m is visible exist(...)==2,
% this means that the DLL isn't avaiable and use_qzdiv is set to 1
if exist('mjdgges','file')==2
options_.use_qzdiv = 1;
else
options_.use_qzdiv = 0;
end
options_.aim_solver = 0; % i.e. by default do not use G.Anderson's AIM solver, use mjdgges instead
options_.k_order_solver=0; % by default do not use k_order_perturbation but mjdgges
options_.partial_information = 0;

73
matlab/qz/mjdgges.m
View File

@ -1,28 +1,25 @@
function [err,ss,tt,w,sdim,eigval,info] = mjdgges(e,d,qz_criterium, fake)
%function [err,ss,tt,w,sdim,eigval,info] = mjdgges(e,d,qz_criterium)
% QZ decomposition, Sims' codes are used.
function [err, ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhreshold)
%
% INPUTS
% e [double] real square (n*n) matrix.
% d [double] real square (n*n) matrix.
% qz_criterium [double] scalar (1+epsilon).
% zhreshold [double] ignored (in the DLL, this parameter is used for
% detecting eigenvalues too close to 0/0)
%
% OUTPUTS
% err [double] scalar: 1 indicates failure, 0 indicates success
% ss [complex] (n*n) matrix.
% tt [complex] (n*n) matrix.
% w [complex] (n*n) matrix.
% sdim [integer] scalar.
% eigval [complex] (n*1) vector.
% err [integer] scalar: 1 indicates failure, 0 indicates success
% ss [double] (n*n) quasi-triangular matrix.
% tt [double] (n*n) quasi-triangular matrix.
% zz [double] (n*n) orthogonal matrix.
% sdim [integer] scalar: number of stable eigenvalues.
% eigval [complex] (n*1) vector of generalized eigenvalues.
% info [integer] scalar.
%
% ALGORITHM
% Sims's qzdiv routine is used.
%
% SPECIAL REQUIREMENTS
% none.
% Copyright (C) 1996-2017 Dynare Team
% Copyright (C) 1996-2018 Dynare Team
%
% This file is part of Dynare.
%
@ -39,49 +36,35 @@ function [err,ss,tt,w,sdim,eigval,info] = mjdgges(e,d,qz_criterium, fake)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Check number of inputs and outputs.
if nargin>4 || nargin<2 || nargout>7 || nargout==0
error('MJDGGES: takes 2 or 3 input arguments and between 1 and 7 output arguments.')
if nargin > 5 || nargin < 2 || nargout > 7 || nargout == 0
error('MJDGGES: takes 2, 3 or 4 input arguments and between 1 and 7 output arguments.')
end
% Check the first two inputs.
[me,ne] = size(e);
[md,nd] = size(d);
if ( ~isreal(e) || ~isreal(d) || me~=ne || md~=nd || me~=nd)
% info should be negative in this case, see dgges.f.
if isoctave
error('Octave unsupported, since it does not have real qz, ordqz and ordeig')
end
[me, ne] = size(e);
[md, nd] = size(d);
if ~isreal(e) || ~isreal(d) || me ~= ne || md ~= nd || me ~= nd
error('MJDGGES requires two square real matrices of the same dimension.')
end
% Set default value of qz_criterium.
if nargin <3
if nargin < 3
qz_criterium = 1 + 1e-6;
end
info = 0;
% Initialization of the output arguments.
ss = zeros(ne,ne);
tt = zeros(ne,ne);
w = zeros(ne,ne);
sdim = 0;
eigval = zeros(ne,1);
% Computational part.
try
if isoctave()
% Force complex QZ factorization.
e = complex(e);
d = complex(d);
end
[ss,tt,qq,w,a,b,c] = qz(e, d);
warning_old_state = warning;
warning off;
[tt,ss,qq,w] = qzdiv(qz_criterium, tt, ss, qq, w);
eigval = diag(ss)./diag(tt);
warning(warning_old_state);
sdim = sum(abs(eigval) < qz_criterium);
[ss, tt, qq, zz] = qz(e, d, 'real');
eigval = ordeig(ss, tt);
select = abs(eigval) < qz_criterium;
sdim = sum(select);
[ss, tt, qq, zz] = ordqz(ss, tt, qq, zz, select);
eigval = ordeig(ss, tt);
catch
info = 1;% Not as precise as lapack's info!
info = 1; % Not as precise as lapack's info!
end
err = 0;
err = 0;