% FIGURE 3.4.5 
      function A = LR(A,N,ERRLIM)
      if (N <= 2) 
         return
      end
%                              USE LR ITERATION TO REDUCE HESSENBERG
%                              MATRIX A TO QUASI-TRIANGULAR FORM
%
%                              PIVOT = 'T' IF PIVOTING ALLOWED
      PIVOT = 'T'; 
      NITER = 1000*N;
      for ITER=1:NITER
%                              REDUCE A TO UPPER TRIANGULAR FORM USING
%                              GAUSSIAN ELIMINATION (PREMULTIPLY BY
%                              Mij^(-1) MATRICES)
         for I=1:N-1
            if (abs(A(I,I)) < ERRLIM & PIVOT == 'T') 
%                              SWITCH ROWS I AND I+1 IF NECESSARY
%                              (PREMULTIPLY BY Pi,i+1)
               IPERM(I) = I+1;
               for K=I:N
                  TEMP = A(I+1,K);
                  A(I+1,K) = A(I,K);
                  A(I,K) = TEMP;
               end
            else 
               IPERM(I) = I;
            end 
            if (abs(A(I+1,I)) < ERRLIM) 
               R = 0.0;
            else 
               if (abs(A(I,I)) < ERRLIM) 
                  error('***** LR ITERATION DOES NOT CONVERGE *****')
               end
               R = A(I+1,I)/A(I,I);
            end
%                              USE SAVE TO SAVE R FOR POST-
%                              MULTIPLICATION PHASE
            SAVE(I) = R;
            if (R == 0.0) 
               continue
            end
%                              IF MATRIX TRIDIAGONAL, AND PIVOTING NOT
%                              DONE, LIMITS ON K CAN BE:  K = I , I+1
            for K=I:N
               A(I+1,K) = A(I+1,K) - R*A(I,K);
            end
         end
%                              NOW POSTMULTIPLY BY Mij MATRICES
         for I=1:N-1
            if (IPERM(I) ~= I) 
%                              SWITCH COLUMNS I AND I+1 IF NECESSARY
%                              (POSTMULTIPLY BY Pi,i+1^(-1) = Pi,i+1)
               for K=1:I+1
                  TEMP = A(K,I+1);
                  A(K,I+1) = A(K,I);
                  A(K,I) = TEMP;
               end
            end
            R = SAVE(I);
            if (R == 0.0)
               continue
            end
%                              IF MATRIX TRIDIAGONAL, AND PIVOTING NOT
%                              DONE, LIMITS ON K CAN BE:  K = I , I+1
            for K=1:I+1
               A(K,I) = A(K,I) + R*A(K,I+1);
            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
      error('***** LR ITERATION DOES NOT CONVERGE *****')
