home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / libs / matlab.lzh / MATLAB / MATLAB.LZH / Source / MatLab / HTRIDI.FOR < prev    next >
Encoding:
Text File  |  1991-04-13  |  4.5 KB  |  157 lines
  1.       SUBROUTINE HTRIDI (NM, N, AR, AI, D, E, E2, TAU)
  2.       IMPLICIT NONE
  3. C
  4. C A TRANSLATION OF A COMPLEX ANALOGUE OF THE ALGOL PROCEDURE TRED1,
  5. C   NUM. MATH. 11, 181-195 (1968) BY MARTIN, REINSCH, AND WILKINSON.
  6. C   HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226 (1971).
  7. C
  8. C REDUCES A COMPLEX HERMITIAN MATRIX TO A REAL SYMMETRIC TRIDIAGONAL MATRIX
  9. C   USING UNITARY SIMILARITY TRANSFORMATIONS.
  10. C
  11. C ON INPUT:
  12. C   NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  13. C      AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  14. C
  15. C   N IS THE ORDER OF THE MATRIX.
  16. C
  17. C   AR & AI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE
  18. C      COMPLEX HERMITIAN INPUT MATRIX.  ONLY THE LOWER TRIANGLE OF THE
  19. C      MATRIX NEED BE SUPPLIED.
  20. C
  21. C ON OUTPUT:
  22. C   AR & AI CONTAIN INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN
  23. C      THE REDUCTION IN THEIR FULL LOWER TRIANGLES.  THEIR STRICT UPPER
  24. C      TRIANGLES AND THE DIAGONAL OF AR ARE UNALTERED.
  25. C
  26. C   D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX.
  27. C
  28. C   E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX IN ITS
  29. C      LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  30. C
  31. C   E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  32. C      E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  33. C
  34. C   TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
  35. C
  36. C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  37. C   APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  38. C
  39.       INTEGER NM, N
  40.       DOUBLE PRECISION AR(NM,N), AI(NM,N), D(N), E(N), E2(N), TAU(2,N)
  41. C
  42.       INTEGER I, J, K, L, II, JP1
  43.       DOUBLE PRECISION F, G, H, FI, GI, HH, SI, SCALE
  44. C
  45.       DOUBLE PRECISION FLOP, PYTHAG
  46. C
  47. C
  48.       TAU(1,N) = 1.0D0
  49.       TAU(2,N) = 0.0D0
  50.       DO 100 I = 1, N
  51.         D(I) = AR(I,I)
  52. 100   CONTINUE
  53. C
  54. C ***      FOR I = N STEP -1 UNTIL 1 DO
  55.       DO 300 II = 1, N
  56.         I = N+1-II
  57.         L = I-1
  58.         H = 0.0D0
  59.         SCALE = 0.0D0
  60.         IF (L.LT.1) GO TO 130
  61. C
  62. C ***      SCALE ROW (ALGOL TOL THEN NOT NEEDED)
  63.         DO 120 K = 1, L
  64.           SCALE = FLOP (SCALE+DABS (AR(I,K))+DABS (AI(I,K)))
  65. 120     CONTINUE
  66.         IF (SCALE.NE.0.0D0) GO TO 140
  67.         TAU(1,L) = 1.0D0
  68.         TAU(2,L) = 0.0D0
  69. 130     CONTINUE
  70.         E(I) = 0.0D0
  71.         E2(I) = 0.0D0
  72.         GO TO 290
  73. C
  74. 140     CONTINUE
  75.         DO 150 K = 1, L
  76.           AR(I,K) = FLOP (AR(I,K)/SCALE)
  77.           AI(I,K) = FLOP (AI(I,K)/SCALE)
  78.           H = FLOP (H+AR(I,K)*AR(I,K)+AI(I,K)*AI(I,K))
  79. 150     CONTINUE
  80.         E2(I) = FLOP (SCALE*SCALE*H)
  81.         G = FLOP (DSQRT (H))
  82.         E(I) = FLOP (SCALE*G)
  83.         F = PYTHAG (AR(I,L), AI(I,L))
  84. C
  85. C ***      FORM NEXT DIAGONAL ELEMENT OF MATRIX T
  86.         IF (F.EQ.0.0D0) GO TO 160
  87.         TAU(1,L) = FLOP ((AI(I,L)*TAU(2,I)-AR(I,L)*TAU(1,I))/F)
  88.         SI = FLOP ((AR(I,L)*TAU(2,I)+AI(I,L)*TAU(1,I))/F)
  89.         H = FLOP (H+F*G)
  90.         G = FLOP (1.0D0+G/F)
  91.         AR(I,L) = FLOP (G*AR(I,L))
  92.         AI(I,L) = FLOP (G*AI(I,L))
  93.         IF (L.EQ.1) GO TO 270
  94.         GO TO 170
  95. C
  96. 160     CONTINUE
  97.         TAU(1,L) = -TAU(1,I)
  98.         SI = TAU(2,I)
  99.         AR(I,L) = G
  100. 170     CONTINUE
  101.         F = 0.0D0
  102.         DO 240 J = 1, L
  103.           G = 0.0D0
  104.           GI = 0.0D0
  105. C
  106. C ***      FORM ELEMENT OF A*U
  107.           DO 180 K = 1, J
  108.             G = FLOP (G+AR(J,K)*AR(I,K)+AI(J,K)*AI(I,K))
  109.             GI = FLOP (GI-AR(J,K)*AI(I,K)+AI(J,K)*AR(I,K))
  110. 180       CONTINUE
  111.           JP1 = J+1
  112.           IF (L.LT.JP1) GO TO 220
  113.           DO 200 K = JP1, L
  114.             G = FLOP (G+AR(K,J)*AR(I,K)-AI(K,J)*AI(I,K))
  115.             GI = FLOP (GI-AR(K,J)*AI(I,K)-AI(K,J)*AR(I,K))
  116. 200       CONTINUE
  117. C
  118. C ***      FORM ELEMENT OF P
  119. 220       CONTINUE
  120.           E(J) = FLOP (G/H)
  121.           TAU(2,J) = FLOP (GI/H)
  122.           F = FLOP (F+E(J)*AR(I,J)-TAU(2,J)*AI(I,J))
  123. 240     CONTINUE
  124.         HH = FLOP (F/(H+H))
  125. C
  126. C ***      FORM REDUCED A
  127.         DO 261 J = 1, L
  128.           F = AR(I,J)
  129.           G = FLOP (E(J)-HH*F)
  130.           E(J) = G
  131.           FI = -AI(I,J)
  132.           GI = FLOP (TAU(2,J)-HH*FI)
  133.           TAU(2,J) = -GI
  134.           DO 260 K = 1, J
  135.             AR(J,K) = FLOP (AR(J,K)-F*E(K)-G*AR(I,K)
  136.      .                      +FI*TAU(2,K)+GI*AI(I,K))
  137.             AI(J,K) = FLOP (AI(J,K)-F*TAU(2,K)-G*AI(I,K)
  138.      .                      -FI*E(K)-GI*AR(I,K))
  139. 260       CONTINUE
  140. 261     CONTINUE
  141. C
  142. 270     CONTINUE
  143.         DO 280 K = 1, L
  144.           AR(I,K) = FLOP (SCALE*AR(I,K))
  145.           AI(I,K) = FLOP (SCALE*AI(I,K))
  146. 280     CONTINUE
  147.         TAU(2,L) = -SI
  148. 290     CONTINUE
  149.         HH = D(I)
  150.         D(I) = AR(I,I)
  151.         AR(I,I) = HH
  152.         AI(I,I) = FLOP (SCALE*DSQRT (H))
  153. 300   CONTINUE
  154. C
  155.       RETURN
  156.       END
  157.