      SUBROUTINE SINVIT(N,D,E,W,Z,IERR,RV1,RV2,RV3,RV4,RV6)

C     Michael B. Porter 7/1/85

C     This subroutine is a modified version of TINVIT in EISPACK
     
C     Employs FORTRAN 77 features in order to improve readability.
C     Also has been specialized to the case when
C     we only want a single vector.

C     WARNING: The modifications have not been done extremely
C     carefully.  Very probably this routine will fail for
C     certain difficult cases for which the original tinvit would
C     not.  For instance, problems with nearly degenerate eigenvectors
C     or separable problems.

C     Refer to the ALGOL procedure TRISTURM by Peters and Wilkinson.
C     Handbook for Auto. Comp., VOL.II-Linear Algebra, 418-439(1971).

C     This subroutine finds the eigenvector of a tridiagonal
C     symmetric matrix corresponding to a specified eigenvalue,
C     using inverse iteration.

C     On input

C        N is the order of the matrix.
C        D contains the diagonal elements of the input matrix.
C        E contains the subdiagonal elements of the input matrix
C          in its last N-1 positions.  E(1) is arbitrary.
C        W contains the eigenvalue

C     On output

C        All input arrays are unaltered.
C        Z contains the eigenvector.
C          Any vector which fails to converge is set to zero.
C        IERR is set to
C           0       For normal return,
C          -1       If the eigenvector fails to converge in 5 iterations.

C        RV1, RV2, RV3, RV4, and RV6 are temporary storage arrays.

      REAL*8 MACHEP
      PARAMETER ( MACHEP = 1.0D-10 )

      INTEGER I, IERR, N
      REAL*8 EPS3,EPS4,NORM,UK
      COMPLEX*16 U,V, XU,D(N),E(N),
     +       RV1(N),RV2(N),RV3(N),RV4(N),RV6(N),W,Z(N)
      COMPLEX*16 ZTEMP

C     --- use DIMAG in f77, AIMAG in f90 below
      RLEN( ZTEMP ) = ABS( DBLE( ZTEMP ) ) + ABS( DIMAG( ZTEMP ) )

C     --- Compute (infinity) norm of matrix

      NORM = RLEN( D(1) )

      DO I = 2, N
         NORM = NORM + RLEN( D(I) ) + RLEN( E(I) )
      END DO

C     EPS3 Replaces zero pivots
C     EPS4 a very small number used in scaling down the iterates

      EPS3 = MACHEP * NORM
      UK   = N
      EPS4 = UK * EPS3
      UK   = EPS4 / DSQRT( UK )

C     *** Elimination with interchanges ***

      XU = 1.0
      U  = D( 1 ) - W
      V  = E( 2 )
      DO I = 2, N

         IF ( RLEN( E(I) ) .GE. RLEN( U ) ) THEN
            XU = U / E(I)
            RV4( I   ) = XU
            RV1( I-1 ) = E( I )
            RV2( I-1 ) = D( I ) - W
            RV3( I-1 ) = 0.0
            IF ( I .NE. N ) RV3( I-1 ) = E( I+1 )
            U = V - XU * RV2( I-1 )
            V =   - XU * RV3( I-1 )
         ELSE
            XU = E(I) / U
            RV4( I   ) = XU
            RV1( I-1 ) = U
            RV2( I-1 ) = V
            RV3( I-1 ) = 0.0
            U = D( I ) - W - XU * V
            IF ( I .NE. N ) V = E( I+1 )
         ENDIF

      END DO

      IF ( U .EQ. 0.0 ) U = EPS3
      RV1( N ) = U
      RV2( N ) = 0.0
      RV3( N ) = 0.0

C     *** Main loop of inverse iteration ***

C     --- Set initial vector to all UK

      DO I = 1, N
         RV6( I ) = UK
      END DO

      DO ITER = 1, 5

C        ---  Back substitution
         DO I = N, 1, -1
            RV6( I ) = ( RV6(I) - U * RV2(I) - V * RV3(I) ) / RV1(I)
            V = U
            U = RV6( I )
         END DO

C        --- Compute norm of vector
         NORM = 0.0
         DO I = 1, N
            NORM = NORM + RLEN( RV6( I ) )
         END DO

C        --- Convergence?
         IF ( NORM .LT. 1.0 ) THEN

C           --- Scale the vector down
            XU = EPS4 / NORM
            DO I = 1, N
               RV6( I ) = RV6( I ) * XU
            END DO

C           --- Forward elimination
            DO I = 2, N
               U = RV6( I )

C              IF RV1(I-1) = E(I), a row interchange
C              was performed earlier in the triangularization process

               IF ( RV1( I-1 ) .EQ. E( I ) ) THEN
                  U          = RV6( I-1 )
                  RV6( I-1 ) = RV6( I   )
               ENDIF
               RV6( I ) = U - RV4( I ) * RV6( I-1 )
            END DO
         ELSE

C           --- Convergence: Set norm to unity and return

            IERR = 0
            XU = 1.0 / NORM

            DO I = 1, N
               Z( I ) = XU * RV6( I )
            END DO

            RETURN

         ENDIF
      END DO   ! next iteration

C     Fall through the loop means failure to converge

      IERR = -1

      RETURN
      END
