      SUBROUTINE HTRIBK (NM, N, AR, AI, TAU, M, ZR, ZI)
      IMPLICIT NONE
C
C A TRANSLATION OF A COMPLEX ANALOGUE OF THE ALGOL PROCEDURE TRBAK1,
C   NUM. MATH. 11, 181-195 (1968) BY MARTIN, REINSCH, AND WILKINSON.
C   HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226 (1971).
C
C FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN MATRIX BY BACK TRANSFORMING
C   THOSE OF THE CORRESPONDING REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED
C   BY HTRIDI.
C
C ON INPUT:
C   NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
C      AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
C
C   N IS THE ORDER OF THE MATRIX.
C
C   AR & AI CONTAIN INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN
C      THE REDUCTION BY HTRIDI IN THEIR FULL LOWER TRIANGLES EXCEPT FOR THE
C      DIAGONAL OF AR.
C
C   TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
C
C   M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
C
C   ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS.
C
C ON OUTPUT:
C   ZR & ZI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE
C      TRANSFORMED EIGENVECTORS IN THEIR FIRST M COLUMNS.
C
C   NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR
C      IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C   APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
      INTEGER NM, N, M
      DOUBLE PRECISION AR(NM,N), AI(NM,N), TAU(2,N), ZR(NM,M), ZI(NM,M)
C
      INTEGER I, J, K, L
      DOUBLE PRECISION H, S, SI
C
      DOUBLE PRECISION FLOP
C
C
      IF (M.EQ.0) GO TO 200
C
C ***      TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC TRIDIAGONAL
C           MATRIX TO THOSE OF THE HERMITIAN TRIDIAGONAL MATRIX
      DO 51 K = 1, N
        DO 50 J = 1, M
          ZI(K,J) = FLOP (-ZR(K,J)*TAU(2,K))
          ZR(K,J) = FLOP (ZR(K,J)*TAU(1,K))
50      CONTINUE
51    CONTINUE
      IF (N.EQ.1) GO TO 200
C
C ***      RECOVER AND APPLY THE HOUSEHOLDER MATRICES
      DO 140 I = 2, N
        L = I-1
        H = AI(I,I)
        IF (H.EQ.0.0D0) GO TO 140
        DO 130 J = 1, M
          S = 0.0D0
          SI = 0.0D0
          DO 110 K = 1, L
            S = FLOP (S+AR(I,K)*ZR(K,J)-AI(I,K)*ZI(K,J))
            SI = FLOP (SI+AR(I,K)*ZI(K,J)+AI(I,K)*ZR(K,J))
110       CONTINUE
C
C ***      DOUBLE DIVISIONS AVOID POSSIBLE UNDERFLOW
          S = FLOP ((S/H)/H)
          SI = FLOP ((SI/H)/H)
          DO 120 K = 1, L
            ZR(K,J) = FLOP (ZR(K,J)-S*AR(I,K)-SI*AI(I,K))
            ZI(K,J) = FLOP (ZI(K,J)-SI*AR(I,K)+S*AI(I,K))
120       CONTINUE
130     CONTINUE
140   CONTINUE
C
200   CONTINUE
      RETURN
      END
