PageRenderTime 59ms CodeModel.GetById 19ms RepoModel.GetById 1ms app.codeStats 0ms

/UAFWLIS2.4/slatec/tsturm.f

http://ua-fwlis.googlecode.com/
FORTRAN Legacy | 405 lines | 222 code | 0 blank | 183 comment | 0 complexity | a1cb637d3bbbf0fb44ceb5f05cf0123f MD5 | raw file
Possible License(s): GPL-3.0
  1. *DECK TSTURM
  2. SUBROUTINE TSTURM (NM, N, EPS1, D, E, E2, LB, UB, MM, M, W, Z,
  3. + IERR, RV1, RV2, RV3, RV4, RV5, RV6)
  4. C***BEGIN PROLOGUE TSTURM
  5. C***PURPOSE Find those eigenvalues of a symmetric tridiagonal matrix
  6. C in a given interval and their associated eigenvectors by
  7. C Sturm sequencing.
  8. C***LIBRARY SLATEC (EISPACK)
  9. C***CATEGORY D4A5, D4C2A
  10. C***TYPE SINGLE PRECISION (TSTURM-S)
  11. C***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
  12. C***AUTHOR Smith, B. T., et al.
  13. C***DESCRIPTION
  14. C
  15. C This subroutine finds those eigenvalues of a TRIDIAGONAL
  16. C SYMMETRIC matrix which lie in a specified interval and their
  17. C associated eigenvectors, using bisection and inverse iteration.
  18. C
  19. C On Input
  20. C
  21. C NM must be set to the row dimension of the two-dimensional
  22. C array parameter, Z, as declared in the calling program
  23. C dimension statement. NM is an INTEGER variable.
  24. C
  25. C N is the order of the matrix. N is an INTEGER variable.
  26. C N must be less than or equal to NM.
  27. C
  28. C EPS1 is an absolute error tolerance for the computed eigen-
  29. C values. It should be chosen so that the accuracy of these
  30. C eigenvalues is commensurate with relative perturbations of
  31. C the order of the relative machine precision in the matrix
  32. C elements. If the input EPS1 is non-positive, it is reset
  33. C for each submatrix to a default value, namely, minus the
  34. C product of the relative machine precision and the 1-norm of
  35. C the submatrix. EPS1 is a REAL variable.
  36. C
  37. C D contains the diagonal elements of the symmetric tridiagonal
  38. C matrix. D is a one-dimensional REAL array, dimensioned D(N).
  39. C
  40. C E contains the subdiagonal elements of the symmetric
  41. C tridiagonal matrix in its last N-1 positions. E(1) is
  42. C arbitrary. E is a one-dimensional REAL array, dimensioned
  43. C E(N).
  44. C
  45. C E2 contains the squares of the corresponding elements of E.
  46. C E2(1) is arbitrary. E2 is a one-dimensional REAL array,
  47. C dimensioned E2(N).
  48. C
  49. C LB and UB define the interval to be searched for eigenvalues.
  50. C If LB is not less than UB, no eigenvalues will be found.
  51. C LB and UB are REAL variables.
  52. C
  53. C MM should be set to an upper bound for the number of
  54. C eigenvalues in the interval. MM is an INTEGER variable.
  55. C WARNING - If more than MM eigenvalues are determined to lie
  56. C in the interval, an error return is made with no values or
  57. C vectors found.
  58. C
  59. C On Output
  60. C
  61. C EPS1 is unaltered unless it has been reset to its
  62. C (last) default value.
  63. C
  64. C D and E are unaltered.
  65. C
  66. C Elements of E2, corresponding to elements of E regarded as
  67. C negligible, have been replaced by zero causing the matrix to
  68. C split into a direct sum of submatrices. E2(1) is also set
  69. C to zero.
  70. C
  71. C M is the number of eigenvalues determined to lie in (LB,UB).
  72. C M is an INTEGER variable.
  73. C
  74. C W contains the M eigenvalues in ascending order if the matrix
  75. C does not split. If the matrix splits, the eigenvalues are
  76. C in ascending order for each submatrix. If a vector error
  77. C exit is made, W contains those values already found. W is a
  78. C one-dimensional REAL array, dimensioned W(MM).
  79. C
  80. C Z contains the associated set of orthonormal eigenvectors.
  81. C If an error exit is made, Z contains those vectors already
  82. C found. Z is a one-dimensional REAL array, dimensioned
  83. C Z(NM,MM).
  84. C
  85. C IERR is an INTEGER flag set to
  86. C Zero for normal return,
  87. C 3*N+1 if M exceeds MM no eigenvalues or eigenvectors
  88. C are computed,
  89. C 4*N+J if the eigenvector corresponding to the J-th
  90. C eigenvalue fails to converge in 5 iterations, then
  91. C the eigenvalues and eigenvectors in W and Z should
  92. C be correct for indices 1, 2, ..., J-1.
  93. C
  94. C RV1, RV2, RV3, RV4, RV5, and RV6 are temporary storage arrays,
  95. C dimensioned RV1(N), RV2(N), RV3(N), RV4(N), RV5(N), and
  96. C RV6(N).
  97. C
  98. C The ALGOL procedure STURMCNT contained in TRISTURM
  99. C appears in TSTURM in-line.
  100. C
  101. C Questions and comments should be directed to B. S. Garbow,
  102. C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  103. C ------------------------------------------------------------------
  104. C
  105. C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
  106. C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
  107. C system Routines - EISPACK Guide, Springer-Verlag,
  108. C 1976.
  109. C***ROUTINES CALLED R1MACH
  110. C***REVISION HISTORY (YYMMDD)
  111. C 760101 DATE WRITTEN
  112. C 890531 Changed all specific intrinsics to generic. (WRB)
  113. C 890531 REVISION DATE from Version 3.2
  114. C 891214 Prologue converted to Version 4.0 format. (BAB)
  115. C 920501 Reformatted the REFERENCES section. (WRB)
  116. C***END PROLOGUE TSTURM
  117. C
  118. INTEGER I,J,K,M,N,P,Q,R,S,II,IP,JJ,MM,M1,M2,NM,ITS
  119. INTEGER IERR,GROUP,ISTURM
  120. REAL D(*),E(*),E2(*),W(*),Z(NM,*)
  121. REAL RV1(*),RV2(*),RV3(*),RV4(*),RV5(*),RV6(*)
  122. REAL U,V,LB,T1,T2,UB,UK,XU,X0,X1,EPS1,EPS2,EPS3,EPS4
  123. REAL NORM,MACHEP,S1,S2
  124. LOGICAL FIRST
  125. C
  126. SAVE FIRST, MACHEP
  127. DATA FIRST /.TRUE./
  128. C***FIRST EXECUTABLE STATEMENT TSTURM
  129. IF (FIRST) THEN
  130. MACHEP = R1MACH(4)
  131. ENDIF
  132. FIRST = .FALSE.
  133. C
  134. IERR = 0
  135. T1 = LB
  136. T2 = UB
  137. C .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES ..........
  138. DO 40 I = 1, N
  139. IF (I .EQ. 1) GO TO 20
  140. S1 = ABS(D(I)) + ABS(D(I-1))
  141. S2 = S1 + ABS(E(I))
  142. IF (S2 .GT. S1) GO TO 40
  143. 20 E2(I) = 0.0E0
  144. 40 CONTINUE
  145. C .......... DETERMINE THE NUMBER OF EIGENVALUES
  146. C IN THE INTERVAL ..........
  147. P = 1
  148. Q = N
  149. X1 = UB
  150. ISTURM = 1
  151. GO TO 320
  152. 60 M = S
  153. X1 = LB
  154. ISTURM = 2
  155. GO TO 320
  156. 80 M = M - S
  157. IF (M .GT. MM) GO TO 980
  158. Q = 0
  159. R = 0
  160. C .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING
  161. C INTERVAL BY THE GERSCHGORIN BOUNDS ..........
  162. 100 IF (R .EQ. M) GO TO 1001
  163. P = Q + 1
  164. XU = D(P)
  165. X0 = D(P)
  166. U = 0.0E0
  167. C
  168. DO 120 Q = P, N
  169. X1 = U
  170. U = 0.0E0
  171. V = 0.0E0
  172. IF (Q .EQ. N) GO TO 110
  173. U = ABS(E(Q+1))
  174. V = E2(Q+1)
  175. 110 XU = MIN(D(Q)-(X1+U),XU)
  176. X0 = MAX(D(Q)+(X1+U),X0)
  177. IF (V .EQ. 0.0E0) GO TO 140
  178. 120 CONTINUE
  179. C
  180. 140 X1 = MAX(ABS(XU),ABS(X0)) * MACHEP
  181. IF (EPS1 .LE. 0.0E0) EPS1 = -X1
  182. IF (P .NE. Q) GO TO 180
  183. C .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL ..........
  184. IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940
  185. R = R + 1
  186. C
  187. DO 160 I = 1, N
  188. 160 Z(I,R) = 0.0E0
  189. C
  190. W(R) = D(P)
  191. Z(P,R) = 1.0E0
  192. GO TO 940
  193. 180 X1 = X1 * (Q-P+1)
  194. LB = MAX(T1,XU-X1)
  195. UB = MIN(T2,X0+X1)
  196. X1 = LB
  197. ISTURM = 3
  198. GO TO 320
  199. 200 M1 = S + 1
  200. X1 = UB
  201. ISTURM = 4
  202. GO TO 320
  203. 220 M2 = S
  204. IF (M1 .GT. M2) GO TO 940
  205. C .......... FIND ROOTS BY BISECTION ..........
  206. X0 = UB
  207. ISTURM = 5
  208. C
  209. DO 240 I = M1, M2
  210. RV5(I) = UB
  211. RV4(I) = LB
  212. 240 CONTINUE
  213. C .......... LOOP FOR K-TH EIGENVALUE
  214. C FOR K=M2 STEP -1 UNTIL M1 DO --
  215. C (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) ..........
  216. K = M2
  217. 250 XU = LB
  218. C .......... FOR I=K STEP -1 UNTIL M1 DO -- ..........
  219. DO 260 II = M1, K
  220. I = M1 + K - II
  221. IF (XU .GE. RV4(I)) GO TO 260
  222. XU = RV4(I)
  223. GO TO 280
  224. 260 CONTINUE
  225. C
  226. 280 IF (X0 .GT. RV5(K)) X0 = RV5(K)
  227. C .......... NEXT BISECTION STEP ..........
  228. 300 X1 = (XU + X0) * 0.5E0
  229. S1 = 2.0E0*(ABS(XU) + ABS(X0) + ABS(EPS1))
  230. S2 = S1 + ABS(X0 - XU)
  231. IF (S2 .EQ. S1) GO TO 420
  232. C .......... IN-LINE PROCEDURE FOR STURM SEQUENCE ..........
  233. 320 S = P - 1
  234. U = 1.0E0
  235. C
  236. DO 340 I = P, Q
  237. IF (U .NE. 0.0E0) GO TO 325
  238. V = ABS(E(I)) / MACHEP
  239. IF (E2(I) .EQ. 0.0E0) V = 0.0E0
  240. GO TO 330
  241. 325 V = E2(I) / U
  242. 330 U = D(I) - X1 - V
  243. IF (U .LT. 0.0E0) S = S + 1
  244. 340 CONTINUE
  245. C
  246. GO TO (60,80,200,220,360), ISTURM
  247. C .......... REFINE INTERVALS ..........
  248. 360 IF (S .GE. K) GO TO 400
  249. XU = X1
  250. IF (S .GE. M1) GO TO 380
  251. RV4(M1) = X1
  252. GO TO 300
  253. 380 RV4(S+1) = X1
  254. IF (RV5(S) .GT. X1) RV5(S) = X1
  255. GO TO 300
  256. 400 X0 = X1
  257. GO TO 300
  258. C .......... K-TH EIGENVALUE FOUND ..........
  259. 420 RV5(K) = X1
  260. K = K - 1
  261. IF (K .GE. M1) GO TO 250
  262. C .......... FIND VECTORS BY INVERSE ITERATION ..........
  263. NORM = ABS(D(P))
  264. IP = P + 1
  265. C
  266. DO 500 I = IP, Q
  267. 500 NORM = MAX(NORM, ABS(D(I)) + ABS(E(I)))
  268. C .......... EPS2 IS THE CRITERION FOR GROUPING,
  269. C EPS3 REPLACES ZERO PIVOTS AND EQUAL
  270. C ROOTS ARE MODIFIED BY EPS3,
  271. C EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW ..........
  272. EPS2 = 1.0E-3 * NORM
  273. UK = SQRT(REAL(Q-P+5))
  274. EPS3 = UK * MACHEP * NORM
  275. EPS4 = UK * EPS3
  276. UK = EPS4 / SQRT(UK)
  277. GROUP = 0
  278. S = P
  279. C
  280. DO 920 K = M1, M2
  281. R = R + 1
  282. ITS = 1
  283. W(R) = RV5(K)
  284. X1 = RV5(K)
  285. C .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
  286. IF (K .EQ. M1) GO TO 520
  287. IF (X1 - X0 .GE. EPS2) GROUP = -1
  288. GROUP = GROUP + 1
  289. IF (X1 .LE. X0) X1 = X0 + EPS3
  290. C .......... ELIMINATION WITH INTERCHANGES AND
  291. C INITIALIZATION OF VECTOR ..........
  292. 520 V = 0.0E0
  293. C
  294. DO 580 I = P, Q
  295. RV6(I) = UK
  296. IF (I .EQ. P) GO TO 560
  297. IF (ABS(E(I)) .LT. ABS(U)) GO TO 540
  298. XU = U / E(I)
  299. RV4(I) = XU
  300. RV1(I-1) = E(I)
  301. RV2(I-1) = D(I) - X1
  302. RV3(I-1) = 0.0E0
  303. IF (I .NE. Q) RV3(I-1) = E(I+1)
  304. U = V - XU * RV2(I-1)
  305. V = -XU * RV3(I-1)
  306. GO TO 580
  307. 540 XU = E(I) / U
  308. RV4(I) = XU
  309. RV1(I-1) = U
  310. RV2(I-1) = V
  311. RV3(I-1) = 0.0E0
  312. 560 U = D(I) - X1 - XU * V
  313. IF (I .NE. Q) V = E(I+1)
  314. 580 CONTINUE
  315. C
  316. IF (U .EQ. 0.0E0) U = EPS3
  317. RV1(Q) = U
  318. RV2(Q) = 0.0E0
  319. RV3(Q) = 0.0E0
  320. C .......... BACK SUBSTITUTION
  321. C FOR I=Q STEP -1 UNTIL P DO -- ..........
  322. 600 DO 620 II = P, Q
  323. I = P + Q - II
  324. RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
  325. V = U
  326. U = RV6(I)
  327. 620 CONTINUE
  328. C .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
  329. C MEMBERS OF GROUP ..........
  330. IF (GROUP .EQ. 0) GO TO 700
  331. C
  332. DO 680 JJ = 1, GROUP
  333. J = R - GROUP - 1 + JJ
  334. XU = 0.0E0
  335. C
  336. DO 640 I = P, Q
  337. 640 XU = XU + RV6(I) * Z(I,J)
  338. C
  339. DO 660 I = P, Q
  340. 660 RV6(I) = RV6(I) - XU * Z(I,J)
  341. C
  342. 680 CONTINUE
  343. C
  344. 700 NORM = 0.0E0
  345. C
  346. DO 720 I = P, Q
  347. 720 NORM = NORM + ABS(RV6(I))
  348. C
  349. IF (NORM .GE. 1.0E0) GO TO 840
  350. C .......... FORWARD SUBSTITUTION ..........
  351. IF (ITS .EQ. 5) GO TO 960
  352. IF (NORM .NE. 0.0E0) GO TO 740
  353. RV6(S) = EPS4
  354. S = S + 1
  355. IF (S .GT. Q) S = P
  356. GO TO 780
  357. 740 XU = EPS4 / NORM
  358. C
  359. DO 760 I = P, Q
  360. 760 RV6(I) = RV6(I) * XU
  361. C .......... ELIMINATION OPERATIONS ON NEXT VECTOR
  362. C ITERATE ..........
  363. 780 DO 820 I = IP, Q
  364. U = RV6(I)
  365. C .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
  366. C WAS PERFORMED EARLIER IN THE
  367. C TRIANGULARIZATION PROCESS ..........
  368. IF (RV1(I-1) .NE. E(I)) GO TO 800
  369. U = RV6(I-1)
  370. RV6(I-1) = RV6(I)
  371. 800 RV6(I) = U - RV4(I) * RV6(I-1)
  372. 820 CONTINUE
  373. C
  374. ITS = ITS + 1
  375. GO TO 600
  376. C .......... NORMALIZE SO THAT SUM OF SQUARES IS
  377. C 1 AND EXPAND TO FULL ORDER ..........
  378. 840 U = 0.0E0
  379. C
  380. DO 860 I = P, Q
  381. 860 U = U + RV6(I)**2
  382. C
  383. XU = 1.0E0 / SQRT(U)
  384. C
  385. DO 880 I = 1, N
  386. 880 Z(I,R) = 0.0E0
  387. C
  388. DO 900 I = P, Q
  389. 900 Z(I,R) = RV6(I) * XU
  390. C
  391. X0 = X1
  392. 920 CONTINUE
  393. C
  394. 940 IF (Q .LT. N) GO TO 100
  395. GO TO 1001
  396. C .......... SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
  397. 960 IERR = 4 * N + R
  398. GO TO 1001
  399. C .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF
  400. C EIGENVALUES IN INTERVAL ..........
  401. 980 IERR = 3 * N + 1
  402. 1001 LB = T1
  403. UB = T2
  404. RETURN
  405. END