Télécharger dsteqr.eso

Retour à la liste

Numérotation des lignes :

  1. C DSTEQR SOURCE BP208322 15/10/13 21:15:57 8670
  2. *> \brief \b DSTEQR
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DSTEQR + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsteqr.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsteqr.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsteqr.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
  23. *
  24. * .. Scalar Arguments ..
  25. * CHARACTER COMPZ
  26. * INTEGER INFO, LDZ, N
  27. * ..
  28. * .. Array Arguments ..
  29. * REAL*8 D( * ), E( * ), WORK( * ), Z( LDZ, * )
  30. * ..
  31. *
  32. *
  33. *> \par Purpose:
  34. * =============
  35. *>
  36. *> \verbatim
  37. *>
  38. *> DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
  39. *> symmetric tridiagonal matrix using the implicit QL or QR method.
  40. *> The eigenvectors of a full or band symmetric matrix can also be found
  41. *> if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
  42. *> tridiagonal form.
  43. *> \endverbatim
  44. *
  45. * Arguments:
  46. * ==========
  47. *
  48. *> \param[in] COMPZ
  49. *> \verbatim
  50. *> COMPZ is CHARACTER*1
  51. *> = 'N': Compute eigenvalues only.
  52. *> = 'V': Compute eigenvalues and eigenvectors of the original
  53. *> symmetric matrix. On entry, Z must contain the
  54. *> orthogonal matrix used to reduce the original matrix
  55. *> to tridiagonal form.
  56. *> = 'I': Compute eigenvalues and eigenvectors of the
  57. *> tridiagonal matrix. Z is initialized to the identity
  58. *> matrix.
  59. *> \endverbatim
  60. *>
  61. *> \param[in] N
  62. *> \verbatim
  63. *> N is INTEGER
  64. *> The order of the matrix. N >= 0.
  65. *> \endverbatim
  66. *>
  67. *> \param[in,out] D
  68. *> \verbatim
  69. *> D is DOUBLE PRECISION array, dimension (N)
  70. *> On entry, the diagonal elements of the tridiagonal matrix.
  71. *> On exit, if INFO = 0, the eigenvalues in ascending order.
  72. *> \endverbatim
  73. *>
  74. *> \param[in,out] E
  75. *> \verbatim
  76. *> E is DOUBLE PRECISION array, dimension (N-1)
  77. *> On entry, the (n-1) subdiagonal elements of the tridiagonal
  78. *> matrix.
  79. *> On exit, E has been destroyed.
  80. *> \endverbatim
  81. *>
  82. *> \param[in,out] Z
  83. *> \verbatim
  84. *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
  85. *> On entry, if COMPZ = 'V', then Z contains the orthogonal
  86. *> matrix used in the reduction to tridiagonal form.
  87. *> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
  88. *> orthonormal eigenvectors of the original symmetric matrix,
  89. *> and if COMPZ = 'I', Z contains the orthonormal eigenvectors
  90. *> of the symmetric tridiagonal matrix.
  91. *> If COMPZ = 'N', then Z is not referenced.
  92. *> \endverbatim
  93. *>
  94. *> \param[in] LDZ
  95. *> \verbatim
  96. *> LDZ is INTEGER
  97. *> The leading dimension of the array Z. LDZ >= 1, and if
  98. *> eigenvectors are desired, then LDZ >= max(1,N).
  99. *> \endverbatim
  100. *>
  101. *> \param[out] WORK
  102. *> \verbatim
  103. *> WORK is DOUBLE PRECISION array, dimension (max(1,2*N-2))
  104. *> If COMPZ = 'N', then WORK is not referenced.
  105. *> \endverbatim
  106. *>
  107. *> \param[out] INFO
  108. *> \verbatim
  109. *> INFO is INTEGER
  110. *> = 0: successful exit
  111. *> < 0: if INFO = -i, the i-th argument had an illegal value
  112. *> > 0: the algorithm has failed to find all the eigenvalues in
  113. *> a total of 30*N iterations; if INFO = i, then i
  114. *> elements of E have not converged to zero; on exit, D
  115. *> and E contain the elements of a symmetric tridiagonal
  116. *> matrix which is orthogonally similar to the original
  117. *> matrix.
  118. *> \endverbatim
  119. *
  120. * Authors:
  121. * ========
  122. *
  123. *> \author Univ. of Tennessee
  124. *> \author Univ. of California Berkeley
  125. *> \author Univ. of Colorado Denver
  126. *> \author NAG Ltd.
  127. *
  128. *> \date November 2011
  129. *
  130. *> \ingroup auxOTHERcomputational
  131. *
  132. * =====================================================================
  133. SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
  134. *
  135. * -- LAPACK computational routine (version 3.4.0) --
  136. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  137. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  138. * November 2011
  139. *
  140. * .. Scalar Arguments ..
  141. CHARACTER COMPZ
  142. INTEGER INFO, LDZ, N
  143. * ..
  144. * .. Array Arguments ..
  145. REAL*8 D( * ), E( * ), WORK( * ), Z( LDZ, * )
  146. * ..
  147. *
  148. * =====================================================================
  149. *
  150. * .. Parameters ..
  151. REAL*8 ZERO, ONE, TWO, THREE
  152. PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
  153. $ THREE = 3.0D0 )
  154. INTEGER MAXIT
  155. PARAMETER ( MAXIT = 30 )
  156. * ..
  157. * .. Local Scalars ..
  158. INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
  159. $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
  160. $ NM1, NMAXIT
  161. REAL*8 ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
  162. $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
  163. * ..
  164. * .. External Functions ..
  165. LOGICAL LSAME
  166. REAL*8 DLAMCH, DLANST, DLAPY2
  167. EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2
  168. * ..
  169. * .. External Subroutines ..
  170. * ..
  171. ** .. Intrinsic Functions ..
  172. * INTRINSIC ABS, MAX, SIGN, SQRT
  173. ** ..
  174. ** .. Executable Statements ..
  175. *
  176. * Test the input parameters.
  177. *
  178. INFO = 0
  179. *
  180. IF( LSAME( COMPZ, 'N' ) ) THEN
  181. ICOMPZ = 0
  182. ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
  183. ICOMPZ = 1
  184. ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
  185. ICOMPZ = 2
  186. ELSE
  187. ICOMPZ = -1
  188. END IF
  189. IF( ICOMPZ.LT.0 ) THEN
  190. INFO = -1
  191. ELSE IF( N.LT.0 ) THEN
  192. INFO = -2
  193. ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
  194. $ N ) ) ) THEN
  195. INFO = -6
  196. END IF
  197. IF( INFO.NE.0 ) THEN
  198. CALL XERBLA( 'DSTEQR', -INFO )
  199. RETURN
  200. END IF
  201. *
  202. * Quick return if possible
  203. *
  204. IF( N.EQ.0 )
  205. $ RETURN
  206. *
  207. IF( N.EQ.1 ) THEN
  208. IF( ICOMPZ.EQ.2 )
  209. $ Z( 1, 1 ) = ONE
  210. RETURN
  211. END IF
  212. *
  213. * Determine the unit roundoff and over/underflow thresholds.
  214. *
  215. EPS = DLAMCH( 'E' )
  216. EPS2 = EPS**2
  217. SAFMIN = DLAMCH( 'S' )
  218. SAFMAX = ONE / SAFMIN
  219. SSFMAX = SQRT( SAFMAX ) / THREE
  220. SSFMIN = SQRT( SAFMIN ) / EPS2
  221. *
  222. * Compute the eigenvalues and eigenvectors of the tridiagonal
  223. * matrix.
  224. *
  225. IF( ICOMPZ.EQ.2 )
  226. $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
  227. *
  228. NMAXIT = N*MAXIT
  229. JTOT = 0
  230. *
  231. * Determine where the matrix splits and choose QL or QR iteration
  232. * for each block, according to whether top or bottom diagonal
  233. * element is smaller.
  234. *
  235. L1 = 1
  236. NM1 = N - 1
  237. *
  238. 10 CONTINUE
  239. IF( L1.GT.N )
  240. $ GO TO 160
  241. IF( L1.GT.1 )
  242. $ E( L1-1 ) = ZERO
  243. IF( L1.LE.NM1 ) THEN
  244. DO 20 M = L1, NM1
  245. TST = ABS( E( M ) )
  246. IF( TST.EQ.ZERO )
  247. $ GO TO 30
  248. IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
  249. $ 1 ) ) ) )*EPS ) THEN
  250. E( M ) = ZERO
  251. GO TO 30
  252. END IF
  253. 20 CONTINUE
  254. END IF
  255. M = N
  256. *
  257. 30 CONTINUE
  258. L = L1
  259. LSV = L
  260. LEND = M
  261. LENDSV = LEND
  262. L1 = M + 1
  263. IF( LEND.EQ.L )
  264. $ GO TO 10
  265. *
  266. * Scale submatrix in rows and columns L to LEND
  267. *
  268. ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) )
  269. ISCALE = 0
  270. IF( ANORM.EQ.ZERO )
  271. $ GO TO 10
  272. IF( ANORM.GT.SSFMAX ) THEN
  273. ISCALE = 1
  274. CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
  275. $ INFO )
  276. CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
  277. $ INFO )
  278. ELSE IF( ANORM.LT.SSFMIN ) THEN
  279. ISCALE = 2
  280. CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
  281. $ INFO )
  282. CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
  283. $ INFO )
  284. END IF
  285. *
  286. * Choose between QL and QR iteration
  287. *
  288. IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
  289. LEND = LSV
  290. L = LENDSV
  291. END IF
  292. *
  293. IF( LEND.GT.L ) THEN
  294. *
  295. * QL Iteration
  296. *
  297. * Look for small subdiagonal element.
  298. *
  299. 40 CONTINUE
  300. IF( L.NE.LEND ) THEN
  301. LENDM1 = LEND - 1
  302. DO 50 M = L, LENDM1
  303. TST = ABS( E( M ) )**2
  304. IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
  305. $ SAFMIN )GO TO 60
  306. 50 CONTINUE
  307. END IF
  308. *
  309. M = LEND
  310. *
  311. 60 CONTINUE
  312. IF( M.LT.LEND )
  313. $ E( M ) = ZERO
  314. P = D( L )
  315. IF( M.EQ.L )
  316. $ GO TO 80
  317. *
  318. * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
  319. * to compute its eigensystem.
  320. *
  321. IF( M.EQ.L+1 ) THEN
  322. IF( ICOMPZ.GT.0 ) THEN
  323. CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
  324. WORK( L ) = C
  325. WORK( N-1+L ) = S
  326. CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ),
  327. $ WORK( N-1+L ), Z( 1, L ), LDZ )
  328. ELSE
  329. CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
  330. END IF
  331. D( L ) = RT1
  332. D( L+1 ) = RT2
  333. E( L ) = ZERO
  334. L = L + 2
  335. IF( L.LE.LEND )
  336. $ GO TO 40
  337. GO TO 140
  338. END IF
  339. *
  340. IF( JTOT.EQ.NMAXIT )
  341. $ GO TO 140
  342. JTOT = JTOT + 1
  343. *
  344. * Form shift.
  345. *
  346. G = ( D( L+1 )-P ) / ( TWO*E( L ) )
  347. R = DLAPY2( G, ONE )
  348. G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
  349. *
  350. S = ONE
  351. C = ONE
  352. P = ZERO
  353. *
  354. * Inner loop
  355. *
  356. MM1 = M - 1
  357. DO 70 I = MM1, L, -1
  358. F = S*E( I )
  359. B = C*E( I )
  360. CALL DLARTG( G, F, C, S, R )
  361. IF( I.NE.M-1 )
  362. $ E( I+1 ) = R
  363. G = D( I+1 ) - P
  364. R = ( D( I )-G )*S + TWO*C*B
  365. P = S*R
  366. D( I+1 ) = G + P
  367. G = C*R - B
  368. *
  369. * If eigenvectors are desired, then save rotations.
  370. *
  371. IF( ICOMPZ.GT.0 ) THEN
  372. WORK( I ) = C
  373. WORK( N-1+I ) = -S
  374. END IF
  375. *
  376. 70 CONTINUE
  377. *
  378. * If eigenvectors are desired, then apply saved rotations.
  379. *
  380. IF( ICOMPZ.GT.0 ) THEN
  381. MM = M - L + 1
  382. CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
  383. $ Z( 1, L ), LDZ )
  384. END IF
  385. *
  386. D( L ) = D( L ) - P
  387. E( L ) = G
  388. GO TO 40
  389. *
  390. * Eigenvalue found.
  391. *
  392. 80 CONTINUE
  393. D( L ) = P
  394. *
  395. L = L + 1
  396. IF( L.LE.LEND )
  397. $ GO TO 40
  398. GO TO 140
  399. *
  400. ELSE
  401. *
  402. * QR Iteration
  403. *
  404. * Look for small superdiagonal element.
  405. *
  406. 90 CONTINUE
  407. IF( L.NE.LEND ) THEN
  408. LENDP1 = LEND + 1
  409. DO 100 M = L, LENDP1, -1
  410. TST = ABS( E( M-1 ) )**2
  411. IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
  412. $ SAFMIN )GO TO 110
  413. 100 CONTINUE
  414. END IF
  415. *
  416. M = LEND
  417. *
  418. 110 CONTINUE
  419. IF( M.GT.LEND )
  420. $ E( M-1 ) = ZERO
  421. P = D( L )
  422. IF( M.EQ.L )
  423. $ GO TO 130
  424. *
  425. * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
  426. * to compute its eigensystem.
  427. *
  428. IF( M.EQ.L-1 ) THEN
  429. IF( ICOMPZ.GT.0 ) THEN
  430. CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
  431. WORK( M ) = C
  432. WORK( N-1+M ) = S
  433. CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ),
  434. $ WORK( N-1+M ), Z( 1, L-1 ), LDZ )
  435. ELSE
  436. CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
  437. END IF
  438. D( L-1 ) = RT1
  439. D( L ) = RT2
  440. E( L-1 ) = ZERO
  441. L = L - 2
  442. IF( L.GE.LEND )
  443. $ GO TO 90
  444. GO TO 140
  445. END IF
  446. *
  447. IF( JTOT.EQ.NMAXIT )
  448. $ GO TO 140
  449. JTOT = JTOT + 1
  450. *
  451. * Form shift.
  452. *
  453. G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
  454. R = DLAPY2( G, ONE )
  455. G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
  456. *
  457. S = ONE
  458. C = ONE
  459. P = ZERO
  460. *
  461. * Inner loop
  462. *
  463. LM1 = L - 1
  464. DO 120 I = M, LM1
  465. F = S*E( I )
  466. B = C*E( I )
  467. CALL DLARTG( G, F, C, S, R )
  468. IF( I.NE.M )
  469. $ E( I-1 ) = R
  470. G = D( I ) - P
  471. R = ( D( I+1 )-G )*S + TWO*C*B
  472. P = S*R
  473. D( I ) = G + P
  474. G = C*R - B
  475. *
  476. * If eigenvectors are desired, then save rotations.
  477. *
  478. IF( ICOMPZ.GT.0 ) THEN
  479. WORK( I ) = C
  480. WORK( N-1+I ) = S
  481. END IF
  482. *
  483. 120 CONTINUE
  484. *
  485. * If eigenvectors are desired, then apply saved rotations.
  486. *
  487. IF( ICOMPZ.GT.0 ) THEN
  488. MM = L - M + 1
  489. CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
  490. $ Z( 1, M ), LDZ )
  491. END IF
  492. *
  493. D( L ) = D( L ) - P
  494. E( LM1 ) = G
  495. GO TO 90
  496. *
  497. * Eigenvalue found.
  498. *
  499. 130 CONTINUE
  500. D( L ) = P
  501. *
  502. L = L - 1
  503. IF( L.GE.LEND )
  504. $ GO TO 90
  505. GO TO 140
  506. *
  507. END IF
  508. *
  509. * Undo scaling if necessary
  510. *
  511. 140 CONTINUE
  512. IF( ISCALE.EQ.1 ) THEN
  513. CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
  514. $ D( LSV ), N, INFO )
  515. CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
  516. $ N, INFO )
  517. ELSE IF( ISCALE.EQ.2 ) THEN
  518. CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
  519. $ D( LSV ), N, INFO )
  520. CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
  521. $ N, INFO )
  522. END IF
  523. *
  524. * Check for no convergence to an eigenvalue after a total
  525. * of N*MAXIT iterations.
  526. *
  527. IF( JTOT.LT.NMAXIT )
  528. $ GO TO 10
  529. DO 150 I = 1, N - 1
  530. IF( E( I ).NE.ZERO )
  531. $ INFO = INFO + 1
  532. 150 CONTINUE
  533. GO TO 190
  534. *
  535. * Order eigenvalues and eigenvectors.
  536. *
  537. 160 CONTINUE
  538. IF( ICOMPZ.EQ.0 ) THEN
  539. *
  540. * Use Quick Sort
  541. *
  542. CALL DLASRT( 'I', N, D, INFO )
  543. *
  544. ELSE
  545. *
  546. * Use Selection Sort to minimize swaps of eigenvectors
  547. *
  548. DO 180 II = 2, N
  549. I = II - 1
  550. K = I
  551. P = D( I )
  552. DO 170 J = II, N
  553. IF( D( J ).LT.P ) THEN
  554. K = J
  555. P = D( J )
  556. END IF
  557. 170 CONTINUE
  558. IF( K.NE.I ) THEN
  559. D( K ) = D( I )
  560. D( I ) = P
  561. CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
  562. END IF
  563. 180 CONTINUE
  564. END IF
  565. *
  566. 190 CONTINUE
  567. RETURN
  568. *
  569. * End of DSTEQR
  570. *
  571. END
  572.  
  573.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales