% FIGURE 3.3.4
      function A = QR(A,N,ERRLIM)
      if (N <= 2)
         return
      end
%                              USE QR ITERATION TO REDUCE HESSENBERG
%                              MATRIX A TO QUASI-TRIANGULAR FORM
      NITER = 1000*N;
      for ITER=1:NITER
%                              REDUCE A TO UPPER TRIANGULAR FORM USING
%                              ORTHOGONAL REDUCTION (PREMULTIPLY BY
%                              Qij^T MATRICES)
         for I=1:N-1
            if (A(I+1,I) == 0.0)
               C = 1.0;
               S = 0.0;
            else
               DEN = sqrt(A(I,I)^2 + A(I+1,I)^2);
               C = A(I,I)/DEN;
               S = A(I+1,I)/DEN;
            end
%                              USE SAVE TO SAVE C,S FOR POST-
%                              MULTIPLICATION PHASE
            SAVE(1,I) = C;
            SAVE(2,I) = S;
            if (S == 0.0)
               continue
            end
%                              IF MATRIX SYMMETRIC, LIMITS ON K
%                              CAN BE:  K = I : min(I+2,N)
            for K=I:N
               PIK = C*A(I,K) + S*A(I+1,K);
               PJK =-S*A(I,K) + C*A(I+1,K);
               A(I,K)   = PIK;
               A(I+1,K) = PJK;
            end
         end
%                              NOW POSTMULTIPLY BY Qij MATRICES
         for I=1:N-1
            C = SAVE(1,I);
            S = SAVE(2,I);
            if (S == 0.0)
               continue
            end
%                              IF MATRIX SYMMETRIC, LIMITS ON K
%                              CAN BE:  K = max(1,I-1) : I+1
            for K=1:I+1
               BKI = C*A(K,I) + S*A(K,I+1);
               BKJ =-S*A(K,I) + C*A(K,I+1);
               A(K,I)   = BKI;
               A(K,I+1) = BKJ;
            end
         end
%                              SET NEARLY ZERO SUBDIAGONALS TO ZERO,
%                              TO AVOID UNDERFLOW.
         for I=1:N-1
            if (abs(A(I+1,I)) < ERRLIM)
               A(I+1,I) = 0.0;
            end
         end
%                              CHECK FOR CONVERGENCE TO "QUASI-
%                              TRIANGULAR" FORM.
         ICONV = 1;
         for I=2:N-1
            if (A(I,I-1) ~= 0.0 & A(I+1,I) ~= 0.0)
               ICONV = 0;
            end
         end
         if (ICONV == 1)
            return
         end
      end
%                              HAS NOT CONVERGED IN NITER ITERATIONS
      error('***** QR ITERATION DOES NOT CONVERGE *****')