1*DECK HTRIDI 2 SUBROUTINE HTRIDI (NM, N, AR, AI, D, E, E2, TAU) 3C***BEGIN PROLOGUE HTRIDI 4C***PURPOSE Reduce a complex Hermitian matrix to a real symmetric 5C tridiagonal matrix using unitary similarity 6C transformations. 7C***LIBRARY SLATEC (EISPACK) 8C***CATEGORY D4C1B1 9C***TYPE SINGLE PRECISION (HTRIDI-S) 10C***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK 11C***AUTHOR Smith, B. T., et al. 12C***DESCRIPTION 13C 14C This subroutine is a translation of a complex analogue of 15C the ALGOL procedure TRED1, NUM. MATH. 11, 181-195(1968) 16C by Martin, Reinsch, and Wilkinson. 17C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). 18C 19C This subroutine reduces a COMPLEX HERMITIAN matrix 20C to a real symmetric tridiagonal matrix using 21C unitary similarity transformations. 22C 23C On INPUT 24C 25C NM must be set to the row dimension of the two-dimensional 26C array parameters, AR and AI, as declared in the calling 27C program dimension statement. NM is an INTEGER variable. 28C 29C N is the order of the matrix A=(AR,AI). N is an INTEGER 30C variable. N must be less than or equal to NM. 31C 32C AR and AI contain the real and imaginary parts, respectively, 33C of the complex Hermitian input matrix. Only the lower 34C triangle of the matrix need be supplied. AR and AI are two- 35C dimensional REAL arrays, dimensioned AR(NM,N) and AI(NM,N). 36C 37C On OUTPUT 38C 39C AR and AI contain some information about the unitary trans- 40C formations used in the reduction in the strict lower triangle 41C of AR and the full lower triangle of AI. The rest of the 42C matrices are unaltered. 43C 44C D contains the diagonal elements of the real symmetric 45C tridiagonal matrix. D is a one-dimensional REAL array, 46C dimensioned D(N). 47C 48C E contains the subdiagonal elements of the real tridiagonal 49C matrix in its last N-1 positions. E(1) is set to zero. 50C E is a one-dimensional REAL array, dimensioned E(N). 51C 52C E2 contains the squares of the corresponding elements of E. 53C E2(1) is set to zero. E2 may coincide with E if the squares 54C are not needed. E2 is a one-dimensional REAL array, 55C dimensioned E2(N). 56C 57C TAU contains further information about the transformations. 58C TAU is a one-dimensional REAL array, dimensioned TAU(2,N). 59C 60C Calls PYTHAG(A,B) for sqrt(A**2 + B**2). 61C 62C Questions and comments should be directed to B. S. Garbow, 63C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 64C ------------------------------------------------------------------ 65C 66C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, 67C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- 68C system Routines - EISPACK Guide, Springer-Verlag, 69C 1976. 70C***ROUTINES CALLED PYTHAG 71C***REVISION HISTORY (YYMMDD) 72C 760101 DATE WRITTEN 73C 890831 Modified array declarations. (WRB) 74C 890831 REVISION DATE from Version 3.2 75C 891214 Prologue converted to Version 4.0 format. (BAB) 76C 920501 Reformatted the REFERENCES section. (WRB) 77C***END PROLOGUE HTRIDI 78C 79 INTEGER I,J,K,L,N,II,NM,JP1 80 REAL AR(NM,*),AI(NM,*),D(*),E(*),E2(*),TAU(2,*) 81 REAL F,G,H,FI,GI,HH,SI,SCALE 82 REAL PYTHAG 83C 84C***FIRST EXECUTABLE STATEMENT HTRIDI 85 TAU(1,N) = 1.0E0 86 TAU(2,N) = 0.0E0 87C 88 DO 100 I = 1, N 89 100 D(I) = AR(I,I) 90C .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... 91 DO 300 II = 1, N 92 I = N + 1 - II 93 L = I - 1 94 H = 0.0E0 95 SCALE = 0.0E0 96 IF (L .LT. 1) GO TO 130 97C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... 98 DO 120 K = 1, L 99 120 SCALE = SCALE + ABS(AR(I,K)) + ABS(AI(I,K)) 100C 101 IF (SCALE .NE. 0.0E0) GO TO 140 102 TAU(1,L) = 1.0E0 103 TAU(2,L) = 0.0E0 104 130 E(I) = 0.0E0 105 E2(I) = 0.0E0 106 GO TO 290 107C 108 140 DO 150 K = 1, L 109 AR(I,K) = AR(I,K) / SCALE 110 AI(I,K) = AI(I,K) / SCALE 111 H = H + AR(I,K) * AR(I,K) + AI(I,K) * AI(I,K) 112 150 CONTINUE 113C 114 E2(I) = SCALE * SCALE * H 115 G = SQRT(H) 116 E(I) = SCALE * G 117 F = PYTHAG(AR(I,L),AI(I,L)) 118C .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T .......... 119 IF (F .EQ. 0.0E0) GO TO 160 120 TAU(1,L) = (AI(I,L) * TAU(2,I) - AR(I,L) * TAU(1,I)) / F 121 SI = (AR(I,L) * TAU(2,I) + AI(I,L) * TAU(1,I)) / F 122 H = H + F * G 123 G = 1.0E0 + G / F 124 AR(I,L) = G * AR(I,L) 125 AI(I,L) = G * AI(I,L) 126 IF (L .EQ. 1) GO TO 270 127 GO TO 170 128 160 TAU(1,L) = -TAU(1,I) 129 SI = TAU(2,I) 130 AR(I,L) = G 131 170 F = 0.0E0 132C 133 DO 240 J = 1, L 134 G = 0.0E0 135 GI = 0.0E0 136C .......... FORM ELEMENT OF A*U .......... 137 DO 180 K = 1, J 138 G = G + AR(J,K) * AR(I,K) + AI(J,K) * AI(I,K) 139 GI = GI - AR(J,K) * AI(I,K) + AI(J,K) * AR(I,K) 140 180 CONTINUE 141C 142 JP1 = J + 1 143 IF (L .LT. JP1) GO TO 220 144C 145 DO 200 K = JP1, L 146 G = G + AR(K,J) * AR(I,K) - AI(K,J) * AI(I,K) 147 GI = GI - AR(K,J) * AI(I,K) - AI(K,J) * AR(I,K) 148 200 CONTINUE 149C .......... FORM ELEMENT OF P .......... 150 220 E(J) = G / H 151 TAU(2,J) = GI / H 152 F = F + E(J) * AR(I,J) - TAU(2,J) * AI(I,J) 153 240 CONTINUE 154C 155 HH = F / (H + H) 156C .......... FORM REDUCED A .......... 157 DO 260 J = 1, L 158 F = AR(I,J) 159 G = E(J) - HH * F 160 E(J) = G 161 FI = -AI(I,J) 162 GI = TAU(2,J) - HH * FI 163 TAU(2,J) = -GI 164C 165 DO 260 K = 1, J 166 AR(J,K) = AR(J,K) - F * E(K) - G * AR(I,K) 167 1 + FI * TAU(2,K) + GI * AI(I,K) 168 AI(J,K) = AI(J,K) - F * TAU(2,K) - G * AI(I,K) 169 1 - FI * E(K) - GI * AR(I,K) 170 260 CONTINUE 171C 172 270 DO 280 K = 1, L 173 AR(I,K) = SCALE * AR(I,K) 174 AI(I,K) = SCALE * AI(I,K) 175 280 CONTINUE 176C 177 TAU(2,L) = -SI 178 290 HH = D(I) 179 D(I) = AR(I,I) 180 AR(I,I) = HH 181 AI(I,I) = SCALE * SQRT(H) 182 300 CONTINUE 183C 184 RETURN 185 END 186