      SUBROUTINE TRED2(A,N,NP,D,E)
C
C     * THIS SUBROUTINE IS TAKEN DIRECTLY FROM NUMERICAL RECIPES WITHOUT ANY CHANGES.
C     * IT TRANSFORMS A SYMMETRIC MATRIX TO A 3-DIAGONAL ONE.
C
C     * HOUSEHOLDER REDUCTION OF A REAL, SYMMETRIC, N BY N MATRIX A, STORED IN AN 
C     * NP BY NP PHYSICAL ARRAY. ON OUTPUT, A IS REPLACED BY THE ORTHOGONAL MATRIX Q
C     * EFFECTING THE TRANSFORMATION. D RETURNS THE DIAGONAL ELEMENTS OF THE 
C     * TRIDIAGONAL MATRIX, AND E THE OFF-ELEMENTS, WITH E(1)=0.
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)
      INTEGER N,NP
      REAL A(NP,NP),D(NP),E(NP)
      INTEGER I,J,K,L
      REAL F,G,H,HH,SCALE
      DO 18 I=N,2,-1
        L=I-1
        H=0.E0
        SCALE=0.E0
        IF(L.GT.1)THEN
          DO 11 K=1,L
            SCALE=SCALE+ABS(A(I,K))
11        CONTINUE
          IF(SCALE.EQ.0.E0)THEN
            E(I)=A(I,L)
          ELSE
            DO 12 K=1,L
              A(I,K)=A(I,K)/SCALE
              H=H+A(I,K)**2
12          CONTINUE
            F=A(I,L)
            G=-SIGN(SQRT(H),F)
            E(I)=SCALE*G
            H=H-F*G
            A(I,L)=F-G
            F=0.E0
            DO 15 J=1,L
C     OMIT FOLLOWING LINE IF FINDING ONLY EIGENVALUES
              A(J,I)=A(I,J)/H
              G=0.E0
              DO 13 K=1,J
                G=G+A(J,K)*A(I,K)
13            CONTINUE
              DO 14 K=J+1,L
                G=G+A(K,J)*A(I,K)
14            CONTINUE
              E(J)=G/H
              F=F+E(J)*A(I,J)
15          CONTINUE
            HH=F/(H+H)
            DO 17 J=1,L
              F=A(I,J)
              G=E(J)-HH*F
              E(J)=G
              DO 16 K=1,J
                A(J,K)=A(J,K)-F*E(K)-G*A(I,K)
16            CONTINUE
17          CONTINUE
          ENDIF
        ELSE
          E(I)=A(I,L)
        ENDIF
        D(I)=H
18    CONTINUE
C     OMIT FOLLOWING LINE IF FINDING ONLY EIGENVALUES.
      D(1)=0.E0
      E(1)=0.E0
      DO 24 I=1,N
C     DELETE LINES FROM HERE ...
        L=I-1
        IF(D(I).NE.0.E0)THEN
          DO 22 J=1,L
            G=0.E0
            DO 19 K=1,L
              G=G+A(I,K)*A(K,J)
19          CONTINUE
            DO 21 K=1,L
              A(K,J)=A(K,J)-G*A(K,I)
21          CONTINUE
22        CONTINUE
        ENDIF
C     ... TO HERE WHEN FINDING ONLY EIGENVALUES.
        D(I)=A(I,I)
C     ALSO DELETE LINES FROM HERE ...
        A(I,I)=1.E0
        DO 23 J=1,L
          A(I,J)=0.E0
          A(J,I)=0.E0
23      CONTINUE
C     ... TO HERE WHEN FINDING ONLY EIGENVALUES.
24    CONTINUE
      RETURN
      END
