1! 2C ****************************************************************** 3C * HERMDP.FOR SOLVES EIGENVALUE PROBLEM * 4C * H X = E X * 5C * WHERE H IS A HERMITIAN MATRIX, X ARE THE EIGENVECTORS AND E * 6C * THE EIGENVALUES. * 7C * OTTO F. SANKEY * 8C * DEPT. OF PHYSICS * 9C * ARIZONA STATE UNIVERSITY * 10C * 602 965-4334 * 11C * * 12C * DOUBLE PRECISION 5/15/85 * 13c * updated 2/10/89 ofs. * 14C ******************************************************************* 15C 16C TO USE: 17C SUBROUTINE HERMDP(H,ZR,ZI,E,IDIM,M,W1,W2,icall) 18C H = INPUT HAMILTONIAN. 19c it contains real and imaginary parts of h as follows: 20c h(i,j) = real (h(i,j)) , [ i.ge.j, put Real(lower triangle) 21c in lower triangle]. 22c h(j,i) = imag (h(i,j)) , [ i.ge.j, put imag(lower triangle) 23c in lower triangle]. 24c 25C IDIM = DIMENSION OF H EXACTLY AS IN CALLING ROUTINE. 26C M = DIMENSION OF MATRIX YOU WANT DIAGONALIZED (M.LE.IDIM) 27C ZR,ZI = OUTPUT EIGENVECTORS (R MEANS REAL PART AND I IMAG. PART). 28c zr(i,j) is the real part of the i'th component of the j'th 29c eigenvector -- zr(component, eigenvector). 30C E = OUTPUT EIGENVALUES 31C icall= 3 for eigenvectors AND eigenvalues, icall.ne.3 for eigenvalues 32c only. 33C W1, W2 ARE WORK VECTORS, W1(IDIM),W2(2,IDIM) 34c aid is a work vector dimensioned 501. for bigger matrices, please 35c redimension. 36c 37c program computes trace & compares with sum of eigenvalues as a check. 38C 39C CALL SUBROUTINES: DPHTRI,DPIMTQ,DPHTBK 40c ============================================================== 41 SUBROUTINE hERMDP(H,ZR,ZI,E,IDIM,M,W1,W2,icall) 42c ============================================================== 43c DOUBLE PRECISION H(IDIM,M),E(M),AID(501) 44 integer :: idim, m, icall 45 DOUBLE PRECISION H(IDIM,M),E(M),AID(5501),w3(5501) 46 DOUBLE PRECISION ZR(IDIM,M),ZI(IDIM,M),W1(M),W2(2,M) 47 DOUBLE PRECISION EIGSUM,TRACE 48c ============================================================== 49 COMMON /AIDIAG/AID 50c ============================================================== 51 integer :: i, j, ier 52 53c IF(M.GT.501)WRITE(*,*)' ERROR IN HERMDP ******* M>501' 54c IF(M.GT.501)GO TO 999 55c write(*,*)' welcome to hermdp.f ......' 56 IF(M.GT.5501)WRITE(*,*)' ERROR IN HERMDP ******* M>5501' 57 IF(M.GT.5501)GO TO 999 58C FIND TRACE OF MATRIX, to test sum of eigenvalues. 59 TRACE=0.D0 60 DO 10 I=1,M 6110 TRACE=TRACE+H(I,I) 62C INITIALIZE ZR,ZI. must set zr to identity. 63 DO 1 I=1,M 64 DO 2 J=1,M 652 ZR(I,J)=0.D0 661 ZR(I,I)=1.D0 67c ============================================================== 68c first call. 69 CALL DPHTRI(IDIM,M,H,E,W1,w3,W2) 70 IF(M.LT.300)GO TO 3 71 WRITE(*,*)' FIRST CALL FINISHED' 723 continue 73c ============================================================== 74c second call. 75 CALL DPIMTQ(IDIM,M,E,W1,ZR,IER) 76 IF(M.LT.300)GO TO 4 77 WRITE(*,*)' SECOND CALL FINISHED' 784 continue 79c ============================================================== 80c third call: for eigenvectors. 81 if(icall.eq.3)then 82 CALL DPHTBK(IDIM,M,H,W2,M,ZR,ZI) 83 IF(M.GT.300)WRITE(*,*)' THIRD CALL FINISHED' 84 else 85 write(*,*)' **caution** in hermdp, eigenvcs not desired' 86 end if 87 IF(M.GT.300.and.icall.eq.3)WRITE(*,*)' THIRD CALL FINISHED' 88c ============================================================== 89C SUM EIGENVALUES AND COMPARE WITH TRACE. 90 EIGSUM=0.D0 91 DO 11 I=1,M 9211 EIGSUM=EIGSUM+E(I) 93 IF(100.D0*DABS(EIGSUM-TRACE).LT.DABS(TRACE)*1.D-9)GOTO 12 94 WRITE(*,*)' ERROR: hERMDP SUM OF EIGENVALUES NOT EQUAL TO TRACE' 95 WRITE(*,*)' EIGSUM=',EIGSUM,' TRACE=',TRACE 96 IF(DABS(EIGSUM-TRACE).GT.DABS(TRACE)*1.D-08) 97 X WRITE(*,*)' BAD ERROR IN hERMDP ******* !!!!!! ******' 9812 IF(IER.NE.0) WRITE(*,*)' ERROR IN HERMDP IER IN HERMDP=',IER 99c ============================================================== 100 RETURN 101999 STOP 102 END 103c ============================================================== 104C 105 SUBROUTINE DPHTRI(NM,N,AR,D,E,E2,TAU) 106C 107 INTEGER I,J,K,L,N,II,NM,JP1 108 DOUBLE PRECISION AR(NM,N),D(N),E(N),E2(N),TAU(2,N) 109c DOUBLE PRECISION F,FI,G,GI,H,HH,SI,SCALE,AID(501) 110 DOUBLE PRECISION F,FI,G,GI,H,HH,SI,SCALE,AID(5501) 111 DOUBLE PRECISION DSQRT,CDABS,DABS 112 double precision ajunk 113 COMMON /AIDIAG/AID 114C COMPLEX*16 DCMPLX 115C 116C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF 117C THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968) 118C BY MARTIN, REINSCH, AND WILKINSON. 119C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). 120C 121C THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX 122C TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING 123C UNITARY SIMILARITY TRANSFORMATIONS. 124C 125C ON INPUT- 126C 127C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 128C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 129C DIMENSION STATEMENT, 130C 131C N IS THE ORDER OF THE MATRIX, 132C 133C AR AND AI CONTAIN THE REAL AND AIMAGINARY PARTS, 134C RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX. 135C ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. 136C 137C ON OUTPUT- 138C 139C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- 140C FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER 141C TRIANGLES. THEIR STRICT UPPER TRIANGLES AND THE 142C DIAGONAL OF AR ARE UNALTERED, 143C 144C D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX, 145C 146C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL 147C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO, 148C 149C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. 150C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED, 151C 152C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. 153C 154C ARITHMETIC IS REAL EXCEPT FOR THE USE OF THE SUBROUTINES 155C CABS AND CMPLX IN COMPUTING COMPLEX ABSOLUTE VALUES. 156C 157C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, 158C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 159C 160C ------------------------------------------------------------------ 161C 162 TAU(1,N) = 1.0D0 163 TAU(2,N) = 0.0D0 164c DO 324 I=1,501 165 DO 324 I=1,5501 166324 AID(I)=0.D0 167C 168C 169 DO 100 I = 1, N 170 100 D(I) = AR(I,I) 171C 172C ********** FOR I=N STEP -1 UNTIL 1 DO -- ********** 173 DO 300 II = 1, N 174 I = N + 1 - II 175 L = I - 1 176 H = 0.0D0 177 SCALE = 0.0D0 178 IF (L .LT. 1) GO TO 130 179C ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) ********** 180 DO 120 K = 1, L 181 120 SCALE = SCALE + DABS(AR(I,K)) + DABS(AR(K,I)) 182C 183 IF (SCALE .NE. 0.0) GO TO 140 184 TAU(1,L) = 1.0D0 185 TAU(2,L) = 0.0D0 186 130 E(I) = 0.0D0 187 E2(I) = 0.0D0 188 GO TO 290 189C 190 140 DO 150 K = 1, L 191 AR(I,K) = AR(I,K) / SCALE 192 AR(K,I) = AR(K,I) / SCALE 193 H = H + AR(I,K) * AR(I,K) + AR(K,I) * AR(K,I) 194 150 CONTINUE 195C 196 E2(I) = SCALE * SCALE * H 197 G = DSQRT(H) 198 E(I) = SCALE * G 199 F=DSQRT(AR(I,L)**2+AR(L,I)**2) 200C F = CDABS(DCMPLX(AR(I,L),AI(I,L))) 201C ********** FORM NEXT DIAGONAL ELEMENT OF MATRIX T ********** 202 IF (F .EQ. 0.0D0) GO TO 160 203 TAU(1,L) = (AR(L,I) * TAU(2,I) - AR(I,L) * TAU(1,I)) / F 204 SI = (AR(I,L) * TAU(2,I) + AR(L,I) * TAU(1,I)) / F 205 H = H + F * G 206 G = 1.0D0 + G / F 207 AR(I,L) = G * AR(I,L) 208 AR(L,I) = G * AR(L,I) 209 IF (L .EQ. 1) GO TO 270 210 GO TO 170 211 160 TAU(1,L) = -TAU(1,I) 212 SI = TAU(2,I) 213 AR(I,L) = G 214 170 F = 0.0D0 215C 216 DO 240 J = 1, L 217 G = 0.0D0 218 GI = 0.0D0 219C ********** FORM ELEMENT OF A*U ********** 220 DO 180 K = 1, J 221 AJUNK=AR(K,J) 222 IF(K.EQ.J)AJUNK=AID(K) 223 G = G + AR(J,K) * AR(I,K) + AJUNK * AR(K,I) 224 GI = GI - AR(J,K) * AR(K,I) + AJUNK * AR(I,K) 225 180 CONTINUE 226C 227 JP1 = J + 1 228 IF (L .LT. JP1) GO TO 220 229C 230 DO 200 K = JP1, L 231 G = G + AR(K,J) * AR(I,K) - AR(J,K) * AR(K,I) 232 GI = GI - AR(K,J) * AR(K,I)- AR(J,K) * AR(I,K) 233 200 CONTINUE 234C ********** FORM ELEMENT OF P ********** 235 220 E(J) = G / H 236 TAU(2,J) = GI / H 237 F = F + E(J) * AR(I,J) - TAU(2,J) * AR(J,I) 238 240 CONTINUE 239C 240 HH = F / (H + H) 241C ********** FORM REDUCED A ********** 242 DO 260 J = 1, L 243 F = AR(I,J) 244 G = E(J) - HH * F 245 E(J) = G 246 FI = -AR(J,I) 247 GI = TAU(2,J) - HH * FI 248 TAU(2,J) = -GI 249C 250 DO 260 K = 1, J 251 AJUNK=AR(K,J) 252 IF(J.EQ.K)AJUNK=AID(K) 253 AR(J,K) = AR(J,K) - F * E(K) - G * AR(I,K) 254 X + FI * TAU(2,K) + GI * AR(K,I) 255 IF(K.EQ.J)GO TO 17 256 AR(K,J) = AR(K,J) - F * TAU(2,K) - G * AR(K,I) 257 X - FI * E(K) - GI * AR(I,K) 258 GO TO 260 25917 AID(K)=AID(K) - F * TAU(2,K) - G * AR(K,I) 260 X - FI * E(K) - GI * AR(I,K) 261 260 CONTINUE 262C 263 270 DO 280 K = 1, L 264 AR(I,K) = SCALE * AR(I,K) 265 AR(K,I) = SCALE * AR(K,I) 266 280 CONTINUE 267C 268 TAU(2,L) = -SI 269 290 HH = D(I) 270 D(I) = AR(I,I) 271 AR(I,I) = HH 272 AID(I) = SCALE*DSQRT(H) 273 300 CONTINUE 274C 275 RETURN 276C ********** LAST CARD OF DPHTRI ********** 277 END 278c ============================================================== 279 SUBROUTINE DPIMTQ(NM,N,D,E,Z,IERR) 280C 281 INTEGER I,J,K,L,M,N,II,NM,MML,IERR 282 DOUBLE PRECISION D(N),E(N),Z(NM,N) 283 DOUBLE PRECISION B,C,F,G,P,R,S,MACHEP 284 DOUBLE PRECISION DSQRT,DABS,DSIGN 285C 286C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE DPIMTQ, 287C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, 288C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. 289C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). 290C 291C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS 292C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. 293C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO 294C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS 295C FULL MATRIX TO TRIDIAGONAL FORM. 296C 297C ON INPUT- 298C 299C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 300C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 301C DIMENSION STATEMENT, 302C 303C N IS THE ORDER OF THE MATRIX, 304C 305C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX, 306C 307C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX 308C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY, 309C 310C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE 311C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS 312C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN 313C THE IDENTITY MATRIX. 314C 315C ON OUTPUT- 316C 317C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN 318C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT 319C UNORDERED FOR INDICES 1,2,...,IERR-1, 320C 321C E HAS BEEN DESTROYED, 322C 323C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC 324C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, 325C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED 326C EIGENVALUES, 327C 328C IERR IS SET TO 329C ZERO FOR NORMAL RETURN, 330C J IF THE J-TH EIGENVALUE HAS NOT BEEN 331C DETERMINED AFTER 30 ITERATIONS. 332C 333C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, 334C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 335C 336C ------------------------------------------------------------------ 337C 338C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING 339C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. 340C MACHEP=2**(-26) SP. 2**(-56) DP. (26 AND 56 ARE MAXIMUMS!) 341C WE USE 54 FOR SAFETY. (1.D0+MACHEP.GT.1.D0) DEFINES MINIMUM MACHEP. 342C ********** 343 MACHEP=2.D0**(-54) 344C 345 IERR = 0 346 IF (N .EQ. 1) GO TO 1001 347C 348 DO 100 I = 2, N 349 100 E(I-1) = E(I) 350C 351 E(N) = 0.0D0 352C 353 DO 240 L = 1, N 354 J = 0 355C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ********** 356 105 DO 110 M = L, N 357 IF (M .EQ. N) GO TO 120 358 IF (DABS(E(M)) .LE. MACHEP * (DABS(D(M)) + DABS(D(M+1)))) 359 X GO TO 120 360 110 CONTINUE 361C 362 120 P = D(L) 363 IF (M .EQ. L) GO TO 240 364 IF (J .EQ. 30) GO TO 1000 365 J = J + 1 366C ********** FORM SHIFT ********** 367 G = (D(L+1) - P) / (2.0D0 * E(L)) 368 R = DSQRT(G*G+1.0D0) 369 G = D(M) - P + E(L) / (G + DSIGN(R,G)) 370 S = 1.0D0 371 C = 1.0D0 372 P = 0.0D0 373 MML = M - L 374C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** 375 DO 200 II = 1, MML 376 I = M - II 377 F = S * E(I) 378 B = C * E(I) 379 IF (DABS(F) .LT. DABS(G)) GO TO 150 380 C = G / F 381 R =DSQRT(C*C+1.0D0) 382 E(I+1) = F * R 383 S = 1.0D0 / R 384 C = C * S 385 GO TO 160 386 150 S = F / G 387 R = DSQRT(S*S+1.0D0) 388 E(I+1) = G * R 389 C = 1.0D0 / R 390 S = S * C 391 160 G = D(I+1) - P 392 R = (D(I) - G) * S + 2.0D0 * C * B 393 P = S * R 394 D(I+1) = G + P 395 G = C * R - B 396C ********** FORM VECTOR ********** 397 DO 180 K = 1, N 398 F = Z(K,I+1) 399 Z(K,I+1) = S * Z(K,I) + C * F 400 Z(K,I) = C * Z(K,I) - S * F 401 180 CONTINUE 402C 403 200 CONTINUE 404C 405 D(L) = D(L) - P 406 E(L) = G 407 E(M) = 0.0D0 408 GO TO 105 409 240 CONTINUE 410C ********** ORDER EIGENVALUES AND EIGENVECTORS ********** 411 DO 300 II = 2, N 412 I = II - 1 413 K = I 414 P = D(I) 415C 416 DO 260 J = II, N 417 IF (D(J) .GE. P) GO TO 260 418 K = J 419 P = D(J) 420 260 CONTINUE 421C 422 IF (K .EQ. I) GO TO 300 423 D(K) = D(I) 424 D(I) = P 425C 426 DO 280 J = 1, N 427 P = Z(J,I) 428 Z(J,I) = Z(J,K) 429 Z(J,K) = P 430 280 CONTINUE 431C 432 300 CONTINUE 433C 434 GO TO 1001 435C ********** SET ERROR -- NO CONVERGENCE TO AN 436C EIGENVALUE AFTER 30 ITERATIONS ********** 437 1000 IERR = L 438 1001 RETURN 439C ********** LAST CARD OF DPIMTQ ********** 440 END 441c ============================================================== 442 SUBROUTINE DPHTBK(NM,N,AR,TAU,M,ZR,ZI) 443C 444 INTEGER I,J,K,L,M,N,NM 445 DOUBLE PRECISION AR(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M) 446c DOUBLE PRECISION H,S,SI,AID(501) 447 DOUBLE PRECISION H,S,SI,AID(5501) 448 COMMON /AIDIAG/AID 449C 450C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF 451C THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968) 452C BY MARTIN, REINSCH, AND WILKINSON. 453C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). 454C 455C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN 456C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING 457C REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY DPHTRI. 458C 459C ON INPUT- 460C 461C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 462C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 463C DIMENSION STATEMENT, 464C 465C N IS THE ORDER OF THE MATRIX, 466C 467C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- 468C FORMATIONS USED IN THE REDUCTION BY DPHTRI IN THEIR 469C FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR. 470C 471C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS, 472C 473C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED, 474C 475C ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED 476C IN ITS FIRST M COLUMNS. 477C 478C ON OUTPUT- 479C 480C ZR AND ZI CONTAIN THE REAL AND AIMAGINARY PARTS, 481C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS 482C IN THEIR FIRST M COLUMNS. 483C 484C NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR 485C IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED. 486C 487C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, 488C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 489C 490C ------------------------------------------------------------------ 491 IF(M .EQ. 0 ) GOTO 200 492C ********** TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC 493C TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN 494C TRIDIAGONAL MATRIX. ********** 495 DO 50 K = 1, N 496C 497 DO 50 J = 1, M 498 ZI(K,J) = - ZR(K,J) * TAU(2,K) 499 ZR(K,J) = ZR(K,J) * TAU(1,K) 500 50 CONTINUE 501C 502 IF (N .EQ. 1) GO TO 200 503C ********** RECOVER AND APPLY THE HOUSEHOLDER MATRICES ********** 504 DO 140 I = 2, N 505 L = I - 1 506 H = AID(I) 507 IF (H .EQ. 0.0D0) GO TO 140 508C 509 DO 130 J = 1, M 510 S = 0.0D0 511 SI = 0.0D0 512C 513 DO 110 K = 1, L 514 S = S + AR(I,K) * ZR(K,J) - AR(K,I) * ZI(K,J) 515 SI = SI + AR(I,K) * ZI(K,J) + AR(K,I) * ZR(K,J) 516 110 CONTINUE 517C 518 S=(S/H)/H 519 SI=(SI/H)/H 520C 521 DO 120 K = 1, L 522 ZR(K,J) = ZR(K,J) - S * AR(I,K) - SI * AR(K,I) 523 ZI(K,J) = ZI(K,J) - SI * AR(I,K) + S * AR(K,I) 524 120 CONTINUE 525C 526 130 CONTINUE 527C 528 140 CONTINUE 529C 530 200 RETURN 531C ********** LAST CARD OF DPHTBK ********** 532 END 533