Télécharger dtrexc.eso

Retour à la liste

Numérotation des lignes :

  1. C DTREXC SOURCE BP208322 15/10/13 21:16:00 8670
  2. *> \brief \b DTREXC
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DTREXC + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrexc.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrexc.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrexc.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
  23. * INFO )
  24. *
  25. * .. Scalar Arguments ..
  26. * CHARACTER COMPQ
  27. * INTEGER IFST, ILST, INFO, LDQ, LDT, N
  28. * ..
  29. * .. Array Arguments ..
  30. * REAL*8 Q( LDQ, * ), T( LDT, * ), WORK( * )
  31. * ..
  32. *
  33. *
  34. *> \par Purpose:
  35. * =============
  36. *>
  37. *> \verbatim
  38. *>
  39. *> DTREXC reorders the real Schur factorization of a real matrix
  40. *> A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
  41. *> moved to row ILST.
  42. *>
  43. *> The real Schur form T is reordered by an orthogonal similarity
  44. *> transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
  45. *> is updated by postmultiplying it with Z.
  46. *>
  47. *> T must be in Schur canonical form (as returned by DHSEQR), that is,
  48. *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
  49. *> 2-by-2 diagonal block has its diagonal elements equal and its
  50. *> off-diagonal elements of opposite sign.
  51. *> \endverbatim
  52. *
  53. * Arguments:
  54. * ==========
  55. *
  56. *> \param[in] COMPQ
  57. *> \verbatim
  58. *> COMPQ is CHARACTER*1
  59. *> = 'V': update the matrix Q of Schur vectors;
  60. *> = 'N': do not update Q.
  61. *> \endverbatim
  62. *>
  63. *> \param[in] N
  64. *> \verbatim
  65. *> N is INTEGER
  66. *> The order of the matrix T. N >= 0.
  67. *> \endverbatim
  68. *>
  69. *> \param[in,out] T
  70. *> \verbatim
  71. *> T is DOUBLE PRECISION array, dimension (LDT,N)
  72. *> On entry, the upper quasi-triangular matrix T, in Schur
  73. *> Schur canonical form.
  74. *> On exit, the reordered upper quasi-triangular matrix, again
  75. *> in Schur canonical form.
  76. *> \endverbatim
  77. *>
  78. *> \param[in] LDT
  79. *> \verbatim
  80. *> LDT is INTEGER
  81. *> The leading dimension of the array T. LDT >= max(1,N).
  82. *> \endverbatim
  83. *>
  84. *> \param[in,out] Q
  85. *> \verbatim
  86. *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
  87. *> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
  88. *> On exit, if COMPQ = 'V', Q has been postmultiplied by the
  89. *> orthogonal transformation matrix Z which reorders T.
  90. *> If COMPQ = 'N', Q is not referenced.
  91. *> \endverbatim
  92. *>
  93. *> \param[in] LDQ
  94. *> \verbatim
  95. *> LDQ is INTEGER
  96. *> The leading dimension of the array Q. LDQ >= max(1,N).
  97. *> \endverbatim
  98. *>
  99. *> \param[in,out] IFST
  100. *> \verbatim
  101. *> IFST is INTEGER
  102. *> \endverbatim
  103. *>
  104. *> \param[in,out] ILST
  105. *> \verbatim
  106. *> ILST is INTEGER
  107. *>
  108. *> Specify the reordering of the diagonal blocks of T.
  109. *> The block with row index IFST is moved to row ILST, by a
  110. *> sequence of transpositions between adjacent blocks.
  111. *> On exit, if IFST pointed on entry to the second row of a
  112. *> 2-by-2 block, it is changed to point to the first row; ILST
  113. *> always points to the first row of the block in its final
  114. *> position (which may differ from its input value by +1 or -1).
  115. *> 1 <= IFST <= N; 1 <= ILST <= N.
  116. *> \endverbatim
  117. *>
  118. *> \param[out] WORK
  119. *> \verbatim
  120. *> WORK is DOUBLE PRECISION array, dimension (N)
  121. *> \endverbatim
  122. *>
  123. *> \param[out] INFO
  124. *> \verbatim
  125. *> INFO is INTEGER
  126. *> = 0: successful exit
  127. *> < 0: if INFO = -i, the i-th argument had an illegal value
  128. *> = 1: two adjacent blocks were too close to swap (the problem
  129. *> is very ill-conditioned); T may have been partially
  130. *> reordered, and ILST points to the first row of the
  131. *> current position of the block being moved.
  132. *> \endverbatim
  133. *
  134. * Authors:
  135. * ========
  136. *
  137. *> \author Univ. of Tennessee
  138. *> \author Univ. of California Berkeley
  139. *> \author Univ. of Colorado Denver
  140. *> \author NAG Ltd.
  141. *
  142. *> \date November 2011
  143. *
  144. *> \ingroup doubleOTHERcomputational
  145. *
  146. * =====================================================================
  147. SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
  148. $ INFO )
  149. *
  150. * -- LAPACK computational routine (version 3.4.0) --
  151. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  152. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  153. * November 2011
  154. *
  155. * .. Scalar Arguments ..
  156. CHARACTER COMPQ
  157. INTEGER IFST, ILST, INFO, LDQ, LDT, N
  158. * ..
  159. * .. Array Arguments ..
  160. REAL*8 Q( LDQ, * ), T( LDT, * ), WORK( * )
  161. * ..
  162. *
  163. * =====================================================================
  164. *
  165. * .. Parameters ..
  166. REAL*8 ZERO
  167. PARAMETER ( ZERO = 0.0D+0 )
  168. * ..
  169. * .. Local Scalars ..
  170. LOGICAL WANTQ
  171. INTEGER HERE, NBF, NBL, NBNEXT
  172. * ..
  173. * .. External Functions ..
  174. LOGICAL LSAME
  175. EXTERNAL LSAME
  176. * ..
  177. * .. External Subroutines ..
  178. EXTERNAL DLAEXC, XERBLA
  179. * ..
  180. ** .. Intrinsic Functions ..
  181. * INTRINSIC MAX
  182. ** ..
  183. ** .. Executable Statements ..
  184. *
  185. * Decode and test the input arguments.
  186. *
  187. INFO = 0
  188. WANTQ = LSAME( COMPQ, 'V' )
  189. IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
  190. INFO = -1
  191. ELSE IF( N.LT.0 ) THEN
  192. INFO = -2
  193. ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  194. INFO = -4
  195. ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
  196. INFO = -6
  197. ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
  198. INFO = -7
  199. ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
  200. INFO = -8
  201. END IF
  202. IF( INFO.NE.0 ) THEN
  203. CALL XERBLA( 'DTREXC', -INFO )
  204. RETURN
  205. END IF
  206. *
  207. * Quick return if possible
  208. *
  209. IF( N.LE.1 )
  210. $ RETURN
  211. *
  212. * Determine the first row of specified block
  213. * and find out it is 1 by 1 or 2 by 2.
  214. *
  215. IF( IFST.GT.1 ) THEN
  216. IF( T( IFST, IFST-1 ).NE.ZERO )
  217. $ IFST = IFST - 1
  218. END IF
  219. NBF = 1
  220. IF( IFST.LT.N ) THEN
  221. IF( T( IFST+1, IFST ).NE.ZERO )
  222. $ NBF = 2
  223. END IF
  224. *
  225. * Determine the first row of the final block
  226. * and find out it is 1 by 1 or 2 by 2.
  227. *
  228. IF( ILST.GT.1 ) THEN
  229. IF( T( ILST, ILST-1 ).NE.ZERO )
  230. $ ILST = ILST - 1
  231. END IF
  232. NBL = 1
  233. IF( ILST.LT.N ) THEN
  234. IF( T( ILST+1, ILST ).NE.ZERO )
  235. $ NBL = 2
  236. END IF
  237. *
  238. IF( IFST.EQ.ILST )
  239. $ RETURN
  240. *
  241. IF( IFST.LT.ILST ) THEN
  242. *
  243. * Update ILST
  244. *
  245. IF( NBF.EQ.2 .AND. NBL.EQ.1 )
  246. $ ILST = ILST - 1
  247. IF( NBF.EQ.1 .AND. NBL.EQ.2 )
  248. $ ILST = ILST + 1
  249. *
  250. HERE = IFST
  251. *
  252. 10 CONTINUE
  253. *
  254. * Swap block with next one below
  255. *
  256. IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
  257. *
  258. * Current block either 1 by 1 or 2 by 2
  259. *
  260. NBNEXT = 1
  261. IF( HERE+NBF+1.LE.N ) THEN
  262. IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
  263. $ NBNEXT = 2
  264. END IF
  265. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
  266. $ WORK, INFO )
  267. IF( INFO.NE.0 ) THEN
  268. ILST = HERE
  269. RETURN
  270. END IF
  271. HERE = HERE + NBNEXT
  272. *
  273. * Test if 2 by 2 block breaks into two 1 by 1 blocks
  274. *
  275. IF( NBF.EQ.2 ) THEN
  276. IF( T( HERE+1, HERE ).EQ.ZERO )
  277. $ NBF = 3
  278. END IF
  279. *
  280. ELSE
  281. *
  282. * Current block consists of two 1 by 1 blocks each of which
  283. * must be swapped individually
  284. *
  285. NBNEXT = 1
  286. IF( HERE+3.LE.N ) THEN
  287. IF( T( HERE+3, HERE+2 ).NE.ZERO )
  288. $ NBNEXT = 2
  289. END IF
  290. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
  291. $ WORK, INFO )
  292. IF( INFO.NE.0 ) THEN
  293. ILST = HERE
  294. RETURN
  295. END IF
  296. IF( NBNEXT.EQ.1 ) THEN
  297. *
  298. * Swap two 1 by 1 blocks, no problems possible
  299. *
  300. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
  301. $ WORK, INFO )
  302. HERE = HERE + 1
  303. ELSE
  304. *
  305. * Recompute NBNEXT in case 2 by 2 split
  306. *
  307. IF( T( HERE+2, HERE+1 ).EQ.ZERO )
  308. $ NBNEXT = 1
  309. IF( NBNEXT.EQ.2 ) THEN
  310. *
  311. * 2 by 2 Block did not split
  312. *
  313. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
  314. $ NBNEXT, WORK, INFO )
  315. IF( INFO.NE.0 ) THEN
  316. ILST = HERE
  317. RETURN
  318. END IF
  319. HERE = HERE + 2
  320. ELSE
  321. *
  322. * 2 by 2 Block did split
  323. *
  324. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
  325. $ WORK, INFO )
  326. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
  327. $ WORK, INFO )
  328. HERE = HERE + 2
  329. END IF
  330. END IF
  331. END IF
  332. IF( HERE.LT.ILST )
  333. $ GO TO 10
  334. *
  335. ELSE
  336. *
  337. HERE = IFST
  338. 20 CONTINUE
  339. *
  340. * Swap block with next one above
  341. *
  342. IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
  343. *
  344. * Current block either 1 by 1 or 2 by 2
  345. *
  346. NBNEXT = 1
  347. IF( HERE.GE.3 ) THEN
  348. IF( T( HERE-1, HERE-2 ).NE.ZERO )
  349. $ NBNEXT = 2
  350. END IF
  351. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
  352. $ NBF, WORK, INFO )
  353. IF( INFO.NE.0 ) THEN
  354. ILST = HERE
  355. RETURN
  356. END IF
  357. HERE = HERE - NBNEXT
  358. *
  359. * Test if 2 by 2 block breaks into two 1 by 1 blocks
  360. *
  361. IF( NBF.EQ.2 ) THEN
  362. IF( T( HERE+1, HERE ).EQ.ZERO )
  363. $ NBF = 3
  364. END IF
  365. *
  366. ELSE
  367. *
  368. * Current block consists of two 1 by 1 blocks each of which
  369. * must be swapped individually
  370. *
  371. NBNEXT = 1
  372. IF( HERE.GE.3 ) THEN
  373. IF( T( HERE-1, HERE-2 ).NE.ZERO )
  374. $ NBNEXT = 2
  375. END IF
  376. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
  377. $ 1, WORK, INFO )
  378. IF( INFO.NE.0 ) THEN
  379. ILST = HERE
  380. RETURN
  381. END IF
  382. IF( NBNEXT.EQ.1 ) THEN
  383. *
  384. * Swap two 1 by 1 blocks, no problems possible
  385. *
  386. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
  387. $ WORK, INFO )
  388. HERE = HERE - 1
  389. ELSE
  390. *
  391. * Recompute NBNEXT in case 2 by 2 split
  392. *
  393. IF( T( HERE, HERE-1 ).EQ.ZERO )
  394. $ NBNEXT = 1
  395. IF( NBNEXT.EQ.2 ) THEN
  396. *
  397. * 2 by 2 Block did not split
  398. *
  399. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
  400. $ WORK, INFO )
  401. IF( INFO.NE.0 ) THEN
  402. ILST = HERE
  403. RETURN
  404. END IF
  405. HERE = HERE - 2
  406. ELSE
  407. *
  408. * 2 by 2 Block did split
  409. *
  410. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
  411. $ WORK, INFO )
  412. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
  413. $ WORK, INFO )
  414. HERE = HERE - 2
  415. END IF
  416. END IF
  417. END IF
  418. IF( HERE.GT.ILST )
  419. $ GO TO 20
  420. END IF
  421. ILST = HERE
  422. *
  423. RETURN
  424. *
  425. * End of DTREXC
  426. *
  427. END
  428.  
  429.  

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