Télécharger dstqrb.eso

Retour à la liste

Numérotation des lignes :

  1. C DSTQRB SOURCE BP208322 15/10/13 21:15:57 8670
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dstqrb
  6. c
  7. c\Description:
  8. c Computes all eigenvalues and the last component of the eigenvectors
  9. c of a symmetric tridiagonal matrix using the implicit QL or QR method.
  10. c
  11. c This is mostly a modification of the LAPACK routine dsteqr.
  12. c See Remarks.
  13. c
  14. c\Usage:
  15. c call dstqrb
  16. c ( N, D, E, Z, WORK, INFO )
  17. c
  18. c\Arguments
  19. c N Integer. (INPUT)
  20. c The number of rows and columns in the matrix. N >= 0.
  21. c
  22. c D REAL*8 array, dimension (N). (INPUT/OUTPUT)
  23. c On entry, D contains the diagonal elements of the
  24. c tridiagonal matrix.
  25. c On exit, D contains the eigenvalues, in ascending order.
  26. c If an error exit is made, the eigenvalues are correct
  27. c for indices 1,2,...,INFO-1, but they are unordered and
  28. c may not be the smallest eigenvalues of the matrix.
  29. c
  30. c E REAL*8 array, dimension (N-1). (INPUT/OUTPUT)
  31. c On entry, E contains the subdiagonal elements of the
  32. c tridiagonal matrix in positions 1 through N-1.
  33. c On exit, E has been destroyed.
  34. c
  35. c Z REAL*8 array, dimension (N). (OUTPUT)
  36. c On exit, Z contains the last row of the orthonormal
  37. c eigenvector matrix of the symmetric tridiagonal matrix.
  38. c If an error exit is made, Z contains the last row of the
  39. c eigenvector matrix associated with the stored eigenvalues.
  40. c
  41. c WORK Double precision array, dimension (max(1,2*N-2)). (WORKSPACE)
  42. c Workspace used in accumulating the transformation for
  43. c computing the last components of the eigenvectors.
  44. c
  45. c INFO Integer. (OUTPUT)
  46. c = 0: normal return.
  47. c < 0: if INFO = -i, the i-th argument had an illegal value.
  48. c > 0: if INFO = +i, the i-th eigenvalue has not converged
  49. c after a total of 30*N iterations.
  50. c
  51. c\Remarks
  52. c 1. None.
  53. c
  54. c-----------------------------------------------------------------------
  55. c
  56. c\BeginLib
  57. c
  58. c\Local variables:
  59. c xxxxxx real
  60. c
  61. c\Routines called:
  62. c daxpy Level 1 BLAS that computes a vector triad.
  63. c dcopy Level 1 BLAS that copies one vector to another.
  64. c dswap Level 1 BLAS that swaps the contents of two vectors.
  65. c lsame LAPACK character comparison routine.
  66. c dlae2 LAPACK routine that computes the eigenvalues of a 2-by-2
  67. c symmetric matrix.
  68. c dlaev2 LAPACK routine that eigendecomposition of a 2-by-2 symmetric
  69. c matrix.
  70. c dlamch LAPACK routine that determines machine constants.
  71. c dlanst LAPACK routine that computes the norm of a matrix.
  72. c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
  73. c dlartg LAPACK Givens rotation construction routine.
  74. c dlascl LAPACK routine for careful scaling of a matrix.
  75. c dlaset LAPACK matrix initialization routine.
  76. c dlasr LAPACK routine that applies an orthogonal transformation to
  77. c a matrix.
  78. c dlasrt LAPACK sorting routine.
  79. c dsteqr LAPACK routine that computes eigenvalues and eigenvectors
  80. c of a symmetric tridiagonal matrix.
  81. c xerbla LAPACK error handler routine.
  82. c
  83. c\Authors
  84. c Danny Sorensen Phuong Vu
  85. c Richard Lehoucq CRPC / Rice University
  86. c Dept. of Computational & Houston, Texas
  87. c Applied Mathematics
  88. c Rice University
  89. c Houston, Texas
  90. c
  91. c\SCCS Information: @(#)
  92. c FILE: stqrb.F SID: 2.5 DATE OF SID: 8/27/96 RELEASE: 2
  93. c
  94. c\Remarks
  95. c 1. Starting with version 2.5, this routine is a modified version
  96. c of LAPACK version 2.0 subroutine SSTEQR. No lines are deleted,
  97. c only commeted out and new lines inserted.
  98. c All lines commented out have "c$$$" at the beginning.
  99. c Note that the LAPACK version 1.0 subroutine SSTEQR contained
  100. c bugs.
  101. c
  102. c\EndLib
  103. c
  104. c-----------------------------------------------------------------------
  105. c
  106. subroutine dstqrb ( n, d, e, z, work, info )
  107. c
  108. c %------------------%
  109. c | Scalar Arguments |
  110. c %------------------%
  111. c
  112. integer info, n
  113. c
  114. c %-----------------%
  115. c | Array Arguments |
  116. c %-----------------%
  117. c
  118. REAL*8
  119. & d( n ), e( n-1 ), z( n ), work( 2*n-2 )
  120. c
  121. c .. parameters ..
  122. REAL*8
  123. & zero, one, two, three
  124. parameter ( zero = 0.0D+0, one = 1.0D+0,
  125. & two = 2.0D+0, three = 3.0D+0 )
  126. integer maxit
  127. parameter ( maxit = 30 )
  128. c ..
  129. c .. local scalars ..
  130. integer i, icompz, ii, iscale, j, jtot, k, l, l1, lend,
  131. & lendm1, lendp1, lendsv, lm1, lsv, m, mm, mm1,
  132. & nm1, nmaxit
  133. REAL*8
  134. & anorm, b, c, eps, eps2, f, g, p, r, rt1, rt2,
  135. & s, safmax, safmin, ssfmax, ssfmin, tst
  136. c ..
  137. c .. external functions ..
  138. logical lsame
  139. REAL*8
  140. external lsame, dlamch, dlanst, dlapy2
  141. c ..
  142. c .. external subroutines ..
  143. c ..
  144. *c .. intrinsic functions ..
  145. * intrinsic abs, max, sign, sqrt
  146. *c ..
  147. *c .. executable statements ..
  148. c
  149. c test the input parameters.
  150. c
  151. info = 0
  152. c
  153. c$$$ IF( LSAME( COMPZ, 'N' ) ) THEN
  154. c$$$ ICOMPZ = 0
  155. c$$$ ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
  156. c$$$ ICOMPZ = 1
  157. c$$$ ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
  158. c$$$ ICOMPZ = 2
  159. c$$$ ELSE
  160. c$$$ ICOMPZ = -1
  161. c$$$ END IF
  162. c$$$ IF( ICOMPZ.LT.0 ) THEN
  163. c$$$ INFO = -1
  164. c$$$ ELSE IF( N.LT.0 ) THEN
  165. c$$$ INFO = -2
  166. c$$$ ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
  167. c$$$ $ N ) ) ) THEN
  168. c$$$ INFO = -6
  169. c$$$ END IF
  170. c$$$ IF( INFO.NE.0 ) THEN
  171. c$$$ CALL XERBLA( 'SSTEQR', -INFO )
  172. c$$$ RETURN
  173. c$$$ END IF
  174. c
  175. c *** New starting with version 2.5 ***
  176. c
  177. icompz = 2
  178. c *************************************
  179. c
  180. c quick return if possible
  181. c
  182. if( n.eq.0 )
  183. $ return
  184. c
  185. if( n.eq.1 ) then
  186. if( icompz.eq.2 ) z( 1 ) = one
  187. return
  188. end if
  189. c
  190. c determine the unit roundoff and over/underflow thresholds.
  191. c
  192. eps = dlamch( 'e' )
  193. eps2 = eps**2
  194. safmin = dlamch( 's' )
  195. safmax = one / safmin
  196. ssfmax = sqrt( safmax ) / three
  197. ssfmin = sqrt( safmin ) / eps2
  198. c
  199. c compute the eigenvalues and eigenvectors of the tridiagonal
  200. c matrix.
  201. c
  202. c$$ if( icompz.eq.2 )
  203. c$$$ $ call dlaset( 'full', n, n, zero, one, z, ldz )
  204. c
  205. c *** New starting with version 2.5 ***
  206. c
  207. if ( icompz .eq. 2 ) then
  208. do 5 j = 1, n-1
  209. z(j) = zero
  210. 5 continue
  211. z( n ) = one
  212. end if
  213. c *************************************
  214. c
  215. nmaxit = n*maxit
  216. jtot = 0
  217. c
  218. c determine where the matrix splits and choose ql or qr iteration
  219. c for each block, according to whether top or bottom diagonal
  220. c element is smaller.
  221. c
  222. l1 = 1
  223. nm1 = n - 1
  224. c
  225. 10 continue
  226. if( l1.gt.n )
  227. $ go to 160
  228. if( l1.gt.1 )
  229. $ e( l1-1 ) = zero
  230. if( l1.le.nm1 ) then
  231. do 20 m = l1, nm1
  232. tst = abs( e( m ) )
  233. if( tst.eq.zero )
  234. $ go to 30
  235. if( tst.le.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
  236. $ 1 ) ) ) )*eps ) then
  237. e( m ) = zero
  238. go to 30
  239. end if
  240. 20 continue
  241. end if
  242. m = n
  243. c
  244. 30 continue
  245. l = l1
  246. lsv = l
  247. lend = m
  248. lendsv = lend
  249. l1 = m + 1
  250. if( lend.eq.l )
  251. $ go to 10
  252. c
  253. c scale submatrix in rows and columns l to lend
  254. c
  255. anorm = dlanst( 'i', lend-l+1, d( l ), e( l ) )
  256. iscale = 0
  257. if( anorm.eq.zero )
  258. $ go to 10
  259. if( anorm.gt.ssfmax ) then
  260. iscale = 1
  261. call dlascl( 'g', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
  262. $ info )
  263. call dlascl( 'g', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
  264. $ info )
  265. else if( anorm.lt.ssfmin ) then
  266. iscale = 2
  267. call dlascl( 'g', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
  268. $ info )
  269. call dlascl( 'g', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
  270. $ info )
  271. end if
  272. c
  273. c choose between ql and qr iteration
  274. c
  275. if( abs( d( lend ) ).lt.abs( d( l ) ) ) then
  276. lend = lsv
  277. l = lendsv
  278. end if
  279. c
  280. if( lend.gt.l ) then
  281. c
  282. c ql iteration
  283. c
  284. c look for small subdiagonal element.
  285. c
  286. 40 continue
  287. if( l.ne.lend ) then
  288. lendm1 = lend - 1
  289. do 50 m = l, lendm1
  290. tst = abs( e( m ) )**2
  291. if( tst.le.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
  292. $ safmin )go to 60
  293. 50 continue
  294. end if
  295. c
  296. m = lend
  297. c
  298. 60 continue
  299. if( m.lt.lend )
  300. $ e( m ) = zero
  301. p = d( l )
  302. if( m.eq.l )
  303. $ go to 80
  304. c
  305. c if remaining matrix is 2-by-2, use dlae2 or dlaev2
  306. c to compute its eigensystem.
  307. c
  308. if( m.eq.l+1 ) then
  309. if( icompz.gt.0 ) then
  310. call dlaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
  311. work( l ) = c
  312. work( n-1+l ) = s
  313. c$$$ call dlasr( 'r', 'v', 'b', n, 2, work( l ),
  314. c$$$ $ work( n-1+l ), z( 1, l ), ldz )
  315. c
  316. c *** New starting with version 2.5 ***
  317. c
  318. tst = z(l+1)
  319. z(l+1) = c*tst - s*z(l)
  320. z(l) = s*tst + c*z(l)
  321. c *************************************
  322. else
  323. call dlae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
  324. end if
  325. d( l ) = rt1
  326. d( l+1 ) = rt2
  327. e( l ) = zero
  328. l = l + 2
  329. if( l.le.lend )
  330. $ go to 40
  331. go to 140
  332. end if
  333. c
  334. if( jtot.eq.nmaxit )
  335. $ go to 140
  336. jtot = jtot + 1
  337. c
  338. c form shift.
  339. c
  340. g = ( d( l+1 )-p ) / ( two*e( l ) )
  341. r = dlapy2( g, one )
  342. g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
  343. c
  344. s = one
  345. c = one
  346. p = zero
  347. c
  348. c inner loop
  349. c
  350. mm1 = m - 1
  351. do 70 i = mm1, l, -1
  352. f = s*e( i )
  353. b = c*e( i )
  354. call dlartg( g, f, c, s, r )
  355. if( i.ne.m-1 )
  356. $ e( i+1 ) = r
  357. g = d( i+1 ) - p
  358. r = ( d( i )-g )*s + two*c*b
  359. p = s*r
  360. d( i+1 ) = g + p
  361. g = c*r - b
  362. c
  363. c if eigenvectors are desired, then save rotations.
  364. c
  365. if( icompz.gt.0 ) then
  366. work( i ) = c
  367. work( n-1+i ) = -s
  368. end if
  369. c
  370. 70 continue
  371. c
  372. c if eigenvectors are desired, then apply saved rotations.
  373. c
  374. if( icompz.gt.0 ) then
  375. mm = m - l + 1
  376. c$$$ call dlasr( 'r', 'v', 'b', n, mm, work( l ), work( n-1+l ),
  377. c$$$ $ z( 1, l ), ldz )
  378. c
  379. c *** New starting with version 2.5 ***
  380. c
  381. call dlasr( 'r', 'v', 'b', 1, mm, work( l ),
  382. & work( n-1+l ), z( l ), 1 )
  383. c *************************************
  384. end if
  385. c
  386. d( l ) = d( l ) - p
  387. e( l ) = g
  388. go to 40
  389. c
  390. c eigenvalue found.
  391. c
  392. 80 continue
  393. d( l ) = p
  394. c
  395. l = l + 1
  396. if( l.le.lend )
  397. $ go to 40
  398. go to 140
  399. c
  400. else
  401. c
  402. c qr iteration
  403. c
  404. c look for small superdiagonal element.
  405. c
  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. c
  416. m = lend
  417. c
  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. c
  425. c if remaining matrix is 2-by-2, use dlae2 or dlaev2
  426. c to compute its eigensystem.
  427. c
  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. c$$$ work( m ) = c
  432. c$$$ work( n-1+m ) = s
  433. c$$$ call dlasr( 'r', 'v', 'f', n, 2, work( m ),
  434. c$$$ $ work( n-1+m ), z( 1, l-1 ), ldz )
  435. c
  436. c *** New starting with version 2.5 ***
  437. c
  438. tst = z(l)
  439. z(l) = c*tst - s*z(l-1)
  440. z(l-1) = s*tst + c*z(l-1)
  441. c *************************************
  442. else
  443. call dlae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
  444. end if
  445. d( l-1 ) = rt1
  446. d( l ) = rt2
  447. e( l-1 ) = zero
  448. l = l - 2
  449. if( l.ge.lend )
  450. $ go to 90
  451. go to 140
  452. end if
  453. c
  454. if( jtot.eq.nmaxit )
  455. $ go to 140
  456. jtot = jtot + 1
  457. c
  458. c form shift.
  459. c
  460. g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
  461. r = dlapy2( g, one )
  462. g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
  463. c
  464. s = one
  465. c = one
  466. p = zero
  467. c
  468. c inner loop
  469. c
  470. lm1 = l - 1
  471. do 120 i = m, lm1
  472. f = s*e( i )
  473. b = c*e( i )
  474. call dlartg( g, f, c, s, r )
  475. if( i.ne.m )
  476. $ e( i-1 ) = r
  477. g = d( i ) - p
  478. r = ( d( i+1 )-g )*s + two*c*b
  479. p = s*r
  480. d( i ) = g + p
  481. g = c*r - b
  482. c
  483. c if eigenvectors are desired, then save rotations.
  484. c
  485. if( icompz.gt.0 ) then
  486. work( i ) = c
  487. work( n-1+i ) = s
  488. end if
  489. c
  490. 120 continue
  491. c
  492. c if eigenvectors are desired, then apply saved rotations.
  493. c
  494. if( icompz.gt.0 ) then
  495. mm = l - m + 1
  496. c$$$ call dlasr( 'r', 'v', 'f', n, mm, work( m ), work( n-1+m ),
  497. c$$$ $ z( 1, m ), ldz )
  498. c
  499. c *** New starting with version 2.5 ***
  500. c
  501. call dlasr( 'r', 'v', 'f', 1, mm, work( m ), work( n-1+m ),
  502. & z( m ), 1 )
  503. c *************************************
  504. end if
  505. c
  506. d( l ) = d( l ) - p
  507. e( lm1 ) = g
  508. go to 90
  509. c
  510. c eigenvalue found.
  511. c
  512. 130 continue
  513. d( l ) = p
  514. c
  515. l = l - 1
  516. if( l.ge.lend )
  517. $ go to 90
  518. go to 140
  519. c
  520. end if
  521. c
  522. c undo scaling if necessary
  523. c
  524. 140 continue
  525. if( iscale.eq.1 ) then
  526. call dlascl( 'g', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
  527. $ d( lsv ), n, info )
  528. call dlascl( 'g', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
  529. $ n, info )
  530. else if( iscale.eq.2 ) then
  531. call dlascl( 'g', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
  532. $ d( lsv ), n, info )
  533. call dlascl( 'g', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
  534. $ n, info )
  535. end if
  536. c
  537. c check for no convergence to an eigenvalue after a total
  538. c of n*maxit iterations.
  539. c
  540. if( jtot.lt.nmaxit )
  541. $ go to 10
  542. do 150 i = 1, n - 1
  543. if( e( i ).ne.zero )
  544. $ info = info + 1
  545. 150 continue
  546. go to 190
  547. c
  548. c order eigenvalues and eigenvectors.
  549. c
  550. 160 continue
  551. if( icompz.eq.0 ) then
  552. c
  553. c use quick sort
  554. c
  555. call dlasrt( 'i', n, d, info )
  556. c
  557. else
  558. c
  559. c use selection sort to minimize swaps of eigenvectors
  560. c
  561. do 180 ii = 2, n
  562. i = ii - 1
  563. k = i
  564. p = d( i )
  565. do 170 j = ii, n
  566. if( d( j ).lt.p ) then
  567. k = j
  568. p = d( j )
  569. end if
  570. 170 continue
  571. if( k.ne.i ) then
  572. d( k ) = d( i )
  573. d( i ) = p
  574. c$$$ call dswap( n, z( 1, i ), 1, z( 1, k ), 1 )
  575. c *** New starting with version 2.5 ***
  576. c
  577. p = z(k)
  578. z(k) = z(i)
  579. z(i) = p
  580. c *************************************
  581. end if
  582. 180 continue
  583. end if
  584. c
  585. 190 continue
  586. return
  587. c
  588. c %---------------%
  589. c | End of dstqrb |
  590. c %---------------%
  591. c
  592. end
  593.  
  594.  

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