      SUBROUTINE TQLI1(D,E,N,NP,Z)
C
C     * THIS SUBROUTINE IS TAKEN FROM NUMERICAL RECIPES.
C     * IT SOLVES EIGENVALUE PROBLEM FOR A 3-DIAGONAL MATRIX.
C     * THE ORIGINAL NAME OF THE SUBROUTINE HAS BEEN CHANGED FROM "TQLI" TO "TQLI1".
C     * "PAUSE" STATEMENT IS REPLACED WITH WITH A CALL TO THE "XIT" SUBROUTINE.
C     * THE MAXIMAL NUMBER OF ITERATIONS IS INCREASED FROM 30 TO 1000.
C
C     * QL ALGORITHM WITH IMPLICIT SHIFTS, TO DETERMINE THE EIGENVALUES AND
C     * EIGENVECTORS OF A REAL, SYMMETRIC, TRIDIAGONAL MATRIX, OR OF A REAL, 
C     * SYMMETRIC MATRIX PREVIOUSLY REDUCED BY TRED2. 
C     * D IS A VECTOR OF LENGTH NP. ON INPUT, ITS FIRST N ELEMENTS ARE
C     * THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. 
C     * ON OUTPUT, IT RETURNS THE EIGENVALUES. THE VECTOR E INPUTS THE SUBDIAGONAL
C     * ELEMENTS OF THE TRIDIAGONAL MATRIX, WITH E(1) ARBITRARY. 
C     * ON OUTPUT E IS DESTROYED.
C     * IF THE EIGENVECTORS ARE DESIRED, THE MATRIX Z (N BY N MATRIX STORED IN 
C     * NP BY NP ARRAY) IS INPUT AS THE IDENTITY MATRIX. IF THE EIGENVECTORS OF A
C     * MATRIX THT HAS BEEN REDUCED BY TRED2 ARE REQUIRED, THEN Z IS INPUT AS THE 
C     * MATRIX OUTPUT BY TRED2. IN EITHER CASE, THE KTH COLUMN IF Z RETURNS
C     * THE NORMALIZED EIGENVECTOR CORRESPONDING TO D(K).
C
C     JUN 28/96 - SLAVA KHARIN
C
C  (C) COPR. 1986-92 NUMERICAL RECIPES SOFTWARE #=!E=#,)]UBCJ.
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      PARAMETER (MAXITER=1000)
      INTEGER N,NP
      REAL D(NP),E(NP),Z(NP,NP)
CU    USES PYTHAG
      INTEGER I,ITER,K,L,M
      REAL B,C,DD,F,G,P,R,S,PYTHAG
      DO 11 I=2,N
        E(I-1)=E(I)
11    CONTINUE
      E(N)=0.E0
      DO 15 L=1,N
        ITER=0
1       DO 12 M=L,N-1
          DD=ABS(D(M))+ABS(D(M+1))
          IF (ABS(E(M))+DD.EQ.DD) GOTO 2
12      CONTINUE
        M=N
2       IF(M.NE.L)THEN
          IF(ITER.EQ.MAXITER)THEN
             WRITE(6,'(A)') ' TOO MANY ITERATIONS IN TQLI1'
             CALL                                  XIT ('TQLI1',-1)
          ENDIF
          ITER=ITER+1
          G=(D(L+1)-D(L))/(2.E0*E(L))
          R=PYTHAG(G,1.E0)
          G=D(M)-D(L)+E(L)/(G+SIGN(R,G))
          S=1.E0
          C=1.E0
          P=0.E0
          DO 14 I=M-1,L,-1
            F=S*E(I)
            B=C*E(I)
            R=PYTHAG(F,G)
            E(I+1)=R
            IF(R.EQ.0.E0)THEN
              D(I+1)=D(I+1)-P
              E(M)=0.E0
              GOTO 1
            ENDIF
            S=F/R
            C=G/R
            G=D(I+1)-P
            R=(D(I)-G)*S+2.E0*C*B
            P=S*R
            D(I+1)=G+P
            G=C*R-B
C     OMIT LINES FROM HERE ...
            DO 13 K=1,N
              F=Z(K,I+1)
              Z(K,I+1)=S*Z(K,I)+C*F
              Z(K,I)=C*Z(K,I)-S*F
13          CONTINUE
C     ... TO HERE WHEN FINDING ONLY EIGENVALUES.
14        CONTINUE
          D(L)=D(L)-P
          E(L)=G
          E(M)=0.E0
          GOTO 1
        ENDIF
15    CONTINUE
      RETURN
      END
