% 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
      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 (A(I+1,I) == 0.0)
               R = 0.0;
            else
%                              IF PIVOTING NECESSARY, GIVE UP
               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, 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
            R = SAVE(I);
            if (R == 0.0)
               continue
            end
%                              IF MATRIX TRIDIAGONAL, 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 *****')