Télécharger dtrevc3.eso

Retour à la liste

Numérotation des lignes :

  1. C DTREVC3 SOURCE BP208322 20/09/18 21:16:12 10718
  2. *> \brief \b DTREVC3
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DTREVC3 + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc3.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc3.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc3.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
  23. * VR, LDVR, MM, M, WORK, LWORK, INFO )
  24. *
  25. * .. Scalar Arguments ..
  26. * CHARACTER HOWMNY, SIDE
  27. * INTEGER INFO, LDT, LDVL, LDVR, LWORK, 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. *> DTREVC3 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**H)*T = w*(y**H)
  50. *>
  51. *> where y**H denotes the conjugate 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. *>
  61. *> This uses a Level 3 BLAS version of the back transformation.
  62. *> \endverbatim
  63. *
  64. * Arguments:
  65. * ==========
  66. *
  67. *> \param[in] SIDE
  68. *> \verbatim
  69. *> SIDE is CHARACTER*1
  70. *> = 'R': compute right eigenvectors only;
  71. *> = 'L': compute left eigenvectors only;
  72. *> = 'B': compute both right and left eigenvectors.
  73. *> \endverbatim
  74. *>
  75. *> \param[in] HOWMNY
  76. *> \verbatim
  77. *> HOWMNY is CHARACTER*1
  78. *> = 'A': compute all right and/or left eigenvectors;
  79. *> = 'B': compute all right and/or left eigenvectors,
  80. *> backtransformed by the matrices in VR and/or VL;
  81. *> = 'S': compute selected right and/or left eigenvectors,
  82. *> as indicated by the logical array SELECT.
  83. *> \endverbatim
  84. *>
  85. *> \param[in,out] SELECT
  86. *> \verbatim
  87. *> SELECT is LOGICAL array, dimension (N)
  88. *> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
  89. *> computed.
  90. *> If w(j) is a real eigenvalue, the corresponding real
  91. *> eigenvector is computed if SELECT(j) is .TRUE..
  92. *> If w(j) and w(j+1) are the real and imaginary parts of a
  93. *> complex eigenvalue, the corresponding complex eigenvector is
  94. *> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
  95. *> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
  96. *> .FALSE..
  97. *> Not referenced if HOWMNY = 'A' or 'B'.
  98. *> \endverbatim
  99. *>
  100. *> \param[in] N
  101. *> \verbatim
  102. *> N is INTEGER
  103. *> The order of the matrix T. N >= 0.
  104. *> \endverbatim
  105. *>
  106. *> \param[in] T
  107. *> \verbatim
  108. *> T is REAL*8 array, dimension (LDT,N)
  109. *> The upper quasi-triangular matrix T in Schur canonical form.
  110. *> \endverbatim
  111. *>
  112. *> \param[in] LDT
  113. *> \verbatim
  114. *> LDT is INTEGER
  115. *> The leading dimension of the array T. LDT >= max(1,N).
  116. *> \endverbatim
  117. *>
  118. *> \param[in,out] VL
  119. *> \verbatim
  120. *> VL is REAL*8 array, dimension (LDVL,MM)
  121. *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
  122. *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
  123. *> of Schur vectors returned by DHSEQR).
  124. *> On exit, if SIDE = 'L' or 'B', VL contains:
  125. *> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
  126. *> if HOWMNY = 'B', the matrix Q*Y;
  127. *> if HOWMNY = 'S', the left eigenvectors of T specified by
  128. *> SELECT, stored consecutively in the columns
  129. *> of VL, in the same order as their
  130. *> eigenvalues.
  131. *> A complex eigenvector corresponding to a complex eigenvalue
  132. *> is stored in two consecutive columns, the first holding the
  133. *> real part, and the second the imaginary part.
  134. *> Not referenced if SIDE = 'R'.
  135. *> \endverbatim
  136. *>
  137. *> \param[in] LDVL
  138. *> \verbatim
  139. *> LDVL is INTEGER
  140. *> The leading dimension of the array VL.
  141. *> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
  142. *> \endverbatim
  143. *>
  144. *> \param[in,out] VR
  145. *> \verbatim
  146. *> VR is REAL*8 array, dimension (LDVR,MM)
  147. *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
  148. *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
  149. *> of Schur vectors returned by DHSEQR).
  150. *> On exit, if SIDE = 'R' or 'B', VR contains:
  151. *> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
  152. *> if HOWMNY = 'B', the matrix Q*X;
  153. *> if HOWMNY = 'S', the right eigenvectors of T specified by
  154. *> SELECT, stored consecutively in the columns
  155. *> of VR, in the same order as their
  156. *> eigenvalues.
  157. *> A complex eigenvector corresponding to a complex eigenvalue
  158. *> is stored in two consecutive columns, the first holding the
  159. *> real part and the second the imaginary part.
  160. *> Not referenced if SIDE = 'L'.
  161. *> \endverbatim
  162. *>
  163. *> \param[in] LDVR
  164. *> \verbatim
  165. *> LDVR is INTEGER
  166. *> The leading dimension of the array VR.
  167. *> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
  168. *> \endverbatim
  169. *>
  170. *> \param[in] MM
  171. *> \verbatim
  172. *> MM is INTEGER
  173. *> The number of columns in the arrays VL and/or VR. MM >= M.
  174. *> \endverbatim
  175. *>
  176. *> \param[out] M
  177. *> \verbatim
  178. *> M is INTEGER
  179. *> The number of columns in the arrays VL and/or VR actually
  180. *> used to store the eigenvectors.
  181. *> If HOWMNY = 'A' or 'B', M is set to N.
  182. *> Each selected real eigenvector occupies one column and each
  183. *> selected complex eigenvector occupies two columns.
  184. *> \endverbatim
  185. *>
  186. *> \param[out] WORK
  187. *> \verbatim
  188. *> WORK is REAL*8 array, dimension (MAX(1,LWORK))
  189. *> \endverbatim
  190. *>
  191. *> \param[in] LWORK
  192. *> \verbatim
  193. *> LWORK is INTEGER
  194. *> The dimension of array WORK. LWORK >= max(1,3*N).
  195. *> For optimum performance, LWORK >= N + 2*N*NB, where NB is
  196. *> the optimal blocksize.
  197. *>
  198. *> If LWORK = -1, then a workspace query is assumed; the routine
  199. *> only calculates the optimal size of the WORK array, returns
  200. *> this value as the first entry of the WORK array, and no error
  201. *> message related to LWORK is issued by XERBLA.
  202. *> \endverbatim
  203. *>
  204. *> \param[out] INFO
  205. *> \verbatim
  206. *> INFO is INTEGER
  207. *> = 0: successful exit
  208. *> < 0: if INFO = -i, the i-th argument had an illegal value
  209. *> \endverbatim
  210. *
  211. * Authors:
  212. * ========
  213. *
  214. *> \author Univ. of Tennessee
  215. *> \author Univ. of California Berkeley
  216. *> \author Univ. of Colorado Denver
  217. *> \author NAG Ltd.
  218. *
  219. *> \date November 2017
  220. *
  221. * @precisions fortran d -> s
  222. *
  223. *> \ingroup doubleOTHERcomputational
  224. *
  225. *> \par Further Details:
  226. * =====================
  227. *>
  228. *> \verbatim
  229. *>
  230. *> The algorithm used in this program is basically backward (forward)
  231. *> substitution, with scaling to make the the code robust against
  232. *> possible overflow.
  233. *>
  234. *> Each eigenvector is normalized so that the element of largest
  235. *> magnitude has magnitude 1; here the magnitude of a complex number
  236. *> (x,y) is taken to be |x| + |y|.
  237. *> \endverbatim
  238. *>
  239. * =====================================================================
  240. SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
  241. $ VR, LDVR, MM, M, WORK, LWORK, INFO )
  242. *
  243. * -- LAPACK computational routine (version 3.8.0) --
  244. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  245. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  246. * November 2017
  247.  
  248. * IMPLICIT NONE
  249. IMPLICIT INTEGER(I-N)
  250. IMPLICIT REAL*8(A-H,O-Z)
  251. *
  252. * .. Scalar Arguments ..
  253. CHARACTER HOWMNY, SIDE
  254. INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
  255. * ..
  256. * .. Array Arguments ..
  257. LOGICAL SELECT( * )
  258. REAL*8 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  259. $ WORK( * )
  260. * ..
  261. *
  262. * =====================================================================
  263. *
  264. * .. Parameters ..
  265. REAL*8 ZERO, ONE
  266. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  267. INTEGER NBMIN, NBMAX
  268. PARAMETER ( NBMIN = 8, NBMAX = 128 )
  269. * ..
  270. * .. Local Scalars ..
  271. LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
  272. $ RIGHTV, SOMEV
  273. INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
  274. $ IV, MAXWRK, NB, KI2
  275. REAL*8 BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
  276. $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
  277. $ XNORM
  278. * ..
  279. * .. External Functions ..
  280. LOGICAL LSAME
  281. INTEGER IDAMAX, ILAENV
  282. REAL*8 DDOT, DLAMCH
  283. * EXTERNAL LSAME, IDAMAX, ILAENV, DDOT, DLAMCH
  284. * ..
  285. * .. External Subroutines ..
  286. * EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA,
  287. * $ DGEMM, DLASET, DLABAD, DLACPY
  288. * ..
  289. * .. Intrinsic Functions ..
  290. * INTRINSIC ABS, MAX, SQRT
  291. * ..
  292. * .. Local Arrays ..
  293. REAL*8 X( 2, 2 )
  294. INTEGER ISCOMPLEX( NBMAX )
  295. * ..
  296. * .. Executable Statements ..
  297. *
  298. * Decode and test the input parameters
  299. *
  300. BOTHV = LSAME( SIDE, 'B' )
  301. RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
  302. LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
  303. *
  304. ALLV = LSAME( HOWMNY, 'A' )
  305. OVER = LSAME( HOWMNY, 'B' )
  306. SOMEV = LSAME( HOWMNY, 'S' )
  307. *
  308. INFO = 0
  309. NB = ILAENV( 1, 'DTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
  310. MAXWRK = N + 2*N*NB
  311. WORK(1) = MAXWRK
  312. LQUERY = ( LWORK.EQ.-1 )
  313. IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
  314. INFO = -1
  315. ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
  316. INFO = -2
  317. ELSE IF( N.LT.0 ) THEN
  318. INFO = -4
  319. ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  320. INFO = -6
  321. ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
  322. INFO = -8
  323. ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
  324. INFO = -10
  325. ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
  326. INFO = -14
  327. ELSE
  328. *
  329. * Set M to the number of columns required to store the selected
  330. * eigenvectors, standardize the array SELECT if necessary, and
  331. * test MM.
  332. *
  333. IF( SOMEV ) THEN
  334. M = 0
  335. PAIR = .FALSE.
  336. DO 10 J = 1, N
  337. IF( PAIR ) THEN
  338. PAIR = .FALSE.
  339. SELECT( J ) = .FALSE.
  340. ELSE
  341. IF( J.LT.N ) THEN
  342. IF( T( J+1, J ).EQ.ZERO ) THEN
  343. IF( SELECT( J ) )
  344. $ M = M + 1
  345. ELSE
  346. PAIR = .TRUE.
  347. IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
  348. SELECT( J ) = .TRUE.
  349. M = M + 2
  350. END IF
  351. END IF
  352. ELSE
  353. IF( SELECT( N ) )
  354. $ M = M + 1
  355. END IF
  356. END IF
  357. 10 CONTINUE
  358. ELSE
  359. M = N
  360. END IF
  361. *
  362. IF( MM.LT.M ) THEN
  363. INFO = -11
  364. END IF
  365. END IF
  366. IF( INFO.NE.0 ) THEN
  367. CALL XERBLA( 'DTREVC3', -INFO )
  368. RETURN
  369. ELSE IF( LQUERY ) THEN
  370. RETURN
  371. END IF
  372. *
  373. * Quick return if possible.
  374. *
  375. IF( N.EQ.0 )
  376. $ RETURN
  377. *
  378. * Use blocked version of back-transformation if sufficient workspace.
  379. * Zero-out the workspace to avoid potential NaN propagation.
  380. *
  381. IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
  382. NB = (LWORK - N) / (2*N)
  383. NB = MIN( NB, NBMAX )
  384. CALL DLASET( 'F', N, 1+2*NB, ZERO, ZERO, WORK, N )
  385. ELSE
  386. NB = 1
  387. END IF
  388. *
  389. * Set the constants to control overflow.
  390. *
  391. UNFL = DLAMCH( 'Safe minimum' )
  392. OVFL = ONE / UNFL
  393. CALL DLABAD( UNFL, OVFL )
  394. ULP = DLAMCH( 'Precision' )
  395. SMLNUM = UNFL*( N / ULP )
  396. BIGNUM = ( ONE-ULP ) / SMLNUM
  397. *
  398. * Compute 1-norm of each column of strictly upper triangular
  399. * part of T to control overflow in triangular solver.
  400. *
  401. WORK( 1 ) = ZERO
  402. DO 30 J = 2, N
  403. WORK( J ) = ZERO
  404. DO 20 I = 1, J - 1
  405. WORK( J ) = WORK( J ) + ABS( T( I, J ) )
  406. 20 CONTINUE
  407. 30 CONTINUE
  408. *
  409. * Index IP is used to specify the real or complex eigenvalue:
  410. * IP = 0, real eigenvalue,
  411. * 1, first of conjugate complex pair: (wr,wi)
  412. * -1, second of conjugate complex pair: (wr,wi)
  413. * ISCOMPLEX array stores IP for each column in current block.
  414. *
  415. IF( RIGHTV ) THEN
  416. *
  417. * ============================================================
  418. * Compute right eigenvectors.
  419. *
  420. * IV is index of column in current block.
  421. * For complex right vector, uses IV-1 for real part and IV for complex part.
  422. * Non-blocked version always uses IV=2;
  423. * blocked version starts with IV=NB, goes down to 1 or 2.
  424. * (Note the "0-th" column is used for 1-norms computed above.)
  425. IV = 2
  426. IF( NB.GT.2 ) THEN
  427. IV = NB
  428. END IF
  429.  
  430. IP = 0
  431. IS = M
  432. DO 140 KI = N, 1, -1
  433. IF( IP.EQ.-1 ) THEN
  434. * previous iteration (ki+1) was second of conjugate pair,
  435. * so this ki is first of conjugate pair; skip to end of loop
  436. IP = 1
  437. GO TO 140
  438. ELSE IF( KI.EQ.1 ) THEN
  439. * last column, so this ki must be real eigenvalue
  440. IP = 0
  441. ELSE IF( T( KI, KI-1 ).EQ.ZERO ) THEN
  442. * zero on sub-diagonal, so this ki is real eigenvalue
  443. IP = 0
  444. ELSE
  445. * non-zero on sub-diagonal, so this ki is second of conjugate pair
  446. IP = -1
  447. END IF
  448.  
  449. IF( SOMEV ) THEN
  450. IF( IP.EQ.0 ) THEN
  451. IF( .NOT.SELECT( KI ) )
  452. $ GO TO 140
  453. ELSE
  454. IF( .NOT.SELECT( KI-1 ) )
  455. $ GO TO 140
  456. END IF
  457. END IF
  458. *
  459. * Compute the KI-th eigenvalue (WR,WI).
  460. *
  461. WR = T( KI, KI )
  462. WI = ZERO
  463. IF( IP.NE.0 )
  464. $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
  465. $ SQRT( ABS( T( KI-1, KI ) ) )
  466. SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
  467. *
  468. IF( IP.EQ.0 ) THEN
  469. *
  470. * --------------------------------------------------------
  471. * Real right eigenvector
  472. *
  473. WORK( KI + IV*N ) = ONE
  474. *
  475. * Form right-hand side.
  476. *
  477. DO 50 K = 1, KI - 1
  478. WORK( K + IV*N ) = -T( K, KI )
  479. 50 CONTINUE
  480. *
  481. * Solve upper quasi-triangular system:
  482. * [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
  483. *
  484. JNXT = KI - 1
  485. DO 60 J = KI - 1, 1, -1
  486. IF( J.GT.JNXT )
  487. $ GO TO 60
  488. J1 = J
  489. J2 = J
  490. JNXT = J - 1
  491. IF( J.GT.1 ) THEN
  492. IF( T( J, J-1 ).NE.ZERO ) THEN
  493. J1 = J - 1
  494. JNXT = J - 2
  495. END IF
  496. END IF
  497. *
  498. IF( J1.EQ.J2 ) THEN
  499. *
  500. * 1-by-1 diagonal block
  501. *
  502. CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
  503. $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
  504. $ ZERO, X, 2, SCALE, XNORM, IERR )
  505. *
  506. * Scale X(1,1) to avoid overflow when updating
  507. * the right-hand side.
  508. *
  509. IF( XNORM.GT.ONE ) THEN
  510. IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
  511. X( 1, 1 ) = X( 1, 1 ) / XNORM
  512. SCALE = SCALE / XNORM
  513. END IF
  514. END IF
  515. *
  516. * Scale if necessary
  517. *
  518. IF( SCALE.NE.ONE )
  519. $ CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
  520. WORK( J+IV*N ) = X( 1, 1 )
  521. *
  522. * Update right-hand side
  523. *
  524. CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
  525. $ WORK( 1+IV*N ), 1 )
  526. *
  527. ELSE
  528. *
  529. * 2-by-2 diagonal block
  530. *
  531. CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
  532. $ T( J-1, J-1 ), LDT, ONE, ONE,
  533. $ WORK( J-1+IV*N ), N, WR, ZERO, X, 2,
  534. $ SCALE, XNORM, IERR )
  535. *
  536. * Scale X(1,1) and X(2,1) to avoid overflow when
  537. * updating the right-hand side.
  538. *
  539. IF( XNORM.GT.ONE ) THEN
  540. BETA = MAX( WORK( J-1 ), WORK( J ) )
  541. IF( BETA.GT.BIGNUM / XNORM ) THEN
  542. X( 1, 1 ) = X( 1, 1 ) / XNORM
  543. X( 2, 1 ) = X( 2, 1 ) / XNORM
  544. SCALE = SCALE / XNORM
  545. END IF
  546. END IF
  547. *
  548. * Scale if necessary
  549. *
  550. IF( SCALE.NE.ONE )
  551. $ CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
  552. WORK( J-1+IV*N ) = X( 1, 1 )
  553. WORK( J +IV*N ) = X( 2, 1 )
  554. *
  555. * Update right-hand side
  556. *
  557. CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
  558. $ WORK( 1+IV*N ), 1 )
  559. CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
  560. $ WORK( 1+IV*N ), 1 )
  561. END IF
  562. 60 CONTINUE
  563. *
  564. * Copy the vector x or Q*x to VR and normalize.
  565. *
  566. IF( .NOT.OVER ) THEN
  567. * ------------------------------
  568. * no back-transform: copy x to VR and normalize.
  569. CALL DCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
  570. *
  571. II = IDAMAX( KI, VR( 1, IS ), 1 )
  572. REMAX = ONE / ABS( VR( II, IS ) )
  573. CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
  574. *
  575. DO 70 K = KI + 1, N
  576. VR( K, IS ) = ZERO
  577. 70 CONTINUE
  578. *
  579. ELSE IF( NB.EQ.1 ) THEN
  580. * ------------------------------
  581. * version 1: back-transform each vector with GEMV, Q*x.
  582. IF( KI.GT.1 )
  583. $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
  584. $ WORK( 1 + IV*N ), 1, WORK( KI + IV*N ),
  585. $ VR( 1, KI ), 1 )
  586. *
  587. II = IDAMAX( N, VR( 1, KI ), 1 )
  588. REMAX = ONE / ABS( VR( II, KI ) )
  589. CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
  590. *
  591. ELSE
  592. * ------------------------------
  593. * version 2: back-transform block of vectors with GEMM
  594. * zero out below vector
  595. DO K = KI + 1, N
  596. WORK( K + IV*N ) = ZERO
  597. END DO
  598. ISCOMPLEX( IV ) = IP
  599. * back-transform and normalization is done below
  600. END IF
  601. ELSE
  602. *
  603. * --------------------------------------------------------
  604. * Complex right eigenvector.
  605. *
  606. * Initial solve
  607. * [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
  608. * [ ( T(KI, KI-1) T(KI, KI) ) ]
  609. *
  610. IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
  611. WORK( KI-1 + (IV-1)*N ) = ONE
  612. WORK( KI + (IV )*N ) = WI / T( KI-1, KI )
  613. ELSE
  614. WORK( KI-1 + (IV-1)*N ) = -WI / T( KI, KI-1 )
  615. WORK( KI + (IV )*N ) = ONE
  616. END IF
  617. WORK( KI + (IV-1)*N ) = ZERO
  618. WORK( KI-1 + (IV )*N ) = ZERO
  619. *
  620. * Form right-hand side.
  621. *
  622. DO 80 K = 1, KI - 2
  623. WORK( K+(IV-1)*N ) = -WORK( KI-1+(IV-1)*N )*T(K,KI-1)
  624. WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(K,KI )
  625. 80 CONTINUE
  626. *
  627. * Solve upper quasi-triangular system:
  628. * [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
  629. *
  630. JNXT = KI - 2
  631. DO 90 J = KI - 2, 1, -1
  632. IF( J.GT.JNXT )
  633. $ GO TO 90
  634. J1 = J
  635. J2 = J
  636. JNXT = J - 1
  637. IF( J.GT.1 ) THEN
  638. IF( T( J, J-1 ).NE.ZERO ) THEN
  639. J1 = J - 1
  640. JNXT = J - 2
  641. END IF
  642. END IF
  643. *
  644. IF( J1.EQ.J2 ) THEN
  645. *
  646. * 1-by-1 diagonal block
  647. *
  648. CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
  649. $ LDT, ONE, ONE, WORK( J+(IV-1)*N ), N,
  650. $ WR, WI, X, 2, SCALE, XNORM, IERR )
  651. *
  652. * Scale X(1,1) and X(1,2) to avoid overflow when
  653. * updating the right-hand side.
  654. *
  655. IF( XNORM.GT.ONE ) THEN
  656. IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
  657. X( 1, 1 ) = X( 1, 1 ) / XNORM
  658. X( 1, 2 ) = X( 1, 2 ) / XNORM
  659. SCALE = SCALE / XNORM
  660. END IF
  661. END IF
  662. *
  663. * Scale if necessary
  664. *
  665. IF( SCALE.NE.ONE ) THEN
  666. CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
  667. CALL DSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
  668. END IF
  669. WORK( J+(IV-1)*N ) = X( 1, 1 )
  670. WORK( J+(IV )*N ) = X( 1, 2 )
  671. *
  672. * Update the right-hand side
  673. *
  674. CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
  675. $ WORK( 1+(IV-1)*N ), 1 )
  676. CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
  677. $ WORK( 1+(IV )*N ), 1 )
  678. *
  679. ELSE
  680. *
  681. * 2-by-2 diagonal block
  682. *
  683. CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
  684. $ T( J-1, J-1 ), LDT, ONE, ONE,
  685. $ WORK( J-1+(IV-1)*N ), N, WR, WI, X, 2,
  686. $ SCALE, XNORM, IERR )
  687. *
  688. * Scale X to avoid overflow when updating
  689. * the right-hand side.
  690. *
  691. IF( XNORM.GT.ONE ) THEN
  692. BETA = MAX( WORK( J-1 ), WORK( J ) )
  693. IF( BETA.GT.BIGNUM / XNORM ) THEN
  694. REC = ONE / XNORM
  695. X( 1, 1 ) = X( 1, 1 )*REC
  696. X( 1, 2 ) = X( 1, 2 )*REC
  697. X( 2, 1 ) = X( 2, 1 )*REC
  698. X( 2, 2 ) = X( 2, 2 )*REC
  699. SCALE = SCALE*REC
  700. END IF
  701. END IF
  702. *
  703. * Scale if necessary
  704. *
  705. IF( SCALE.NE.ONE ) THEN
  706. CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
  707. CALL DSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
  708. END IF
  709. WORK( J-1+(IV-1)*N ) = X( 1, 1 )
  710. WORK( J +(IV-1)*N ) = X( 2, 1 )
  711. WORK( J-1+(IV )*N ) = X( 1, 2 )
  712. WORK( J +(IV )*N ) = X( 2, 2 )
  713. *
  714. * Update the right-hand side
  715. *
  716. CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
  717. $ WORK( 1+(IV-1)*N ), 1 )
  718. CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
  719. $ WORK( 1+(IV-1)*N ), 1 )
  720. CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
  721. $ WORK( 1+(IV )*N ), 1 )
  722. CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
  723. $ WORK( 1+(IV )*N ), 1 )
  724. END IF
  725. 90 CONTINUE
  726. *
  727. * Copy the vector x or Q*x to VR and normalize.
  728. *
  729. IF( .NOT.OVER ) THEN
  730. * ------------------------------
  731. * no back-transform: copy x to VR and normalize.
  732. CALL DCOPY( KI, WORK( 1+(IV-1)*N ), 1, VR(1,IS-1), 1 )
  733. CALL DCOPY( KI, WORK( 1+(IV )*N ), 1, VR(1,IS ), 1 )
  734. *
  735. EMAX = ZERO
  736. DO 100 K = 1, KI
  737. EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
  738. $ ABS( VR( K, IS ) ) )
  739. 100 CONTINUE
  740. REMAX = ONE / EMAX
  741. CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
  742. CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
  743. *
  744. DO 110 K = KI + 1, N
  745. VR( K, IS-1 ) = ZERO
  746. VR( K, IS ) = ZERO
  747. 110 CONTINUE
  748. *
  749. ELSE IF( NB.EQ.1 ) THEN
  750. * ------------------------------
  751. * version 1: back-transform each vector with GEMV, Q*x.
  752. IF( KI.GT.2 ) THEN
  753. CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
  754. $ WORK( 1 + (IV-1)*N ), 1,
  755. $ WORK( KI-1 + (IV-1)*N ), VR(1,KI-1), 1)
  756. CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
  757. $ WORK( 1 + (IV)*N ), 1,
  758. $ WORK( KI + (IV)*N ), VR( 1, KI ), 1 )
  759. ELSE
  760. CALL DSCAL( N, WORK(KI-1+(IV-1)*N), VR(1,KI-1), 1)
  761. CALL DSCAL( N, WORK(KI +(IV )*N), VR(1,KI ), 1)
  762. END IF
  763. *
  764. EMAX = ZERO
  765. DO 120 K = 1, N
  766. EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
  767. $ ABS( VR( K, KI ) ) )
  768. 120 CONTINUE
  769. REMAX = ONE / EMAX
  770. CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
  771. CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
  772. *
  773. ELSE
  774. * ------------------------------
  775. * version 2: back-transform block of vectors with GEMM
  776. * zero out below vector
  777. DO K = KI + 1, N
  778. WORK( K + (IV-1)*N ) = ZERO
  779. WORK( K + (IV )*N ) = ZERO
  780. END DO
  781. ISCOMPLEX( IV-1 ) = -IP
  782. ISCOMPLEX( IV ) = IP
  783. IV = IV - 1
  784. * back-transform and normalization is done below
  785. END IF
  786. END IF
  787.  
  788. IF( NB.GT.1 ) THEN
  789. * --------------------------------------------------------
  790. * Blocked version of back-transform
  791. * For complex case, KI2 includes both vectors (KI-1 and KI)
  792. IF( IP.EQ.0 ) THEN
  793. KI2 = KI
  794. ELSE
  795. KI2 = KI - 1
  796. END IF
  797.  
  798. * Columns IV:NB of work are valid vectors.
  799. * When the number of vectors stored reaches NB-1 or NB,
  800. * or if this was last vector, do the GEMM
  801. IF( (IV.LE.2) .OR. (KI2.EQ.1) ) THEN
  802. CALL DGEMM( 'N', 'N', N, NB-IV+1, KI2+NB-IV, ONE,
  803. $ VR, LDVR,
  804. $ WORK( 1 + (IV)*N ), N,
  805. $ ZERO,
  806. $ WORK( 1 + (NB+IV)*N ), N )
  807. * normalize vectors
  808. DO K = IV, NB
  809. IF( ISCOMPLEX(K).EQ.0 ) THEN
  810. * real eigenvector
  811. II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
  812. REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
  813. ELSE IF( ISCOMPLEX(K).EQ.1 ) THEN
  814. * first eigenvector of conjugate pair
  815. EMAX = ZERO
  816. DO II = 1, N
  817. EMAX = MAX( EMAX,
  818. $ ABS( WORK( II + (NB+K )*N ) )+
  819. $ ABS( WORK( II + (NB+K+1)*N ) ) )
  820. END DO
  821. REMAX = ONE / EMAX
  822. * else if ISCOMPLEX(K).EQ.-1
  823. * second eigenvector of conjugate pair
  824. * reuse same REMAX as previous K
  825. END IF
  826. CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
  827. END DO
  828. CALL DLACPY( 'F', N, NB-IV+1,
  829. $ WORK( 1 + (NB+IV)*N ), N,
  830. $ VR( 1, KI2 ), LDVR )
  831. IV = NB
  832. ELSE
  833. IV = IV - 1
  834. END IF
  835. END IF
  836. *! blocked back-transform
  837. *
  838. IS = IS - 1
  839. IF( IP.NE.0 )
  840. $ IS = IS - 1
  841. 140 CONTINUE
  842. END IF
  843.  
  844. IF( LEFTV ) THEN
  845. *
  846. * ============================================================
  847. * Compute left eigenvectors.
  848. *
  849. * IV is index of column in current block.
  850. * For complex left vector, uses IV for real part and IV+1 for complex part.
  851. * Non-blocked version always uses IV=1;
  852. * blocked version starts with IV=1, goes up to NB-1 or NB.
  853. * (Note the "0-th" column is used for 1-norms computed above.)
  854. IV = 1
  855. IP = 0
  856. IS = 1
  857. DO 260 KI = 1, N
  858. IF( IP.EQ.1 ) THEN
  859. * previous iteration (ki-1) was first of conjugate pair,
  860. * so this ki is second of conjugate pair; skip to end of loop
  861. IP = -1
  862. GO TO 260
  863. ELSE IF( KI.EQ.N ) THEN
  864. * last column, so this ki must be real eigenvalue
  865. IP = 0
  866. ELSE IF( T( KI+1, KI ).EQ.ZERO ) THEN
  867. * zero on sub-diagonal, so this ki is real eigenvalue
  868. IP = 0
  869. ELSE
  870. * non-zero on sub-diagonal, so this ki is first of conjugate pair
  871. IP = 1
  872. END IF
  873. *
  874. IF( SOMEV ) THEN
  875. IF( .NOT.SELECT( KI ) )
  876. $ GO TO 260
  877. END IF
  878. *
  879. * Compute the KI-th eigenvalue (WR,WI).
  880. *
  881. WR = T( KI, KI )
  882. WI = ZERO
  883. IF( IP.NE.0 )
  884. $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
  885. $ SQRT( ABS( T( KI+1, KI ) ) )
  886. SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
  887. *
  888. IF( IP.EQ.0 ) THEN
  889. *
  890. * --------------------------------------------------------
  891. * Real left eigenvector
  892. *
  893. WORK( KI + IV*N ) = ONE
  894. *
  895. * Form right-hand side.
  896. *
  897. DO 160 K = KI + 1, N
  898. WORK( K + IV*N ) = -T( KI, K )
  899. 160 CONTINUE
  900. *
  901. * Solve transposed quasi-triangular system:
  902. * [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
  903. *
  904. VMAX = ONE
  905. VCRIT = BIGNUM
  906. *
  907. JNXT = KI + 1
  908. DO 170 J = KI + 1, N
  909. IF( J.LT.JNXT )
  910. $ GO TO 170
  911. J1 = J
  912. J2 = J
  913. JNXT = J + 1
  914. IF( J.LT.N ) THEN
  915. IF( T( J+1, J ).NE.ZERO ) THEN
  916. J2 = J + 1
  917. JNXT = J + 2
  918. END IF
  919. END IF
  920. *
  921. IF( J1.EQ.J2 ) THEN
  922. *
  923. * 1-by-1 diagonal block
  924. *
  925. * Scale if necessary to avoid overflow when forming
  926. * the right-hand side.
  927. *
  928. IF( WORK( J ).GT.VCRIT ) THEN
  929. REC = ONE / VMAX
  930. CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
  931. VMAX = ONE
  932. VCRIT = BIGNUM
  933. END IF
  934. *
  935. WORK( J+IV*N ) = WORK( J+IV*N ) -
  936. $ DDOT( J-KI-1, T( KI+1, J ), 1,
  937. $ WORK( KI+1+IV*N ), 1 )
  938. *
  939. * Solve [ T(J,J) - WR ]**T * X = WORK
  940. *
  941. CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
  942. $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
  943. $ ZERO, X, 2, SCALE, XNORM, IERR )
  944. *
  945. * Scale if necessary
  946. *
  947. IF( SCALE.NE.ONE )
  948. $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
  949. WORK( J+IV*N ) = X( 1, 1 )
  950. VMAX = MAX( ABS( WORK( J+IV*N ) ), VMAX )
  951. VCRIT = BIGNUM / VMAX
  952. *
  953. ELSE
  954. *
  955. * 2-by-2 diagonal block
  956. *
  957. * Scale if necessary to avoid overflow when forming
  958. * the right-hand side.
  959. *
  960. BETA = MAX( WORK( J ), WORK( J+1 ) )
  961. IF( BETA.GT.VCRIT ) THEN
  962. REC = ONE / VMAX
  963. CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
  964. VMAX = ONE
  965. VCRIT = BIGNUM
  966. END IF
  967. *
  968. WORK( J+IV*N ) = WORK( J+IV*N ) -
  969. $ DDOT( J-KI-1, T( KI+1, J ), 1,
  970. $ WORK( KI+1+IV*N ), 1 )
  971. *
  972. WORK( J+1+IV*N ) = WORK( J+1+IV*N ) -
  973. $ DDOT( J-KI-1, T( KI+1, J+1 ), 1,
  974. $ WORK( KI+1+IV*N ), 1 )
  975. *
  976. * Solve
  977. * [ T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
  978. * [ T(J+1,J) T(J+1,J+1)-WR ] ( WORK2 )
  979. *
  980. CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
  981. $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
  982. $ ZERO, X, 2, SCALE, XNORM, IERR )
  983. *
  984. * Scale if necessary
  985. *
  986. IF( SCALE.NE.ONE )
  987. $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
  988. WORK( J +IV*N ) = X( 1, 1 )
  989. WORK( J+1+IV*N ) = X( 2, 1 )
  990. *
  991. VMAX = MAX( ABS( WORK( J +IV*N ) ),
  992. $ ABS( WORK( J+1+IV*N ) ), VMAX )
  993. VCRIT = BIGNUM / VMAX
  994. *
  995. END IF
  996. 170 CONTINUE
  997. *
  998. * Copy the vector x or Q*x to VL and normalize.
  999. *
  1000. IF( .NOT.OVER ) THEN
  1001. * ------------------------------
  1002. * no back-transform: copy x to VL and normalize.
  1003. CALL DCOPY( N-KI+1, WORK( KI + IV*N ), 1,
  1004. $ VL( KI, IS ), 1 )
  1005. *
  1006. II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
  1007. REMAX = ONE / ABS( VL( II, IS ) )
  1008. CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
  1009. *
  1010. DO 180 K = 1, KI - 1
  1011. VL( K, IS ) = ZERO
  1012. 180 CONTINUE
  1013. *
  1014. ELSE IF( NB.EQ.1 ) THEN
  1015. * ------------------------------
  1016. * version 1: back-transform each vector with GEMV, Q*x.
  1017. IF( KI.LT.N )
  1018. $ CALL DGEMV( 'N', N, N-KI, ONE,
  1019. $ VL( 1, KI+1 ), LDVL,
  1020. $ WORK( KI+1 + IV*N ), 1,
  1021. $ WORK( KI + IV*N ), VL( 1, KI ), 1 )
  1022. *
  1023. II = IDAMAX( N, VL( 1, KI ), 1 )
  1024. REMAX = ONE / ABS( VL( II, KI ) )
  1025. CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
  1026. *
  1027. ELSE
  1028. * ------------------------------
  1029. * version 2: back-transform block of vectors with GEMM
  1030. * zero out above vector
  1031. * could go from KI-NV+1 to KI-1
  1032. DO K = 1, KI - 1
  1033. WORK( K + IV*N ) = ZERO
  1034. END DO
  1035. ISCOMPLEX( IV ) = IP
  1036. * back-transform and normalization is done below
  1037. END IF
  1038. ELSE
  1039. *
  1040. * --------------------------------------------------------
  1041. * Complex left eigenvector.
  1042. *
  1043. * Initial solve:
  1044. * [ ( T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI) ]*X = 0.
  1045. * [ ( T(KI+1,KI) T(KI+1,KI+1) ) ]
  1046. *
  1047. IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
  1048. WORK( KI + (IV )*N ) = WI / T( KI, KI+1 )
  1049. WORK( KI+1 + (IV+1)*N ) = ONE
  1050. ELSE
  1051. WORK( KI + (IV )*N ) = ONE
  1052. WORK( KI+1 + (IV+1)*N ) = -WI / T( KI+1, KI )
  1053. END IF
  1054. WORK( KI+1 + (IV )*N ) = ZERO
  1055. WORK( KI + (IV+1)*N ) = ZERO
  1056. *
  1057. * Form right-hand side.
  1058. *
  1059. DO 190 K = KI + 2, N
  1060. WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(KI, K)
  1061. WORK( K+(IV+1)*N ) = -WORK( KI+1+(IV+1)*N )*T(KI+1,K)
  1062. 190 CONTINUE
  1063. *
  1064. * Solve transposed quasi-triangular system:
  1065. * [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
  1066. *
  1067. VMAX = ONE
  1068. VCRIT = BIGNUM
  1069. *
  1070. JNXT = KI + 2
  1071. DO 200 J = KI + 2, N
  1072. IF( J.LT.JNXT )
  1073. $ GO TO 200
  1074. J1 = J
  1075. J2 = J
  1076. JNXT = J + 1
  1077. IF( J.LT.N ) THEN
  1078. IF( T( J+1, J ).NE.ZERO ) THEN
  1079. J2 = J + 1
  1080. JNXT = J + 2
  1081. END IF
  1082. END IF
  1083. *
  1084. IF( J1.EQ.J2 ) THEN
  1085. *
  1086. * 1-by-1 diagonal block
  1087. *
  1088. * Scale if necessary to avoid overflow when
  1089. * forming the right-hand side elements.
  1090. *
  1091. IF( WORK( J ).GT.VCRIT ) THEN
  1092. REC = ONE / VMAX
  1093. CALL DSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
  1094. CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
  1095. VMAX = ONE
  1096. VCRIT = BIGNUM
  1097. END IF
  1098. *
  1099. WORK( J+(IV )*N ) = WORK( J+(IV)*N ) -
  1100. $ DDOT( J-KI-2, T( KI+2, J ), 1,
  1101. $ WORK( KI+2+(IV)*N ), 1 )
  1102. WORK( J+(IV+1)*N ) = WORK( J+(IV+1)*N ) -
  1103. $ DDOT( J-KI-2, T( KI+2, J ), 1,
  1104. $ WORK( KI+2+(IV+1)*N ), 1 )
  1105. *
  1106. * Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
  1107. *
  1108. CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
  1109. $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
  1110. $ -WI, X, 2, SCALE, XNORM, IERR )
  1111. *
  1112. * Scale if necessary
  1113. *
  1114. IF( SCALE.NE.ONE ) THEN
  1115. CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
  1116. CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
  1117. END IF
  1118. WORK( J+(IV )*N ) = X( 1, 1 )
  1119. WORK( J+(IV+1)*N ) = X( 1, 2 )
  1120. VMAX = MAX( ABS( WORK( J+(IV )*N ) ),
  1121. $ ABS( WORK( J+(IV+1)*N ) ), VMAX )
  1122. VCRIT = BIGNUM / VMAX
  1123. *
  1124. ELSE
  1125. *
  1126. * 2-by-2 diagonal block
  1127. *
  1128. * Scale if necessary to avoid overflow when forming
  1129. * the right-hand side elements.
  1130. *
  1131. BETA = MAX( WORK( J ), WORK( J+1 ) )
  1132. IF( BETA.GT.VCRIT ) THEN
  1133. REC = ONE / VMAX
  1134. CALL DSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
  1135. CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
  1136. VMAX = ONE
  1137. VCRIT = BIGNUM
  1138. END IF
  1139. *
  1140. WORK( J +(IV )*N ) = WORK( J+(IV)*N ) -
  1141. $ DDOT( J-KI-2, T( KI+2, J ), 1,
  1142. $ WORK( KI+2+(IV)*N ), 1 )
  1143. *
  1144. WORK( J +(IV+1)*N ) = WORK( J+(IV+1)*N ) -
  1145. $ DDOT( J-KI-2, T( KI+2, J ), 1,
  1146. $ WORK( KI+2+(IV+1)*N ), 1 )
  1147. *
  1148. WORK( J+1+(IV )*N ) = WORK( J+1+(IV)*N ) -
  1149. $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
  1150. $ WORK( KI+2+(IV)*N ), 1 )
  1151. *
  1152. WORK( J+1+(IV+1)*N ) = WORK( J+1+(IV+1)*N ) -
  1153. $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
  1154. $ WORK( KI+2+(IV+1)*N ), 1 )
  1155. *
  1156. * Solve 2-by-2 complex linear equation
  1157. * [ (T(j,j) T(j,j+1) )**T - (wr-i*wi)*I ]*X = SCALE*B
  1158. * [ (T(j+1,j) T(j+1,j+1)) ]
  1159. *
  1160. CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
  1161. $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
  1162. $ -WI, X, 2, SCALE, XNORM, IERR )
  1163. *
  1164. * Scale if necessary
  1165. *
  1166. IF( SCALE.NE.ONE ) THEN
  1167. CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
  1168. CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
  1169. END IF
  1170. WORK( J +(IV )*N ) = X( 1, 1 )
  1171. WORK( J +(IV+1)*N ) = X( 1, 2 )
  1172. WORK( J+1+(IV )*N ) = X( 2, 1 )
  1173. WORK( J+1+(IV+1)*N ) = X( 2, 2 )
  1174. VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
  1175. $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ),
  1176. $ VMAX )
  1177. VCRIT = BIGNUM / VMAX
  1178. *
  1179. END IF
  1180. 200 CONTINUE
  1181. *
  1182. * Copy the vector x or Q*x to VL and normalize.
  1183. *
  1184. IF( .NOT.OVER ) THEN
  1185. * ------------------------------
  1186. * no back-transform: copy x to VL and normalize.
  1187. CALL DCOPY( N-KI+1, WORK( KI + (IV )*N ), 1,
  1188. $ VL( KI, IS ), 1 )
  1189. CALL DCOPY( N-KI+1, WORK( KI + (IV+1)*N ), 1,
  1190. $ VL( KI, IS+1 ), 1 )
  1191. *
  1192. EMAX = ZERO
  1193. DO 220 K = KI, N
  1194. EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
  1195. $ ABS( VL( K, IS+1 ) ) )
  1196. 220 CONTINUE
  1197. REMAX = ONE / EMAX
  1198. CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
  1199. CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
  1200. *
  1201. DO 230 K = 1, KI - 1
  1202. VL( K, IS ) = ZERO
  1203. VL( K, IS+1 ) = ZERO
  1204. 230 CONTINUE
  1205. *
  1206. ELSE IF( NB.EQ.1 ) THEN
  1207. * ------------------------------
  1208. * version 1: back-transform each vector with GEMV, Q*x.
  1209. IF( KI.LT.N-1 ) THEN
  1210. CALL DGEMV( 'N', N, N-KI-1, ONE,
  1211. $ VL( 1, KI+2 ), LDVL,
  1212. $ WORK( KI+2 + (IV)*N ), 1,
  1213. $ WORK( KI + (IV)*N ),
  1214. $ VL( 1, KI ), 1 )
  1215. CALL DGEMV( 'N', N, N-KI-1, ONE,
  1216. $ VL( 1, KI+2 ), LDVL,
  1217. $ WORK( KI+2 + (IV+1)*N ), 1,
  1218. $ WORK( KI+1 + (IV+1)*N ),
  1219. $ VL( 1, KI+1 ), 1 )
  1220. ELSE
  1221. CALL DSCAL( N, WORK(KI+ (IV )*N), VL(1, KI ), 1)
  1222. CALL DSCAL( N, WORK(KI+1+(IV+1)*N), VL(1, KI+1), 1)
  1223. END IF
  1224. *
  1225. EMAX = ZERO
  1226. DO 240 K = 1, N
  1227. EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
  1228. $ ABS( VL( K, KI+1 ) ) )
  1229. 240 CONTINUE
  1230. REMAX = ONE / EMAX
  1231. CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
  1232. CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
  1233. *
  1234. ELSE
  1235. * ------------------------------
  1236. * version 2: back-transform block of vectors with GEMM
  1237. * zero out above vector
  1238. * could go from KI-NV+1 to KI-1
  1239. DO K = 1, KI - 1
  1240. WORK( K + (IV )*N ) = ZERO
  1241. WORK( K + (IV+1)*N ) = ZERO
  1242. END DO
  1243. ISCOMPLEX( IV ) = IP
  1244. ISCOMPLEX( IV+1 ) = -IP
  1245. IV = IV + 1
  1246. * back-transform and normalization is done below
  1247. END IF
  1248. END IF
  1249.  
  1250. IF( NB.GT.1 ) THEN
  1251. * --------------------------------------------------------
  1252. * Blocked version of back-transform
  1253. * For complex case, KI2 includes both vectors (KI and KI+1)
  1254. IF( IP.EQ.0 ) THEN
  1255. KI2 = KI
  1256. ELSE
  1257. KI2 = KI + 1
  1258. END IF
  1259.  
  1260. * Columns 1:IV of work are valid vectors.
  1261. * When the number of vectors stored reaches NB-1 or NB,
  1262. * or if this was last vector, do the GEMM
  1263. IF( (IV.GE.NB-1) .OR. (KI2.EQ.N) ) THEN
  1264. CALL DGEMM( 'N', 'N', N, IV, N-KI2+IV, ONE,
  1265. $ VL( 1, KI2-IV+1 ), LDVL,
  1266. $ WORK( KI2-IV+1 + (1)*N ), N,
  1267. $ ZERO,
  1268. $ WORK( 1 + (NB+1)*N ), N )
  1269. * normalize vectors
  1270. DO K = 1, IV
  1271. IF( ISCOMPLEX(K).EQ.0) THEN
  1272. * real eigenvector
  1273. II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
  1274. REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
  1275. ELSE IF( ISCOMPLEX(K).EQ.1) THEN
  1276. * first eigenvector of conjugate pair
  1277. EMAX = ZERO
  1278. DO II = 1, N
  1279. EMAX = MAX( EMAX,
  1280. $ ABS( WORK( II + (NB+K )*N ) )+
  1281. $ ABS( WORK( II + (NB+K+1)*N ) ) )
  1282. END DO
  1283. REMAX = ONE / EMAX
  1284. * else if ISCOMPLEX(K).EQ.-1
  1285. * second eigenvector of conjugate pair
  1286. * reuse same REMAX as previous K
  1287. END IF
  1288. CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
  1289. END DO
  1290. CALL DLACPY( 'F', N, IV,
  1291. $ WORK( 1 + (NB+1)*N ), N,
  1292. $ VL( 1, KI2-IV+1 ), LDVL )
  1293. IV = 1
  1294. ELSE
  1295. IV = IV + 1
  1296. END IF
  1297. END IF
  1298. *! blocked back-transform
  1299. *
  1300. IS = IS + 1
  1301. IF( IP.NE.0 )
  1302. $ IS = IS + 1
  1303. 260 CONTINUE
  1304. END IF
  1305. *
  1306. RETURN
  1307. *
  1308. * End of DTREVC3
  1309. *
  1310. END
  1311.  
  1312.  
  1313.  

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