qzswitch

 from Chris Sims web site
 http://eco-072399b.princeton.edu/yftp/gensys/mfiles/QZSWITCH.M

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