Télécharger dtrsen.eso

Retour à la liste

Numérotation des lignes :

  1. C DTRSEN SOURCE BP208322 15/10/13 21:16:02 8670
  2. *> \brief \b DTRSEN
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DTRSEN + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsen.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsen.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsen.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
  23. * M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
  24. *
  25. * .. Scalar Arguments ..
  26. * CHARACTER COMPQ, JOB
  27. * INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N
  28. * REAL*8 S, SEP
  29. * ..
  30. * .. Array Arguments ..
  31. * LOGICAL SELECT( * )
  32. * INTEGER IWORK( * )
  33. * REAL*8 Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
  34. * $ WR( * )
  35. * ..
  36. *
  37. *
  38. *> \par Purpose:
  39. * =============
  40. *>
  41. *> \verbatim
  42. *>
  43. *> DTRSEN reorders the real Schur factorization of a real matrix
  44. *> A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in
  45. *> the leading diagonal blocks of the upper quasi-triangular matrix T,
  46. *> and the leading columns of Q form an orthonormal basis of the
  47. *> corresponding right invariant subspace.
  48. *>
  49. *> Optionally the routine computes the reciprocal condition numbers of
  50. *> the cluster of eigenvalues and/or the invariant subspace.
  51. *>
  52. *> T must be in Schur canonical form (as returned by DHSEQR), that is,
  53. *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
  54. *> 2-by-2 diagonal block has its diagonal elements equal and its
  55. *> off-diagonal elements of opposite sign.
  56. *> \endverbatim
  57. *
  58. * Arguments:
  59. * ==========
  60. *
  61. *> \param[in] JOB
  62. *> \verbatim
  63. *> JOB is CHARACTER*1
  64. *> Specifies whether condition numbers are required for the
  65. *> cluster of eigenvalues (S) or the invariant subspace (SEP):
  66. *> = 'N': none;
  67. *> = 'E': for eigenvalues only (S);
  68. *> = 'V': for invariant subspace only (SEP);
  69. *> = 'B': for both eigenvalues and invariant subspace (S and
  70. *> SEP).
  71. *> \endverbatim
  72. *>
  73. *> \param[in] COMPQ
  74. *> \verbatim
  75. *> COMPQ is CHARACTER*1
  76. *> = 'V': update the matrix Q of Schur vectors;
  77. *> = 'N': do not update Q.
  78. *> \endverbatim
  79. *>
  80. *> \param[in] SELECT
  81. *> \verbatim
  82. *> SELECT is LOGICAL array, dimension (N)
  83. *> SELECT specifies the eigenvalues in the selected cluster. To
  84. *> select a real eigenvalue w(j), SELECT(j) must be set to
  85. *> .TRUE.. To select a complex conjugate pair of eigenvalues
  86. *> w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
  87. *> either SELECT(j) or SELECT(j+1) or both must be set to
  88. *> .TRUE.; a complex conjugate pair of eigenvalues must be
  89. *> either both included in the cluster or both excluded.
  90. *> \endverbatim
  91. *>
  92. *> \param[in] N
  93. *> \verbatim
  94. *> N is INTEGER
  95. *> The order of the matrix T. N >= 0.
  96. *> \endverbatim
  97. *>
  98. *> \param[in,out] T
  99. *> \verbatim
  100. *> T is DOUBLE PRECISION array, dimension (LDT,N)
  101. *> On entry, the upper quasi-triangular matrix T, in Schur
  102. *> canonical form.
  103. *> On exit, T is overwritten by the reordered matrix T, again in
  104. *> Schur canonical form, with the selected eigenvalues in the
  105. *> leading diagonal blocks.
  106. *> \endverbatim
  107. *>
  108. *> \param[in] LDT
  109. *> \verbatim
  110. *> LDT is INTEGER
  111. *> The leading dimension of the array T. LDT >= max(1,N).
  112. *> \endverbatim
  113. *>
  114. *> \param[in,out] Q
  115. *> \verbatim
  116. *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
  117. *> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
  118. *> On exit, if COMPQ = 'V', Q has been postmultiplied by the
  119. *> orthogonal transformation matrix which reorders T; the
  120. *> leading M columns of Q form an orthonormal basis for the
  121. *> specified invariant subspace.
  122. *> If COMPQ = 'N', Q is not referenced.
  123. *> \endverbatim
  124. *>
  125. *> \param[in] LDQ
  126. *> \verbatim
  127. *> LDQ is INTEGER
  128. *> The leading dimension of the array Q.
  129. *> LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
  130. *> \endverbatim
  131. *>
  132. *> \param[out] WR
  133. *> \verbatim
  134. *> WR is DOUBLE PRECISION array, dimension (N)
  135. *> \endverbatim
  136. *> \param[out] WI
  137. *> \verbatim
  138. *> WI is DOUBLE PRECISION array, dimension (N)
  139. *>
  140. *> The real and imaginary parts, respectively, of the reordered
  141. *> eigenvalues of T. The eigenvalues are stored in the same
  142. *> order as on the diagonal of T, with WR(i) = T(i,i) and, if
  143. *> T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and
  144. *> WI(i+1) = -WI(i). Note that if a complex eigenvalue is
  145. *> sufficiently ill-conditioned, then its value may differ
  146. *> significantly from its value before reordering.
  147. *> \endverbatim
  148. *>
  149. *> \param[out] M
  150. *> \verbatim
  151. *> M is INTEGER
  152. *> The dimension of the specified invariant subspace.
  153. *> 0 < = M <= N.
  154. *> \endverbatim
  155. *>
  156. *> \param[out] S
  157. *> \verbatim
  158. *> S is DOUBLE PRECISION
  159. *> If JOB = 'E' or 'B', S is a lower bound on the reciprocal
  160. *> condition number for the selected cluster of eigenvalues.
  161. *> S cannot underestimate the true reciprocal condition number
  162. *> by more than a factor of sqrt(N). If M = 0 or N, S = 1.
  163. *> If JOB = 'N' or 'V', S is not referenced.
  164. *> \endverbatim
  165. *>
  166. *> \param[out] SEP
  167. *> \verbatim
  168. *> SEP is DOUBLE PRECISION
  169. *> If JOB = 'V' or 'B', SEP is the estimated reciprocal
  170. *> condition number of the specified invariant subspace. If
  171. *> M = 0 or N, SEP = norm(T).
  172. *> If JOB = 'N' or 'E', SEP is not referenced.
  173. *> \endverbatim
  174. *>
  175. *> \param[out] WORK
  176. *> \verbatim
  177. *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  178. *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  179. *> \endverbatim
  180. *>
  181. *> \param[in] LWORK
  182. *> \verbatim
  183. *> LWORK is INTEGER
  184. *> The dimension of the array WORK.
  185. *> If JOB = 'N', LWORK >= max(1,N);
  186. *> if JOB = 'E', LWORK >= max(1,M*(N-M));
  187. *> if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
  188. *>
  189. *> If LWORK = -1, then a workspace query is assumed; the routine
  190. *> only calculates the optimal size of the WORK array, returns
  191. *> this value as the first entry of the WORK array, and no error
  192. *> message related to LWORK is issued by XERBLA.
  193. *> \endverbatim
  194. *>
  195. *> \param[out] IWORK
  196. *> \verbatim
  197. *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
  198. *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
  199. *> \endverbatim
  200. *>
  201. *> \param[in] LIWORK
  202. *> \verbatim
  203. *> LIWORK is INTEGER
  204. *> The dimension of the array IWORK.
  205. *> If JOB = 'N' or 'E', LIWORK >= 1;
  206. *> if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).
  207. *>
  208. *> If LIWORK = -1, then a workspace query is assumed; the
  209. *> routine only calculates the optimal size of the IWORK array,
  210. *> returns this value as the first entry of the IWORK array, and
  211. *> no error message related to LIWORK is issued by XERBLA.
  212. *> \endverbatim
  213. *>
  214. *> \param[out] INFO
  215. *> \verbatim
  216. *> INFO is INTEGER
  217. *> = 0: successful exit
  218. *> < 0: if INFO = -i, the i-th argument had an illegal value
  219. *> = 1: reordering of T failed because some eigenvalues are too
  220. *> close to separate (the problem is very ill-conditioned);
  221. *> T may have been partially reordered, and WR and WI
  222. *> contain the eigenvalues in the same order as in T; S and
  223. *> SEP (if requested) are set to zero.
  224. *> \endverbatim
  225. *
  226. * Authors:
  227. * ========
  228. *
  229. *> \author Univ. of Tennessee
  230. *> \author Univ. of California Berkeley
  231. *> \author Univ. of Colorado Denver
  232. *> \author NAG Ltd.
  233. *
  234. *> \date April 2012
  235. *
  236. *> \ingroup doubleOTHERcomputational
  237. *
  238. *> \par Further Details:
  239. * =====================
  240. *>
  241. *> \verbatim
  242. *>
  243. *> DTRSEN first collects the selected eigenvalues by computing an
  244. *> orthogonal transformation Z to move them to the top left corner of T.
  245. *> In other words, the selected eigenvalues are the eigenvalues of T11
  246. *> in:
  247. *>
  248. *> Z**T * T * Z = ( T11 T12 ) n1
  249. *> ( 0 T22 ) n2
  250. *> n1 n2
  251. *>
  252. *> where N = n1+n2 and Z**T means the transpose of Z. The first n1 columns
  253. *> of Z span the specified invariant subspace of T.
  254. *>
  255. *> If T has been obtained from the real Schur factorization of a matrix
  256. *> A = Q*T*Q**T, then the reordered real Schur factorization of A is given
  257. *> by A = (Q*Z)*(Z**T*T*Z)*(Q*Z)**T, and the first n1 columns of Q*Z span
  258. *> the corresponding invariant subspace of A.
  259. *>
  260. *> The reciprocal condition number of the average of the eigenvalues of
  261. *> T11 may be returned in S. S lies between 0 (very badly conditioned)
  262. *> and 1 (very well conditioned). It is computed as follows. First we
  263. *> compute R so that
  264. *>
  265. *> P = ( I R ) n1
  266. *> ( 0 0 ) n2
  267. *> n1 n2
  268. *>
  269. *> is the projector on the invariant subspace associated with T11.
  270. *> R is the solution of the Sylvester equation:
  271. *>
  272. *> T11*R - R*T22 = T12.
  273. *>
  274. *> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
  275. *> the two-norm of M. Then S is computed as the lower bound
  276. *>
  277. *> (1 + F-norm(R)**2)**(-1/2)
  278. *>
  279. *> on the reciprocal of 2-norm(P), the true reciprocal condition number.
  280. *> S cannot underestimate 1 / 2-norm(P) by more than a factor of
  281. *> sqrt(N).
  282. *>
  283. *> An approximate error bound for the computed average of the
  284. *> eigenvalues of T11 is
  285. *>
  286. *> EPS * norm(T) / S
  287. *>
  288. *> where EPS is the machine precision.
  289. *>
  290. *> The reciprocal condition number of the right invariant subspace
  291. *> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
  292. *> SEP is defined as the separation of T11 and T22:
  293. *>
  294. *> sep( T11, T22 ) = sigma-min( C )
  295. *>
  296. *> where sigma-min(C) is the smallest singular value of the
  297. *> n1*n2-by-n1*n2 matrix
  298. *>
  299. *> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
  300. *>
  301. *> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
  302. *> product. We estimate sigma-min(C) by the reciprocal of an estimate of
  303. *> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
  304. *> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
  305. *>
  306. *> When SEP is small, small changes in T can cause large changes in
  307. *> the invariant subspace. An approximate bound on the maximum angular
  308. *> error in the computed right invariant subspace is
  309. *>
  310. *> EPS * norm(T) / SEP
  311. *> \endverbatim
  312. *>
  313. * =====================================================================
  314. SUBROUTINE DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
  315. $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
  316. *
  317. * -- LAPACK computational routine (version 3.4.1) --
  318. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  319. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  320. * April 2012
  321. *
  322. * .. Scalar Arguments ..
  323. CHARACTER COMPQ, JOB
  324. INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N
  325. REAL*8 S, SEP
  326. * ..
  327. * .. Array Arguments ..
  328. LOGICAL SELECT( * )
  329. INTEGER IWORK( * )
  330. REAL*8 Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
  331. $ WR( * )
  332. * ..
  333. *
  334. * =====================================================================
  335. *
  336. * .. Parameters ..
  337. REAL*8 ZERO, ONE
  338. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  339. * ..
  340. * .. Local Scalars ..
  341. LOGICAL LQUERY, PAIR, SWAP, WANTBH, WANTQ, WANTS,
  342. $ WANTSP
  343. INTEGER IERR, K, KASE, KK, KS, LIWMIN, LWMIN, N1, N2,
  344. $ NN
  345. REAL*8 EST, RNORM, SCALE
  346. * ..
  347. * .. Local Arrays ..
  348. INTEGER ISAVE( 3 )
  349. * ..
  350. * .. External Functions ..
  351. LOGICAL LSAME
  352. REAL*8 DLANGE
  353. EXTERNAL LSAME, DLANGE
  354. * ..
  355. * .. External Subroutines ..
  356. * ..
  357. ** .. Intrinsic Functions ..
  358. * INTRINSIC ABS, MAX, SQRT
  359. ** ..
  360. ** .. Executable Statements ..
  361. *
  362. * Decode and test the input parameters
  363. *
  364. WANTBH = LSAME( JOB, 'B' )
  365. WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
  366. WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
  367. WANTQ = LSAME( COMPQ, 'V' )
  368. *
  369. INFO = 0
  370. LQUERY = ( LWORK.EQ.-1 )
  371. IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
  372. $ THEN
  373. INFO = -1
  374. ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
  375. INFO = -2
  376. ELSE IF( N.LT.0 ) THEN
  377. INFO = -4
  378. ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  379. INFO = -6
  380. ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
  381. INFO = -8
  382. ELSE
  383. *
  384. * Set M to the dimension of the specified invariant subspace,
  385. * and test LWORK and LIWORK.
  386. *
  387. M = 0
  388. PAIR = .FALSE.
  389. DO 10 K = 1, N
  390. IF( PAIR ) THEN
  391. PAIR = .FALSE.
  392. ELSE
  393. IF( K.LT.N ) THEN
  394. IF( T( K+1, K ).EQ.ZERO ) THEN
  395. IF( SELECT( K ) )
  396. $ M = M + 1
  397. ELSE
  398. PAIR = .TRUE.
  399. IF( SELECT( K ) .OR. SELECT( K+1 ) )
  400. $ M = M + 2
  401. END IF
  402. ELSE
  403. IF( SELECT( N ) )
  404. $ M = M + 1
  405. END IF
  406. END IF
  407. 10 CONTINUE
  408. *
  409. N1 = M
  410. N2 = N - M
  411. NN = N1*N2
  412. *
  413. IF( WANTSP ) THEN
  414. LWMIN = MAX( 1, 2*NN )
  415. LIWMIN = MAX( 1, NN )
  416. ELSE IF( LSAME( JOB, 'N' ) ) THEN
  417. LWMIN = MAX( 1, N )
  418. LIWMIN = 1
  419. ELSE IF( LSAME( JOB, 'E' ) ) THEN
  420. LWMIN = MAX( 1, NN )
  421. LIWMIN = 1
  422. END IF
  423. *
  424. IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  425. INFO = -15
  426. ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
  427. INFO = -17
  428. END IF
  429. END IF
  430. *
  431. IF( INFO.EQ.0 ) THEN
  432. WORK( 1 ) = LWMIN
  433. IWORK( 1 ) = LIWMIN
  434. END IF
  435. *
  436. IF( INFO.NE.0 ) THEN
  437. CALL XERBLA( 'DTRSEN', -INFO )
  438. RETURN
  439. ELSE IF( LQUERY ) THEN
  440. RETURN
  441. END IF
  442. *
  443. * Quick return if possible.
  444. *
  445. IF( M.EQ.N .OR. M.EQ.0 ) THEN
  446. IF( WANTS )
  447. $ S = ONE
  448. IF( WANTSP )
  449. $ SEP = DLANGE( '1', N, N, T, LDT, WORK )
  450. GO TO 40
  451. END IF
  452. *
  453. * Collect the selected blocks at the top-left corner of T.
  454. *
  455. KS = 0
  456. PAIR = .FALSE.
  457. DO 20 K = 1, N
  458. IF( PAIR ) THEN
  459. PAIR = .FALSE.
  460. ELSE
  461. SWAP = SELECT( K )
  462. IF( K.LT.N ) THEN
  463. IF( T( K+1, K ).NE.ZERO ) THEN
  464. PAIR = .TRUE.
  465. SWAP = SWAP .OR. SELECT( K+1 )
  466. END IF
  467. END IF
  468. IF( SWAP ) THEN
  469. KS = KS + 1
  470. *
  471. * Swap the K-th block to position KS.
  472. *
  473. IERR = 0
  474. KK = K
  475. IF( K.NE.KS )
  476. $ CALL DTREXC( COMPQ, N, T, LDT, Q, LDQ, KK, KS, WORK,
  477. $ IERR )
  478. IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
  479. *
  480. * Blocks too close to swap: exit.
  481. *
  482. INFO = 1
  483. IF( WANTS )
  484. $ S = ZERO
  485. IF( WANTSP )
  486. $ SEP = ZERO
  487. GO TO 40
  488. END IF
  489. IF( PAIR )
  490. $ KS = KS + 1
  491. END IF
  492. END IF
  493. 20 CONTINUE
  494. *
  495. IF( WANTS ) THEN
  496. *
  497. * Solve Sylvester equation for R:
  498. *
  499. * T11*R - R*T22 = scale*T12
  500. *
  501. CALL DLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
  502. CALL DTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
  503. $ LDT, WORK, N1, SCALE, IERR )
  504. *
  505. * Estimate the reciprocal of the condition number of the cluster
  506. * of eigenvalues.
  507. *
  508. RNORM = DLANGE( 'F', N1, N2, WORK, N1, WORK )
  509. IF( RNORM.EQ.ZERO ) THEN
  510. S = ONE
  511. ELSE
  512. S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
  513. $ SQRT( RNORM ) )
  514. END IF
  515. END IF
  516. *
  517. IF( WANTSP ) THEN
  518. *
  519. * Estimate sep(T11,T22).
  520. *
  521. EST = ZERO
  522. KASE = 0
  523. 30 CONTINUE
  524. CALL DLACN2( NN, WORK( NN+1 ), WORK, IWORK, EST, KASE, ISAVE )
  525. IF( KASE.NE.0 ) THEN
  526. IF( KASE.EQ.1 ) THEN
  527. *
  528. * Solve T11*R - R*T22 = scale*X.
  529. *
  530. CALL DTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
  531. $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
  532. $ IERR )
  533. ELSE
  534. *
  535. * Solve T11**T*R - R*T22**T = scale*X.
  536. *
  537. CALL DTRSYL( 'T', 'T', -1, N1, N2, T, LDT,
  538. $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
  539. $ IERR )
  540. END IF
  541. GO TO 30
  542. END IF
  543. *
  544. SEP = SCALE / EST
  545. END IF
  546. *
  547. 40 CONTINUE
  548. *
  549. * Store the output eigenvalues in WR and WI.
  550. *
  551. DO 50 K = 1, N
  552. WR( K ) = T( K, K )
  553. WI( K ) = ZERO
  554. 50 CONTINUE
  555. DO 60 K = 1, N - 1
  556. IF( T( K+1, K ).NE.ZERO ) THEN
  557. WI( K ) = SQRT( ABS( T( K, K+1 ) ) )*
  558. $ SQRT( ABS( T( K+1, K ) ) )
  559. WI( K+1 ) = -WI( K )
  560. END IF
  561. 60 CONTINUE
  562. *
  563. WORK( 1 ) = LWMIN
  564. IWORK( 1 ) = LIWMIN
  565. *
  566. RETURN
  567. *
  568. * End of DTRSEN
  569. *
  570. END
  571.  
  572.  

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