% FIGURE 3.6.1      
      function [A,B] = HESSQZ(A,B,N)
%
%  FUNCTION HESSQZ REDUCES THE GENERALIZED EIGENVALUE PROBLEM
%
%                A*X = LAMBDA*B*X
%
%  TO A SIMILAR PROBLEM WITH THE SAME EIGENVALUES, WHERE A IS 
%  UPPER HESSENBERG AND B IS UPPER TRIANGULAR. 
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    A      - THE N BY N A MATRIX.              A IS UPPER HESSENBERG. 
%
%    B      - THE N BY N B MATRIX.              B IS UPPER TRIANGULAR.  
%
%    N      - THE SIZE OF MATRICES A AND B.
%
%-----------------------------------------------------------------------
%                              PREMULTIPLY A AND B BY ORTHOGONAL MATRIX
%                              (PRODUCT OF GIVENS MATRICES) Q, SUCH   
%                              THAT QB IS UPPER TRIANGULAR.
      for I=1:N-1
         for J=I+1:N
            if (B(J,I) == 0.0)
               continue 
            end
            DEN = sqrt(B(I,I)^2+B(J,I)^2);
            S = -B(J,I)/DEN;
            C =  B(I,I)/DEN;
            for K=I:N
               BIK = C*B(I,K)-S*B(J,K); 
               BJK = S*B(I,K)+C*B(J,K);
               B(I,K) = BIK;
               B(J,K) = BJK;
            end
            for K=1:N
               AIK = C*A(I,K)-S*A(J,K); 
               AJK = S*A(I,K)+C*A(J,K);
               A(I,K) = AIK;
               A(J,K) = AJK;
            end
         end
      end
%                              PREMULTIPLY A AND B BY ORTHOGONAL MATRIX
%                              Q, AND POSTMULTIPLY BY ORTHOGONAL MATRIX
%                              Z, SUCH THAT QAZ IS UPPER HESSENBERG AND  
%                              QBZ IS STILL UPPER TRIANGULAR      
      for I=1:N-2
         for J=N:-1:I+2
            if (A(J,I) == 0.0) 
               continue
            end
%                              PREMULTIPLY A TO ZERO A(J,I)   
            DEN = sqrt(A(J-1,I)^2+A(J,I)^2);
            S = -A(J,I)/DEN;
            C =  A(J-1,I)/DEN;
            for K=I:N
               A1K = C*A(J-1,K) - S*A(J,K);
               A2K = S*A(J-1,K) + C*A(J,K);
               A(J-1,K) = A1K;
               A(J,K) = A2K;
            end
%                              PREMULTIPLY B BY SAME MATRIX, CREATING 
%                              NEW NONZERO B(J,J-1)
            for K=J-1:N
               B1K = C*B(J-1,K) - S*B(J,K);
               B2K = S*B(J-1,K) + C*B(J,K);
               B(J-1,K) = B1K;
               B(J,K) = B2K;
            end
            if (B(J,J-1) == 0.0)
               continue
            end
%                              POSTMULTIPLY B TO ZERO B(J,J-1)     
            DEN = sqrt(B(J,J-1)^2+B(J,J)^2);
            S = -B(J,J-1)/DEN;
            C =  B(J,J)/DEN;
            for K=1:J
               BK1 = C*B(K,J-1) + S*B(K,J);
               BK2 = -S*B(K,J-1) + C*B(K,J);  
               B(K,J-1) = BK1;
               B(K,J) = BK2;
            end
%                              POSTMULTIPLY A BY SAME MATRIX
            for K=1:N
               AK1 = C*A(K,J-1) + S*A(K,J);
               AK2 = -S*A(K,J-1) + C*A(K,J);  
               A(K,J-1) = AK1;
               A(K,J) = AK2;
            end
         end
      end
