C/ FIGURE 3.4.5
      SUBROUTINE LR(A,N,ERRLIM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                              DECLARATIONS FOR ARGUMENTS
      DOUBLE PRECISION A(N,N),ERRLIM
      INTEGER N
C                              DECLARATIONS FOR LOCAL VARIABLES
      DOUBLE PRECISION SAVE(N)
      INTEGER IPERM(N)
      LOGICAL PIVOT
      IF (N.LE.2) RETURN
C                              USE LR ITERATION TO REDUCE HESSENBERG
C                              MATRIX A TO QUASI-TRIANGULAR FORM
C
C                              PIVOT = .TRUE. IF PIVOTING ALLOWED
      PIVOT = .TRUE.
      NITER = 1000*N
      DO 45 ITER=1,NITER
C                              REDUCE A TO UPPER TRIANGULAR FORM USING
C                              GAUSSIAN ELIMINATION (PREMULTIPLY BY
C                              Mij**(-1) MATRICES)
         DO 15 I=1,N-1
            IF (ABS(A(I,I)).LT.ERRLIM .AND. PIVOT) THEN
C                              SWITCH ROWS I AND I+1 IF NECESSARY
C                              (PREMULTIPLY BY Pi,i+1)
               IPERM(I) = I+1
               DO 5 K=I,N
                  TEMP = A(I+1,K)
                  A(I+1,K) = A(I,K)
                  A(I,K) = TEMP
    5          CONTINUE
            ELSE
               IPERM(I) = I
            ENDIF
            IF (ABS(A(I+1,I)).LT.ERRLIM) THEN
               R = 0
            ELSE
               IF (ABS(A(I,I)).LT.ERRLIM) GO TO 50
               R = A(I+1,I)/A(I,I)
            ENDIF
C                              USE SAVE TO SAVE R FOR POST-
C                              MULTIPLICATION PHASE
            SAVE(I) = R
            IF (R.EQ.0.0) GO TO 15
C                              IF MATRIX TRIDIAGONAL, AND PIVOTING NOT
C                              DONE, LIMITS ON K CAN BE:  K = I , I+1
            DO 10 K=I,N
               A(I+1,K) = A(I+1,K) - R*A(I,K)
   10       CONTINUE
   15    CONTINUE
C                              NOW POSTMULTIPLY BY Mij MATRICES
         DO 30 I=1,N-1
            IF (IPERM(I).NE.I) THEN
C                              SWITCH COLUMNS I AND I+1 IF NECESSARY
C                              (POSTMULTIPLY BY Pi,i+1**(-1) = Pi,i+1)
               DO 20 K=1,I+1
                  TEMP = A(K,I+1)
                  A(K,I+1) = A(K,I)
                  A(K,I) = TEMP
   20          CONTINUE
            ENDIF
            R = SAVE(I)
            IF (R.EQ.0.0) GO TO 30
C                              IF MATRIX TRIDIAGONAL, AND PIVOTING NOT
C                              DONE, LIMITS ON K CAN BE:  K = I , I+1
            DO 25 K=1,I+1
               A(K,I) = A(K,I) + R*A(K,I+1)
   25       CONTINUE
   30    CONTINUE
C                              SET NEARLY ZERO SUBDIAGONALS TO ZERO,
C                              TO AVOID UNDERFLOW.
         DO 35 I=1,N-1
            IF (ABS(A(I+1,I)).LT.ERRLIM) A(I+1,I) = 0.0
   35    CONTINUE
C                              CHECK FOR CONVERGENCE TO "QUASI-
C                              TRIANGULAR" FORM.
         ICONV = 1
         DO 40 I=2,N-1
            IF (A(I,I-1).NE.0.0 .AND. A(I+1,I).NE.0.0) ICONV = 0
   40    CONTINUE
         IF (ICONV.EQ.1) RETURN
   45 CONTINUE
C                              HAS NOT CONVERGED IN NITER ITERATIONS
   50 PRINT 55
   55 FORMAT (' ***** LR ITERATION DOES NOT CONVERGE *****')
      RETURN
      END