Télécharger dtrevc.eso

Retour à la liste

Numérotation des lignes :

  1. C DTREVC SOURCE BP208322 15/10/13 21:15:59 8670
  2. *> \brief \b DTREVC
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DTREVC + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  23. * LDVR, MM, M, WORK, INFO )
  24. *
  25. * .. Scalar Arguments ..
  26. * CHARACTER HOWMNY, SIDE
  27. * INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
  28. * ..
  29. * .. Array Arguments ..
  30. * LOGICAL SELECT( * )
  31. * REAL*8 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  32. * $ WORK( * )
  33. * ..
  34. *
  35. *
  36. *> \par Purpose:
  37. * =============
  38. *>
  39. *> \verbatim
  40. *>
  41. *> DTREVC computes some or all of the right and/or left eigenvectors of
  42. *> a real upper quasi-triangular matrix T.
  43. *> Matrices of this type are produced by the Schur factorization of
  44. *> a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
  45. *>
  46. *> The right eigenvector x and the left eigenvector y of T corresponding
  47. *> to an eigenvalue w are defined by:
  48. *>
  49. *> T*x = w*x, (y**T)*T = w*(y**T)
  50. *>
  51. *> where y**T denotes the transpose of y.
  52. *> The eigenvalues are not input to this routine, but are read directly
  53. *> from the diagonal blocks of T.
  54. *>
  55. *> This routine returns the matrices X and/or Y of right and left
  56. *> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
  57. *> input matrix. If Q is the orthogonal factor that reduces a matrix
  58. *> A to Schur form T, then Q*X and Q*Y are the matrices of right and
  59. *> left eigenvectors of A.
  60. *> \endverbatim
  61. *
  62. * Arguments:
  63. * ==========
  64. *
  65. *> \param[in] SIDE
  66. *> \verbatim
  67. *> SIDE is CHARACTER*1
  68. *> = 'R': compute right eigenvectors only;
  69. *> = 'L': compute left eigenvectors only;
  70. *> = 'B': compute both right and left eigenvectors.
  71. *> \endverbatim
  72. *>
  73. *> \param[in] HOWMNY
  74. *> \verbatim
  75. *> HOWMNY is CHARACTER*1
  76. *> = 'A': compute all right and/or left eigenvectors;
  77. *> = 'B': compute all right and/or left eigenvectors,
  78. *> backtransformed by the matrices in VR and/or VL;
  79. *> = 'S': compute selected right and/or left eigenvectors,
  80. *> as indicated by the logical array SELECT.
  81. *> \endverbatim
  82. *>
  83. *> \param[in,out] SELECT
  84. *> \verbatim
  85. *> SELECT is LOGICAL array, dimension (N)
  86. *> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
  87. *> computed.
  88. *> If w(j) is a real eigenvalue, the corresponding real
  89. *> eigenvector is computed if SELECT(j) is .TRUE..
  90. *> If w(j) and w(j+1) are the real and imaginary parts of a
  91. *> complex eigenvalue, the corresponding complex eigenvector is
  92. *> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
  93. *> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
  94. *> .FALSE..
  95. *> Not referenced if HOWMNY = 'A' or 'B'.
  96. *> \endverbatim
  97. *>
  98. *> \param[in] N
  99. *> \verbatim
  100. *> N is INTEGER
  101. *> The order of the matrix T. N >= 0.
  102. *> \endverbatim
  103. *>
  104. *> \param[in] T
  105. *> \verbatim
  106. *> T is REAL*8 array, dimension (LDT,N)
  107. *> The upper quasi-triangular matrix T in Schur canonical form.
  108. *> \endverbatim
  109. *>
  110. *> \param[in] LDT
  111. *> \verbatim
  112. *> LDT is INTEGER
  113. *> The leading dimension of the array T. LDT >= max(1,N).
  114. *> \endverbatim
  115. *>
  116. *> \param[in,out] VL
  117. *> \verbatim
  118. *> VL is REAL*8 array, dimension (LDVL,MM)
  119. *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
  120. *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
  121. *> of Schur vectors returned by DHSEQR).
  122. *> On exit, if SIDE = 'L' or 'B', VL contains:
  123. *> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
  124. *> if HOWMNY = 'B', the matrix Q*Y;
  125. *> if HOWMNY = 'S', the left eigenvectors of T specified by
  126. *> SELECT, stored consecutively in the columns
  127. *> of VL, in the same order as their
  128. *> eigenvalues.
  129. *> A complex eigenvector corresponding to a complex eigenvalue
  130. *> is stored in two consecutive columns, the first holding the
  131. *> real part, and the second the imaginary part.
  132. *> Not referenced if SIDE = 'R'.
  133. *> \endverbatim
  134. *>
  135. *> \param[in] LDVL
  136. *> \verbatim
  137. *> LDVL is INTEGER
  138. *> The leading dimension of the array VL. LDVL >= 1, and if
  139. *> SIDE = 'L' or 'B', LDVL >= N.
  140. *> \endverbatim
  141. *>
  142. *> \param[in,out] VR
  143. *> \verbatim
  144. *> VR is REAL*8 array, dimension (LDVR,MM)
  145. *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
  146. *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
  147. *> of Schur vectors returned by DHSEQR).
  148. *> On exit, if SIDE = 'R' or 'B', VR contains:
  149. *> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
  150. *> if HOWMNY = 'B', the matrix Q*X;
  151. *> if HOWMNY = 'S', the right eigenvectors of T specified by
  152. *> SELECT, stored consecutively in the columns
  153. *> of VR, in the same order as their
  154. *> eigenvalues.
  155. *> A complex eigenvector corresponding to a complex eigenvalue
  156. *> is stored in two consecutive columns, the first holding the
  157. *> real part and the second the imaginary part.
  158. *> Not referenced if SIDE = 'L'.
  159. *> \endverbatim
  160. *>
  161. *> \param[in] LDVR
  162. *> \verbatim
  163. *> LDVR is INTEGER
  164. *> The leading dimension of the array VR. LDVR >= 1, and if
  165. *> SIDE = 'R' or 'B', LDVR >= N.
  166. *> \endverbatim
  167. *>
  168. *> \param[in] MM
  169. *> \verbatim
  170. *> MM is INTEGER
  171. *> The number of columns in the arrays VL and/or VR. MM >= M.
  172. *> \endverbatim
  173. *>
  174. *> \param[out] M
  175. *> \verbatim
  176. *> M is INTEGER
  177. *> The number of columns in the arrays VL and/or VR actually
  178. *> used to store the eigenvectors.
  179. *> If HOWMNY = 'A' or 'B', M is set to N.
  180. *> Each selected real eigenvector occupies one column and each
  181. *> selected complex eigenvector occupies two columns.
  182. *> \endverbatim
  183. *>
  184. *> \param[out] WORK
  185. *> \verbatim
  186. *> WORK is REAL*8 array, dimension (3*N)
  187. *> \endverbatim
  188. *>
  189. *> \param[out] INFO
  190. *> \verbatim
  191. *> INFO is INTEGER
  192. *> = 0: successful exit
  193. *> < 0: if INFO = -i, the i-th argument had an illegal value
  194. *> \endverbatim
  195. *
  196. * Authors:
  197. * ========
  198. *
  199. *> \author Univ. of Tennessee
  200. *> \author Univ. of California Berkeley
  201. *> \author Univ. of Colorado Denver
  202. *> \author NAG Ltd.
  203. *
  204. *> \date November 2011
  205. *
  206. *> \ingroup doubleOTHERcomputational
  207. *
  208. *> \par Further Details:
  209. * =====================
  210. *>
  211. *> \verbatim
  212. *>
  213. *> The algorithm used in this program is basically backward (forward)
  214. *> substitution, with scaling to make the the code robust against
  215. *> possible overflow.
  216. *>
  217. *> Each eigenvector is normalized so that the element of largest
  218. *> magnitude has magnitude 1; here the magnitude of a complex number
  219. *> (x,y) is taken to be |x| + |y|.
  220. *> \endverbatim
  221. *>
  222. * =====================================================================
  223. SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  224. $ LDVR, MM, M, WORK, INFO )
  225. *
  226. * -- LAPACK computational routine (version 3.4.0) --
  227. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  228. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  229. * November 2011
  230. *
  231. * .. Scalar Arguments ..
  232. CHARACTER HOWMNY, SIDE
  233. INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
  234. * ..
  235. * .. Array Arguments ..
  236. LOGICAL SELECT( * )
  237. REAL*8 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  238. $ WORK( * )
  239. * ..
  240. *
  241. * =====================================================================
  242. *
  243. * .. Parameters ..
  244. REAL*8 ZERO, ONE
  245. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  246. * ..
  247. * .. Local Scalars ..
  248. LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
  249. INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
  250. REAL*8 BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
  251. $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
  252. $ XNORM
  253. * ..
  254. * .. External Functions ..
  255. LOGICAL LSAME
  256. INTEGER IDAMAX
  257. REAL*8 DDOT, DLAMCH
  258. EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
  259. * ..
  260. * .. External Subroutines ..
  261. * ..
  262. ** .. Intrinsic Functions ..
  263. * INTRINSIC ABS, MAX, SQRT
  264. ** ..
  265. ** .. Local Arrays ..
  266. REAL*8 X( 2, 2 )
  267. ** ..
  268. ** .. Executable Statements ..
  269. *
  270. * Decode and test the input parameters
  271. *
  272. BOTHV = LSAME( SIDE, 'B' )
  273. RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
  274. LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
  275. *
  276. ALLV = LSAME( HOWMNY, 'A' )
  277. OVER = LSAME( HOWMNY, 'B' )
  278. SOMEV = LSAME( HOWMNY, 'S' )
  279. *
  280. INFO = 0
  281. IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
  282. INFO = -1
  283. ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
  284. INFO = -2
  285. ELSE IF( N.LT.0 ) THEN
  286. INFO = -4
  287. ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  288. INFO = -6
  289. ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
  290. INFO = -8
  291. ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
  292. INFO = -10
  293. ELSE
  294. *
  295. * Set M to the number of columns required to store the selected
  296. * eigenvectors, standardize the array SELECT if necessary, and
  297. * test MM.
  298. *
  299. IF( SOMEV ) THEN
  300. M = 0
  301. PAIR = .FALSE.
  302. DO 10 J = 1, N
  303. IF( PAIR ) THEN
  304. PAIR = .FALSE.
  305. SELECT( J ) = .FALSE.
  306. ELSE
  307. IF( J.LT.N ) THEN
  308. IF( T( J+1, J ).EQ.ZERO ) THEN
  309. IF( SELECT( J ) )
  310. $ M = M + 1
  311. ELSE
  312. PAIR = .TRUE.
  313. IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
  314. SELECT( J ) = .TRUE.
  315. M = M + 2
  316. END IF
  317. END IF
  318. ELSE
  319. IF( SELECT( N ) )
  320. $ M = M + 1
  321. END IF
  322. END IF
  323. 10 CONTINUE
  324. ELSE
  325. M = N
  326. END IF
  327. *
  328. IF( MM.LT.M ) THEN
  329. INFO = -11
  330. END IF
  331. END IF
  332. IF( INFO.NE.0 ) THEN
  333. CALL XERBLA( 'DTREVC', -INFO )
  334. RETURN
  335. END IF
  336. *
  337. * Quick return if possible.
  338. *
  339. IF( N.EQ.0 )
  340. $ RETURN
  341. *
  342. * Set the constants to control overflow.
  343. *
  344. UNFL = DLAMCH( 'Safe minimum' )
  345. OVFL = ONE / UNFL
  346. CALL DLABAD( UNFL, OVFL )
  347. ULP = DLAMCH( 'Precision' )
  348. SMLNUM = UNFL*( N / ULP )
  349. BIGNUM = ( ONE-ULP ) / SMLNUM
  350. *
  351. * Compute 1-norm of each column of strictly upper triangular
  352. * part of T to control overflow in triangular solver.
  353. *
  354. WORK( 1 ) = ZERO
  355. DO 30 J = 2, N
  356. WORK( J ) = ZERO
  357. DO 20 I = 1, J - 1
  358. WORK( J ) = WORK( J ) + ABS( T( I, J ) )
  359. 20 CONTINUE
  360. 30 CONTINUE
  361. *
  362. * Index IP is used to specify the real or complex eigenvalue:
  363. * IP = 0, real eigenvalue,
  364. * 1, first of conjugate complex pair: (wr,wi)
  365. * -1, second of conjugate complex pair: (wr,wi)
  366. *
  367. N2 = 2*N
  368. *
  369. IF( RIGHTV ) THEN
  370. *
  371. * Compute right eigenvectors.
  372. *
  373. IP = 0
  374. IS = M
  375. DO 140 KI = N, 1, -1
  376. *
  377. IF( IP.EQ.1 )
  378. $ GO TO 130
  379. IF( KI.EQ.1 )
  380. $ GO TO 40
  381. IF( T( KI, KI-1 ).EQ.ZERO )
  382. $ GO TO 40
  383. IP = -1
  384. *
  385. 40 CONTINUE
  386. IF( SOMEV ) THEN
  387. IF( IP.EQ.0 ) THEN
  388. IF( .NOT.SELECT( KI ) )
  389. $ GO TO 130
  390. ELSE
  391. IF( .NOT.SELECT( KI-1 ) )
  392. $ GO TO 130
  393. END IF
  394. END IF
  395. *
  396. * Compute the KI-th eigenvalue (WR,WI).
  397. *
  398. WR = T( KI, KI )
  399. WI = ZERO
  400. IF( IP.NE.0 )
  401. $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
  402. $ SQRT( ABS( T( KI-1, KI ) ) )
  403. SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
  404. *
  405. IF( IP.EQ.0 ) THEN
  406. *
  407. * Real right eigenvector
  408. *
  409. WORK( KI+N ) = ONE
  410. *
  411. * Form right-hand side
  412. *
  413. DO 50 K = 1, KI - 1
  414. WORK( K+N ) = -T( K, KI )
  415. 50 CONTINUE
  416. *
  417. * Solve the upper quasi-triangular system:
  418. * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
  419. *
  420. JNXT = KI - 1
  421. DO 60 J = KI - 1, 1, -1
  422. IF( J.GT.JNXT )
  423. $ GO TO 60
  424. J1 = J
  425. J2 = J
  426. JNXT = J - 1
  427. IF( J.GT.1 ) THEN
  428. IF( T( J, J-1 ).NE.ZERO ) THEN
  429. J1 = J - 1
  430. JNXT = J - 2
  431. END IF
  432. END IF
  433. *
  434. IF( J1.EQ.J2 ) THEN
  435. *
  436. * 1-by-1 diagonal block
  437. *
  438. CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
  439. $ LDT, ONE, ONE, WORK( J+N ), N, WR,
  440. $ ZERO, X, 2, SCALE, XNORM, IERR )
  441. *
  442. * Scale X(1,1) to avoid overflow when updating
  443. * the right-hand side.
  444. *
  445. IF( XNORM.GT.ONE ) THEN
  446. IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
  447. X( 1, 1 ) = X( 1, 1 ) / XNORM
  448. SCALE = SCALE / XNORM
  449. END IF
  450. END IF
  451. *
  452. * Scale if necessary
  453. *
  454. IF( SCALE.NE.ONE )
  455. $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
  456. WORK( J+N ) = X( 1, 1 )
  457. *
  458. * Update right-hand side
  459. *
  460. CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
  461. $ WORK( 1+N ), 1 )
  462. *
  463. ELSE
  464. *
  465. * 2-by-2 diagonal block
  466. *
  467. CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
  468. $ T( J-1, J-1 ), LDT, ONE, ONE,
  469. $ WORK( J-1+N ), N, WR, ZERO, X, 2,
  470. $ SCALE, XNORM, IERR )
  471. *
  472. * Scale X(1,1) and X(2,1) to avoid overflow when
  473. * updating the right-hand side.
  474. *
  475. IF( XNORM.GT.ONE ) THEN
  476. BETA = MAX( WORK( J-1 ), WORK( J ) )
  477. IF( BETA.GT.BIGNUM / XNORM ) THEN
  478. X( 1, 1 ) = X( 1, 1 ) / XNORM
  479. X( 2, 1 ) = X( 2, 1 ) / XNORM
  480. SCALE = SCALE / XNORM
  481. END IF
  482. END IF
  483. *
  484. * Scale if necessary
  485. *
  486. IF( SCALE.NE.ONE )
  487. $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
  488. WORK( J-1+N ) = X( 1, 1 )
  489. WORK( J+N ) = X( 2, 1 )
  490. *
  491. * Update right-hand side
  492. *
  493. CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
  494. $ WORK( 1+N ), 1 )
  495. CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
  496. $ WORK( 1+N ), 1 )
  497. END IF
  498. 60 CONTINUE
  499. *
  500. * Copy the vector x or Q*x to VR and normalize.
  501. *
  502. IF( .NOT.OVER ) THEN
  503. CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
  504. *
  505. II = IDAMAX( KI, VR( 1, IS ), 1 )
  506. REMAX = ONE / ABS( VR( II, IS ) )
  507. CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
  508. *
  509. DO 70 K = KI + 1, N
  510. VR( K, IS ) = ZERO
  511. 70 CONTINUE
  512. ELSE
  513. IF( KI.GT.1 )
  514. $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
  515. $ WORK( 1+N ), 1, WORK( KI+N ),
  516. $ VR( 1, KI ), 1 )
  517. *
  518. II = IDAMAX( N, VR( 1, KI ), 1 )
  519. REMAX = ONE / ABS( VR( II, KI ) )
  520. CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
  521. END IF
  522. *
  523. ELSE
  524. *
  525. * Complex right eigenvector.
  526. *
  527. * Initial solve
  528. * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
  529. * [ (T(KI,KI-1) T(KI,KI) ) ]
  530. *
  531. IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
  532. WORK( KI-1+N ) = ONE
  533. WORK( KI+N2 ) = WI / T( KI-1, KI )
  534. ELSE
  535. WORK( KI-1+N ) = -WI / T( KI, KI-1 )
  536. WORK( KI+N2 ) = ONE
  537. END IF
  538. WORK( KI+N ) = ZERO
  539. WORK( KI-1+N2 ) = ZERO
  540. *
  541. * Form right-hand side
  542. *
  543. DO 80 K = 1, KI - 2
  544. WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
  545. WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
  546. 80 CONTINUE
  547. *
  548. * Solve upper quasi-triangular system:
  549. * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
  550. *
  551. JNXT = KI - 2
  552. DO 90 J = KI - 2, 1, -1
  553. IF( J.GT.JNXT )
  554. $ GO TO 90
  555. J1 = J
  556. J2 = J
  557. JNXT = J - 1
  558. IF( J.GT.1 ) THEN
  559. IF( T( J, J-1 ).NE.ZERO ) THEN
  560. J1 = J - 1
  561. JNXT = J - 2
  562. END IF
  563. END IF
  564. *
  565. IF( J1.EQ.J2 ) THEN
  566. *
  567. * 1-by-1 diagonal block
  568. *
  569. CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
  570. $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
  571. $ X, 2, SCALE, XNORM, IERR )
  572. *
  573. * Scale X(1,1) and X(1,2) to avoid overflow when
  574. * updating the right-hand side.
  575. *
  576. IF( XNORM.GT.ONE ) THEN
  577. IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
  578. X( 1, 1 ) = X( 1, 1 ) / XNORM
  579. X( 1, 2 ) = X( 1, 2 ) / XNORM
  580. SCALE = SCALE / XNORM
  581. END IF
  582. END IF
  583. *
  584. * Scale if necessary
  585. *
  586. IF( SCALE.NE.ONE ) THEN
  587. CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
  588. CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
  589. END IF
  590. WORK( J+N ) = X( 1, 1 )
  591. WORK( J+N2 ) = X( 1, 2 )
  592. *
  593. * Update the right-hand side
  594. *
  595. CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
  596. $ WORK( 1+N ), 1 )
  597. CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
  598. $ WORK( 1+N2 ), 1 )
  599. *
  600. ELSE
  601. *
  602. * 2-by-2 diagonal block
  603. *
  604. CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
  605. $ T( J-1, J-1 ), LDT, ONE, ONE,
  606. $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
  607. $ XNORM, IERR )
  608. *
  609. * Scale X to avoid overflow when updating
  610. * the right-hand side.
  611. *
  612. IF( XNORM.GT.ONE ) THEN
  613. BETA = MAX( WORK( J-1 ), WORK( J ) )
  614. IF( BETA.GT.BIGNUM / XNORM ) THEN
  615. REC = ONE / XNORM
  616. X( 1, 1 ) = X( 1, 1 )*REC
  617. X( 1, 2 ) = X( 1, 2 )*REC
  618. X( 2, 1 ) = X( 2, 1 )*REC
  619. X( 2, 2 ) = X( 2, 2 )*REC
  620. SCALE = SCALE*REC
  621. END IF
  622. END IF
  623. *
  624. * Scale if necessary
  625. *
  626. IF( SCALE.NE.ONE ) THEN
  627. CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
  628. CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
  629. END IF
  630. WORK( J-1+N ) = X( 1, 1 )
  631. WORK( J+N ) = X( 2, 1 )
  632. WORK( J-1+N2 ) = X( 1, 2 )
  633. WORK( J+N2 ) = X( 2, 2 )
  634. *
  635. * Update the right-hand side
  636. *
  637. CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
  638. $ WORK( 1+N ), 1 )
  639. CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
  640. $ WORK( 1+N ), 1 )
  641. CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
  642. $ WORK( 1+N2 ), 1 )
  643. CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
  644. $ WORK( 1+N2 ), 1 )
  645. END IF
  646. 90 CONTINUE
  647. *
  648. * Copy the vector x or Q*x to VR and normalize.
  649. *
  650. IF( .NOT.OVER ) THEN
  651. CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
  652. CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
  653. *
  654. EMAX = ZERO
  655. DO 100 K = 1, KI
  656. EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
  657. $ ABS( VR( K, IS ) ) )
  658. 100 CONTINUE
  659. *
  660. REMAX = ONE / EMAX
  661. CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
  662. CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
  663. *
  664. DO 110 K = KI + 1, N
  665. VR( K, IS-1 ) = ZERO
  666. VR( K, IS ) = ZERO
  667. 110 CONTINUE
  668. *
  669. ELSE
  670. *
  671. IF( KI.GT.2 ) THEN
  672. CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
  673. $ WORK( 1+N ), 1, WORK( KI-1+N ),
  674. $ VR( 1, KI-1 ), 1 )
  675. CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
  676. $ WORK( 1+N2 ), 1, WORK( KI+N2 ),
  677. $ VR( 1, KI ), 1 )
  678. ELSE
  679. CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
  680. CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
  681. END IF
  682. *
  683. EMAX = ZERO
  684. DO 120 K = 1, N
  685. EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
  686. $ ABS( VR( K, KI ) ) )
  687. 120 CONTINUE
  688. REMAX = ONE / EMAX
  689. CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
  690. CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
  691. END IF
  692. END IF
  693. *
  694. IS = IS - 1
  695. IF( IP.NE.0 )
  696. $ IS = IS - 1
  697. 130 CONTINUE
  698. IF( IP.EQ.1 )
  699. $ IP = 0
  700. IF( IP.EQ.-1 )
  701. $ IP = 1
  702. 140 CONTINUE
  703. END IF
  704. *
  705. IF( LEFTV ) THEN
  706. *
  707. * Compute left eigenvectors.
  708. *
  709. IP = 0
  710. IS = 1
  711. DO 260 KI = 1, N
  712. *
  713. IF( IP.EQ.-1 )
  714. $ GO TO 250
  715. IF( KI.EQ.N )
  716. $ GO TO 150
  717. IF( T( KI+1, KI ).EQ.ZERO )
  718. $ GO TO 150
  719. IP = 1
  720. *
  721. 150 CONTINUE
  722. IF( SOMEV ) THEN
  723. IF( .NOT.SELECT( KI ) )
  724. $ GO TO 250
  725. END IF
  726. *
  727. * Compute the KI-th eigenvalue (WR,WI).
  728. *
  729. WR = T( KI, KI )
  730. WI = ZERO
  731. IF( IP.NE.0 )
  732. $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
  733. $ SQRT( ABS( T( KI+1, KI ) ) )
  734. SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
  735. *
  736. IF( IP.EQ.0 ) THEN
  737. *
  738. * Real left eigenvector.
  739. *
  740. WORK( KI+N ) = ONE
  741. *
  742. * Form right-hand side
  743. *
  744. DO 160 K = KI + 1, N
  745. WORK( K+N ) = -T( KI, K )
  746. 160 CONTINUE
  747. *
  748. * Solve the quasi-triangular system:
  749. * (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
  750. *
  751. VMAX = ONE
  752. VCRIT = BIGNUM
  753. *
  754. JNXT = KI + 1
  755. DO 170 J = KI + 1, N
  756. IF( J.LT.JNXT )
  757. $ GO TO 170
  758. J1 = J
  759. J2 = J
  760. JNXT = J + 1
  761. IF( J.LT.N ) THEN
  762. IF( T( J+1, J ).NE.ZERO ) THEN
  763. J2 = J + 1
  764. JNXT = J + 2
  765. END IF
  766. END IF
  767. *
  768. IF( J1.EQ.J2 ) THEN
  769. *
  770. * 1-by-1 diagonal block
  771. *
  772. * Scale if necessary to avoid overflow when forming
  773. * the right-hand side.
  774. *
  775. IF( WORK( J ).GT.VCRIT ) THEN
  776. REC = ONE / VMAX
  777. CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
  778. VMAX = ONE
  779. VCRIT = BIGNUM
  780. END IF
  781. *
  782. WORK( J+N ) = WORK( J+N ) -
  783. $ DDOT( J-KI-1, T( KI+1, J ), 1,
  784. $ WORK( KI+1+N ), 1 )
  785. *
  786. * Solve (T(J,J)-WR)**T*X = WORK
  787. *
  788. CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
  789. $ LDT, ONE, ONE, WORK( J+N ), N, WR,
  790. $ ZERO, X, 2, SCALE, XNORM, IERR )
  791. *
  792. * Scale if necessary
  793. *
  794. IF( SCALE.NE.ONE )
  795. $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
  796. WORK( J+N ) = X( 1, 1 )
  797. VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
  798. VCRIT = BIGNUM / VMAX
  799. *
  800. ELSE
  801. *
  802. * 2-by-2 diagonal block
  803. *
  804. * Scale if necessary to avoid overflow when forming
  805. * the right-hand side.
  806. *
  807. BETA = MAX( WORK( J ), WORK( J+1 ) )
  808. IF( BETA.GT.VCRIT ) THEN
  809. REC = ONE / VMAX
  810. CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
  811. VMAX = ONE
  812. VCRIT = BIGNUM
  813. END IF
  814. *
  815. WORK( J+N ) = WORK( J+N ) -
  816. $ DDOT( J-KI-1, T( KI+1, J ), 1,
  817. $ WORK( KI+1+N ), 1 )
  818. *
  819. WORK( J+1+N ) = WORK( J+1+N ) -
  820. $ DDOT( J-KI-1, T( KI+1, J+1 ), 1,
  821. $ WORK( KI+1+N ), 1 )
  822. *
  823. * Solve
  824. * [T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
  825. * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
  826. *
  827. CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
  828. $ LDT, ONE, ONE, WORK( J+N ), N, WR,
  829. $ ZERO, X, 2, SCALE, XNORM, IERR )
  830. *
  831. * Scale if necessary
  832. *
  833. IF( SCALE.NE.ONE )
  834. $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
  835. WORK( J+N ) = X( 1, 1 )
  836. WORK( J+1+N ) = X( 2, 1 )
  837. *
  838. VMAX = MAX( ABS( WORK( J+N ) ),
  839. $ ABS( WORK( J+1+N ) ), VMAX )
  840. VCRIT = BIGNUM / VMAX
  841. *
  842. END IF
  843. 170 CONTINUE
  844. *
  845. * Copy the vector x or Q*x to VL and normalize.
  846. *
  847. IF( .NOT.OVER ) THEN
  848. CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
  849. *
  850. II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
  851. REMAX = ONE / ABS( VL( II, IS ) )
  852. CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
  853. *
  854. DO 180 K = 1, KI - 1
  855. VL( K, IS ) = ZERO
  856. 180 CONTINUE
  857. *
  858. ELSE
  859. *
  860. IF( KI.LT.N )
  861. $ CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
  862. $ WORK( KI+1+N ), 1, WORK( KI+N ),
  863. $ VL( 1, KI ), 1 )
  864. *
  865. II = IDAMAX( N, VL( 1, KI ), 1 )
  866. REMAX = ONE / ABS( VL( II, KI ) )
  867. CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
  868. *
  869. END IF
  870. *
  871. ELSE
  872. *
  873. * Complex left eigenvector.
  874. *
  875. * Initial solve:
  876. * ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
  877. * ((T(KI+1,KI) T(KI+1,KI+1)) )
  878. *
  879. IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
  880. WORK( KI+N ) = WI / T( KI, KI+1 )
  881. WORK( KI+1+N2 ) = ONE
  882. ELSE
  883. WORK( KI+N ) = ONE
  884. WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
  885. END IF
  886. WORK( KI+1+N ) = ZERO
  887. WORK( KI+N2 ) = ZERO
  888. *
  889. * Form right-hand side
  890. *
  891. DO 190 K = KI + 2, N
  892. WORK( K+N ) = -WORK( KI+N )*T( KI, K )
  893. WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
  894. 190 CONTINUE
  895. *
  896. * Solve complex quasi-triangular system:
  897. * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
  898. *
  899. VMAX = ONE
  900. VCRIT = BIGNUM
  901. *
  902. JNXT = KI + 2
  903. DO 200 J = KI + 2, N
  904. IF( J.LT.JNXT )
  905. $ GO TO 200
  906. J1 = J
  907. J2 = J
  908. JNXT = J + 1
  909. IF( J.LT.N ) THEN
  910. IF( T( J+1, J ).NE.ZERO ) THEN
  911. J2 = J + 1
  912. JNXT = J + 2
  913. END IF
  914. END IF
  915. *
  916. IF( J1.EQ.J2 ) THEN
  917. *
  918. * 1-by-1 diagonal block
  919. *
  920. * Scale if necessary to avoid overflow when
  921. * forming the right-hand side elements.
  922. *
  923. IF( WORK( J ).GT.VCRIT ) THEN
  924. REC = ONE / VMAX
  925. CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
  926. CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
  927. VMAX = ONE
  928. VCRIT = BIGNUM
  929. END IF
  930. *
  931. WORK( J+N ) = WORK( J+N ) -
  932. $ DDOT( J-KI-2, T( KI+2, J ), 1,
  933. $ WORK( KI+2+N ), 1 )
  934. WORK( J+N2 ) = WORK( J+N2 ) -
  935. $ DDOT( J-KI-2, T( KI+2, J ), 1,
  936. $ WORK( KI+2+N2 ), 1 )
  937. *
  938. * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
  939. *
  940. CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
  941. $ LDT, ONE, ONE, WORK( J+N ), N, WR,
  942. $ -WI, X, 2, SCALE, XNORM, IERR )
  943. *
  944. * Scale if necessary
  945. *
  946. IF( SCALE.NE.ONE ) THEN
  947. CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
  948. CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
  949. END IF
  950. WORK( J+N ) = X( 1, 1 )
  951. WORK( J+N2 ) = X( 1, 2 )
  952. VMAX = MAX( ABS( WORK( J+N ) ),
  953. $ ABS( WORK( J+N2 ) ), VMAX )
  954. VCRIT = BIGNUM / VMAX
  955. *
  956. ELSE
  957. *
  958. * 2-by-2 diagonal block
  959. *
  960. * Scale if necessary to avoid overflow when forming
  961. * the right-hand side elements.
  962. *
  963. BETA = MAX( WORK( J ), WORK( J+1 ) )
  964. IF( BETA.GT.VCRIT ) THEN
  965. REC = ONE / VMAX
  966. CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
  967. CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
  968. VMAX = ONE
  969. VCRIT = BIGNUM
  970. END IF
  971. *
  972. WORK( J+N ) = WORK( J+N ) -
  973. $ DDOT( J-KI-2, T( KI+2, J ), 1,
  974. $ WORK( KI+2+N ), 1 )
  975. *
  976. WORK( J+N2 ) = WORK( J+N2 ) -
  977. $ DDOT( J-KI-2, T( KI+2, J ), 1,
  978. $ WORK( KI+2+N2 ), 1 )
  979. *
  980. WORK( J+1+N ) = WORK( J+1+N ) -
  981. $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
  982. $ WORK( KI+2+N ), 1 )
  983. *
  984. WORK( J+1+N2 ) = WORK( J+1+N2 ) -
  985. $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
  986. $ WORK( KI+2+N2 ), 1 )
  987. *
  988. * Solve 2-by-2 complex linear equation
  989. * ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B
  990. * ([T(j+1,j) T(j+1,j+1)] )
  991. *
  992. CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
  993. $ LDT, ONE, ONE, WORK( J+N ), N, WR,
  994. $ -WI, X, 2, SCALE, XNORM, IERR )
  995. *
  996. * Scale if necessary
  997. *
  998. IF( SCALE.NE.ONE ) THEN
  999. CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
  1000. CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
  1001. END IF
  1002. WORK( J+N ) = X( 1, 1 )
  1003. WORK( J+N2 ) = X( 1, 2 )
  1004. WORK( J+1+N ) = X( 2, 1 )
  1005. WORK( J+1+N2 ) = X( 2, 2 )
  1006. VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
  1007. $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
  1008. VCRIT = BIGNUM / VMAX
  1009. *
  1010. END IF
  1011. 200 CONTINUE
  1012. *
  1013. * Copy the vector x or Q*x to VL and normalize.
  1014. *
  1015. IF( .NOT.OVER ) THEN
  1016. CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
  1017. CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
  1018. $ 1 )
  1019. *
  1020. EMAX = ZERO
  1021. DO 220 K = KI, N
  1022. EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
  1023. $ ABS( VL( K, IS+1 ) ) )
  1024. 220 CONTINUE
  1025. REMAX = ONE / EMAX
  1026. CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
  1027. CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
  1028. *
  1029. DO 230 K = 1, KI - 1
  1030. VL( K, IS ) = ZERO
  1031. VL( K, IS+1 ) = ZERO
  1032. 230 CONTINUE
  1033. ELSE
  1034. IF( KI.LT.N-1 ) THEN
  1035. CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
  1036. $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
  1037. $ VL( 1, KI ), 1 )
  1038. CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
  1039. $ LDVL, WORK( KI+2+N2 ), 1,
  1040. $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
  1041. ELSE
  1042. CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
  1043. CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
  1044. END IF
  1045. *
  1046. EMAX = ZERO
  1047. DO 240 K = 1, N
  1048. EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
  1049. $ ABS( VL( K, KI+1 ) ) )
  1050. 240 CONTINUE
  1051. REMAX = ONE / EMAX
  1052. CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
  1053. CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
  1054. *
  1055. END IF
  1056. *
  1057. END IF
  1058. *
  1059. IS = IS + 1
  1060. IF( IP.NE.0 )
  1061. $ IS = IS + 1
  1062. 250 CONTINUE
  1063. IF( IP.EQ.-1 )
  1064. $ IP = 0
  1065. IF( IP.EQ.1 )
  1066. $ IP = -1
  1067. *
  1068. 260 CONTINUE
  1069. *
  1070. END IF
  1071. *
  1072. RETURN
  1073. *
  1074. * End of DTREVC
  1075. *
  1076. END
  1077.  
  1078.  

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