1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
function [A,B,Q,Z] = qzdivct(stake,A,B,Q,Z)
%function [A,B,Q,Z] = qzdivct(stake,A,B,Q,Z)
%
% Takes U.T. matrices A, B, orthonormal matrices Q,Z, rearranges them
% so that all cases of real(B(i,i)/A(i,i))>stake are in lower right
% corner, while preserving U.T. and orthonormal properties and Q'AZ' and
% Q'BZ'. abs(A(i,i))<1e-11 is interpreted as a zero and as generating
% an infinitely positive real part of the ratio. All i's for which this
% criterion are satisfied are grouped together in the lower right corner
% of the lower right corner, with the non-zero roots above them. This
% version differs from
% qzdiv in that it works on the real part's value, as is appropriate for
% continuous time models, instead of on the absolute value, as is
% appropriate for discrete time models.
%
% Copyright (C) 1997-2012 Tao Zha
%
% This 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.
%
% It 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.
%
% If you did not received a copy of the GNU General Public License
% with this software, see <http://www.gnu.org/licenses/>.
%
realsmall=sqrt(eps)*10;
%realsmall=1e-3;
[n jnk] = size(A);
root = [diag(A) diag(B)];
% first sort on the non-zero root criterion
xdown0 = abs(root(:,1))<realsmall;
xdown = (xdown0 | (real(root(:,2)./(xdown0+root(:,1))) > stake));
for i = n:-1:1
m=0;
for j=i:-1:1
if xdown0(j)
m=j;
break
end
end
if (m==0)
break
end
for k=m:1:i-1
[A B Q Z] = qzswitch(k,A,B,Q,Z);
root=[diag(A) diag(B)];
xdown0(k:k+1)=flipud(xdown0(k:k+1));
xdown(k:k+1)=flipud(xdown(k:k+1));
if any(xdown(k:k+1)~=(xdown0(k:k+1) | (real(root(k:k+1,2)./(xdown0(k:k+1)+root(k:k+1,1)))) > stake))
disp('xdown shift during 0 pack at i,k:')
disp([i k])
end
end
end
% now repeat, using the stake criterion
for i = n:-1:1
m=0;
for j=i:-1:1
if xdown(j)
m=j;
break
end
end
if (m==0)
return
end
for k=m:1:i-1
gevOld=root(k:k+1,:);
[A B Q Z] = qzswitch(k,A,B,Q,Z);
root=[diag(A) diag(B)];
xdown0(k:k+1)=flipud(xdown0(k:k+1));
xdown(k:k+1)=flipud(xdown(k:k+1));
if any(xdown(k:k+1)~=(xdown0(k:k+1) | (real(root(k:k+1,2)./(xdown0(k:k+1)+root(k:k+1,1)))) > stake))
disp('xdown shift during pos pack at i,k:')
disp([i k])
gev=root(k:k+1,:);
[gevOld gevOld(:,1).\gevOld(:,2);gev gev(:,1).\gev(:,2)]
end
end
end
|
