function [A,B,Q,Z] = qzswitch(i,A,B,Q,Z) Takes U.T. matrices A, B, orthonormal matrices Q,Z, interchanges diagonal elements i and i+1 of both A and B, while maintaining Q'AZ' and Q'BZ' unchanged. If diagonal elements of A and B are zero at matching positions, the returned A will have zeros at both positions on the diagonal. This is natural behavior if this routine is used to drive all zeros on the diagonal of A to the lower right, but in this case the qz transformation is not unique and it is not possible simply to switch the positions of the diagonal elements of both A and B.
0001 function [A,B,Q,Z] = qzswitch(i,A,B,Q,Z) 0002 %function [A,B,Q,Z] = qzswitch(i,A,B,Q,Z) 0003 % 0004 % Takes U.T. matrices A, B, orthonormal matrices Q,Z, interchanges 0005 % diagonal elements i and i+1 of both A and B, while maintaining 0006 % Q'AZ' and Q'BZ' unchanged. If diagonal elements of A and B 0007 % are zero at matching positions, the returned A will have zeros at both 0008 % positions on the diagonal. This is natural behavior if this routine is used 0009 % to drive all zeros on the diagonal of A to the lower right, but in this case 0010 % the qz transformation is not unique and it is not possible simply to switch 0011 % the positions of the diagonal elements of both A and B. 0012 realsmall=sqrt(eps)*10; 0013 %realsmall=1e-3; 0014 a = A(i,i); d = B(i,i); b = A(i,i+1); e = B(i,i+1); 0015 c = A(i+1,i+1); f = B(i+1,i+1); 0016 % A(i:i+1,i:i+1)=[a b; 0 c]; 0017 % B(i:i+1,i:i+1)=[d e; 0 f]; 0018 if (abs(c)<realsmall & abs(f)<realsmall) 0019 if abs(a)<realsmall 0020 % l.r. coincident 0's with u.l. of A=0; do nothing 0021 return 0022 else 0023 % l.r. coincident zeros; put 0 in u.l. of a 0024 wz=[b; -a]; 0025 wz=wz/sqrt(wz'*wz); 0026 wz=[wz [wz(2)';-wz(1)'] ]; 0027 xy=eye(2); 0028 end 0029 elseif (abs(a)<realsmall & abs(d)<realsmall) 0030 if abs(c)<realsmall 0031 % u.l. coincident zeros with l.r. of A=0; do nothing 0032 return 0033 else 0034 % u.l. coincident zeros; put 0 in l.r. of A 0035 wz=eye(2); 0036 xy=[c -b]; 0037 xy=xy/sqrt(xy*xy'); 0038 xy=[[xy(2)' -xy(1)'];xy]; 0039 end 0040 else 0041 % usual case 0042 wz = [c*e-f*b, (c*d-f*a)']; 0043 xy = [(b*d-e*a)', (c*d-f*a)']; 0044 n = sqrt(wz*wz'); 0045 m = sqrt(xy*xy'); 0046 if m<eps*100 0047 % all elements of A and B proportional 0048 return 0049 end 0050 wz = n\wz; 0051 xy = m\xy; 0052 wz = [wz; -wz(2)', wz(1)']; 0053 xy = [xy;-xy(2)', xy(1)']; 0054 end 0055 A(i:i+1,:) = xy*A(i:i+1,:); 0056 B(i:i+1,:) = xy*B(i:i+1,:); 0057 A(:,i:i+1) = A(:,i:i+1)*wz; 0058 B(:,i:i+1) = B(:,i:i+1)*wz; 0059 Z(:,i:i+1) = Z(:,i:i+1)*wz; 0060 Q(i:i+1,:) = xy*Q(i:i+1,:);