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
|
function [A,B,P,R] = nqzdiv(stake,A,B,V,P)
%function [A,B,Q,Z] = nqzdiv(stake,A,B,V,P)
%
% Takes U.T. matrices A, B, V, and P and rearranges them
% so that all cases of abs(B(i,i)/A(i,i))>stake are in lower right
% corner, while the output A and B are diagonal, and P*A*R=G0, P*B*R=G1.
% Note: the input A and B are U.T.; the output A and B are diagonal.
%
% Modified by T.Zha, 5/26/97
%
% 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/>.
%
[n jnk] = size(A);
ald = diag(diag(A)); % d: diagoanl
bed = diag(diag(B));
R = inv(V);
root = abs([diag(ald) diag(bed)]);
root(:,1) = root(:,1)-(root(:,1)<1.e-13).*(root(:,1)+root(:,2));
root(:,2) = root(:,2)./root(:,1);
for i = n:-1:1
m=0;
for j=i:-1:1
if (root(j,2) > stake | root(j,2) < -.1)
m=j;
break
end
end
if (m==0)
return
elseif (m<i-1)
pindx = 1:n; % T.Z., 5/26/97
pindx(m)=i-1;
pindx(i-1)=m;
%*** permutation begins
ald = ald(pindx,pindx);
bed = bed(pindx,pindx);
P = P(:,pindx); % column permutation, T.Z., 5/26/97
R = R(pindx,:); % row permutation
end
end
%*** return the new Q, Z, A, B
A = diag(1./diag(bed));
B = diag(1./diag(ald));
|
