home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tests / treds.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  283 lines
  1.       SUBROUTINE tred1(NM,N,A,D,E,E2)
  2. C
  3.       INTEGER I,J,K,L,N,II,NM,JP1
  4.       REAL A(NM,N),D(N),E(N),E2(N)
  5.       REAL c,F,G,H,SCALE
  6. C
  7. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1,
  8. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  9. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  10. C
  11. C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX
  12. C     TO A SYMMETRIC TRIDIAGONAL MATRIX USING
  13. C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  14. C
  15. C     ON INPUT
  16. C
  17. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  18. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  19. C          DIMENSION STATEMENT.
  20. C
  21. C        N IS THE ORDER OF THE MATRIX.
  22. C
  23. C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
  24. C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
  25. C
  26. C     ON OUTPUT
  27. C
  28. C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
  29. C          FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
  30. C          TRIANGLE.  THE FULL UPPER TRIANGLE OF A IS UNALTERED.
  31. C
  32. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
  33. C
  34. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  35. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  36. C
  37. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  38. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  39. C
  40. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  41. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  42. C
  43. C     ------------------------------------------------------------------
  44. C
  45.       DO 100 I = 1, N
  46.   100 D(I) = A(n,I)
  47.       do 110 i = 1, n
  48.   110 a(n,i) = a(i,i)
  49. C     .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
  50.       DO 300 II = 1, N
  51.          I = N + 1 - II
  52.          L = I - 1
  53.          H = 0.0E0
  54.          SCALE = 0.0E0
  55.          IF (L .LT. 1) GO TO 130
  56. C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
  57.          DO 120 K = 1, L
  58.   120    SCALE = SCALE + ABS(d(k))
  59. C
  60.          IF (SCALE .NE. 0.0E0) GO TO 140
  61.   130    E(I) = 0.0E0
  62.          E2(I) = 0.0E0
  63.          if (l .lt. 1) go to 300
  64.          do 135 k = 1, l
  65.             f = d(k)
  66.             d(k) = a(l,k)
  67.             a(l,k) = a(i,k)
  68.             a(i,k) = f
  69.   135    continue 
  70.          go to 300
  71. C
  72.   140    DO 150 K = 1, L
  73.             d(k) = d(k) / SCALE
  74.             H = H + d(k) * d(k)
  75.   150    CONTINUE
  76. C
  77.          E2(I) = SCALE * SCALE * H
  78.          f = d(l)
  79.          G = -SIGN(SQRT(H),F)
  80.          E(I) = SCALE * G
  81.          H = H - F * G
  82.          d(l) = f - g
  83.          if (l .gt. 1) go to 160
  84.          f = d(1)
  85.          d(1) = a(1,1)
  86.          a(1,1) = a(2,1)
  87.          a(2,1) = f * scale
  88.          go to 300
  89.   160    continue
  90.          F = 0.0E0
  91. C
  92.          do 170 k = 1, l
  93.   170    e(k) = 0.0e0
  94.          DO 240 J = 1, L
  95.             c = d(j)
  96.             g = e(j) + a(j,j) * c
  97.             JP1 = J + 1
  98.             IF (L .LT. JP1) GO TO 220
  99. C
  100.             DO 200 K = JP1, L
  101.                e(k) = e(k) + a(k,j)*c
  102.                g = g + a(k,j) * d(k)
  103.   200       continue
  104. C     .......... FORM ELEMENT OF P ..........
  105.   220       E(J) = G / H
  106.             f = f + e(j) * c
  107.   240    CONTINUE
  108. C
  109.          H = F / (H + H)
  110. C     .......... FORM REDUCED A ..........
  111.          do 250 k = 1, l
  112.   250    e(k) = e(k) - h * d(k)          
  113.          do 280 j = 1, l
  114.             f = d(j)
  115.             g = e(j)
  116. C
  117.             do 260 k = j, l
  118.   260       a(k,j) = a(k,j) - g * d(k) - f * e(k)
  119.             d(j) = a(l,j)
  120.             a(l,j) = a(i,j)
  121.             a(i,j) = f * scale
  122.   280    continue
  123. C
  124.   300 CONTINUE
  125. C
  126.       RETURN
  127.       END
  128.       SUBROUTINE TRED2(NM,N,A,D,E,Z)
  129. C
  130.       INTEGER I,J,K,L,N,II,NM,JP1
  131.       REAL A(NM,N),D(N),E(N),Z(NM,N)
  132.       REAL c,F,G,H,HH,SCALE
  133. C
  134. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
  135. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  136. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  137. C
  138. C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A
  139. C     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
  140. C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  141. C
  142. C     ON INPUT
  143. C
  144. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  145. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  146. C          DIMENSION STATEMENT.
  147. C
  148. C        N IS THE ORDER OF THE MATRIX.
  149. C
  150. C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
  151. C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
  152. C
  153. C     ON OUTPUT
  154. C
  155. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
  156. C
  157. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  158. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
  159. C
  160. C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
  161. C          PRODUCED IN THE REDUCTION.
  162. C
  163. C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
  164. C
  165. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  166. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  167. C
  168. C     ------------------------------------------------------------------
  169. C
  170.       DO 100 I = 1, N
  171. C
  172.          DO 80 J = i, n
  173.    80    z(j,i) = a(j,i)
  174.          d(i) = a(n,i)
  175.   100 CONTINUE
  176. C
  177.       IF (N .EQ. 1) GO TO 320
  178. C     .......... FOR I=N STEP -1 UNTIL 2 DO -- ..........
  179.       DO 300 II = 2, N
  180.          I = N + 2 - II
  181.          L = I - 1
  182.          H = 0.0E0
  183.          SCALE = 0.0E0
  184.          IF (L .LT. 2) GO TO 130
  185. C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
  186.          DO 120 K = 1, L
  187.   120    SCALE = SCALE + ABS(d(k))
  188. C
  189.          IF (SCALE .NE. 0.0E0) GO TO 140
  190.   130    E(I) = d(l)
  191.          h = 0.0e0
  192.          do 135 j = 1,l
  193.             z(i,j) = 0.0e0
  194.             z(j,i) = d(j)
  195.             d(j) = z(l,j)
  196.   135    continue
  197.          GO TO 290
  198. C
  199.   140    DO 150 K = 1, L
  200.             d(k) = d(k) / SCALE
  201.             H = H + d(k) * d(k)
  202.   150    CONTINUE
  203. C
  204.          F = d(l)
  205.          G = -SIGN(SQRT(H),F)
  206.          E(I) = SCALE * G
  207.          H = H - F * G
  208.          d(l) = F - G
  209.          F = 0.0E0
  210. C
  211.          do 170 k = 1, l
  212.   170    e(k) = 0.0e0           
  213.          DO 240 J = 1, L
  214.             c = d(j)
  215.             z(j,i) = c
  216.             g = e(j) + z(j,j)*c
  217.             JP1 = J + 1
  218.             IF (L .LT. JP1) GO TO 220
  219. C
  220.             DO 200 K = JP1, L
  221.                e(k) = e(k) + z(k,j) * c
  222.                G = G + Z(K,J) * d(k)
  223.   200       continue
  224. C     .......... FORM ELEMENT OF P ..........
  225.   220       E(J) = G / H
  226.             F = F + E(J) * c
  227.   240    CONTINUE
  228. C
  229.          HH = F / (H + H)
  230. C     .......... FORM REDUCED A ..........
  231.          do 250 k = 1, l
  232.   250    e(k) = e(k) - hh * d(k)          
  233.          DO 280 J = 1, L
  234.             F = d(j)
  235.             G = E(J)
  236. C
  237.             DO 260 K = j, l
  238.   260       Z(k,j) = Z(k,j) - g * d(k) - f * e(k)
  239.             d(j) = z(l,j)
  240.             z(i,j) = 0.0e0
  241.   280    continue
  242. C
  243.   290    D(I) = H
  244.   300 CONTINUE
  245. C
  246.   320 e(1) = 0.0E0
  247. C     .......... ACCUMULATION OF TRANSFORMATION MATRICES ..........
  248.       if (n .eq. 1) go to 510
  249.       DO 500 I = 2, N
  250.          L = I - 1
  251.          z(n,l) = z(l,l) 
  252.          z(l,l) = 1.0e0
  253.          h = d(i)
  254.          IF (h .EQ. 0.0E0) GO TO 380
  255. C
  256.          do 330 k = 1,l
  257.   330    d(k) = z(k,i) / h         
  258.          DO 360 J = 1, L
  259.             G = 0.0E0
  260. C
  261.             DO 340 K = 1, L
  262.   340       G = G + Z(k,i) * Z(K,J)
  263. C
  264.             DO 360 K = 1, L
  265.                Z(K,J) = Z(K,J) - G * d(k)
  266.   360    CONTINUE
  267. C
  268.   380    do 400 k = 1, l
  269.   400    Z(k,I) = 0.0E0
  270. C
  271.   500 CONTINUE
  272. C
  273.   510 do 520 i = 1,n
  274.          d(i) = z(n,i)
  275.          z(n,i) = 0.0e0
  276.   520 continue
  277.       z(n,n) = 1.0e0
  278. c
  279.       RETURN
  280.       END
  281.  
  282.  
  283.