Télécharger dtrexc.eso

Retour à la liste

Numérotation des lignes :

dtrexc
  1. C DTREXC SOURCE BP208322 18/07/10 21:15:24 9872
  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 >= 1, and if
  97. *> COMPQ = 'V', LDQ >= max(1,N).
  98. *> \endverbatim
  99. *>
  100. *> \param[in,out] IFST
  101. *> \verbatim
  102. *> IFST is INTEGER
  103. *> \endverbatim
  104. *>
  105. *> \param[in,out] ILST
  106. *> \verbatim
  107. *> ILST is INTEGER
  108. *>
  109. *> Specify the reordering of the diagonal blocks of T.
  110. *> The block with row index IFST is moved to row ILST, by a
  111. *> sequence of transpositions between adjacent blocks.
  112. *> On exit, if IFST pointed on entry to the second row of a
  113. *> 2-by-2 block, it is changed to point to the first row; ILST
  114. *> always points to the first row of the block in its final
  115. *> position (which may differ from its input value by +1 or -1).
  116. *> 1 <= IFST <= N; 1 <= ILST <= N.
  117. *> \endverbatim
  118. *>
  119. *> \param[out] WORK
  120. *> \verbatim
  121. *> WORK is DOUBLE PRECISION array, dimension (N)
  122. *> \endverbatim
  123. *>
  124. *> \param[out] INFO
  125. *> \verbatim
  126. *> INFO is INTEGER
  127. *> = 0: successful exit
  128. *> < 0: if INFO = -i, the i-th argument had an illegal value
  129. *> = 1: two adjacent blocks were too close to swap (the problem
  130. *> is very ill-conditioned); T may have been partially
  131. *> reordered, and ILST points to the first row of the
  132. *> current position of the block being moved.
  133. *> \endverbatim
  134. *
  135. * Authors:
  136. * ========
  137. *
  138. *> \author Univ. of Tennessee
  139. *> \author Univ. of California Berkeley
  140. *> \author Univ. of Colorado Denver
  141. *> \author NAG Ltd.
  142. *
  143. *> \date December 2016
  144. *
  145. *> \ingroup doubleOTHERcomputational
  146. *
  147. * =====================================================================
  148. SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
  149. $ INFO )
  150. *
  151. * -- LAPACK computational routine (version 3.7.0) --
  152. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  153. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  154. * December 2016
  155. *
  156. * .. Scalar Arguments ..
  157. CHARACTER COMPQ
  158. INTEGER IFST, ILST, INFO, LDQ, LDT, N
  159. * ..
  160. * .. Array Arguments ..
  161. REAL*8 Q( LDQ, * ), T( LDT, * ), WORK( * )
  162. * ..
  163. *
  164. * =====================================================================
  165. *
  166. * .. Parameters ..
  167. REAL*8 ZERO
  168. PARAMETER ( ZERO = 0.0D+0 )
  169. * ..
  170. * .. Local Scalars ..
  171. LOGICAL WANTQ
  172. INTEGER HERE, NBF, NBL, NBNEXT
  173. * ..
  174. * .. External Functions ..
  175. LOGICAL LSAME
  176. EXTERNAL LSAME
  177. * ..
  178. * .. External Subroutines ..
  179. EXTERNAL DLAEXC, XERBLA
  180. * ..
  181. ** .. Intrinsic Functions ..
  182. * INTRINSIC MAX
  183. ** ..
  184. ** .. Executable Statements ..
  185. *
  186. * Decode and test the input arguments.
  187. *
  188. INFO = 0
  189. WANTQ = LSAME( COMPQ, 'V' )
  190. IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
  191. INFO = -1
  192. ELSE IF( N.LT.0 ) THEN
  193. INFO = -2
  194. ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  195. INFO = -4
  196. ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
  197. INFO = -6
  198. ELSE IF(( IFST.LT.1 .OR. IFST.GT.N ).AND.( N.GT.0 )) THEN
  199. INFO = -7
  200. ELSE IF(( ILST.LT.1 .OR. ILST.GT.N ).AND.( N.GT.0 )) THEN
  201. INFO = -8
  202. END IF
  203. IF( INFO.NE.0 ) THEN
  204. CALL XERBLA( 'DTREXC', -INFO )
  205. RETURN
  206. END IF
  207. *
  208. * Quick return if possible
  209. *
  210. IF( N.LE.1 )
  211. $ RETURN
  212. *
  213. * Determine the first row of specified block
  214. * and find out it is 1 by 1 or 2 by 2.
  215. *
  216. IF( IFST.GT.1 ) THEN
  217. IF( T( IFST, IFST-1 ).NE.ZERO )
  218. $ IFST = IFST - 1
  219. END IF
  220. NBF = 1
  221. IF( IFST.LT.N ) THEN
  222. IF( T( IFST+1, IFST ).NE.ZERO )
  223. $ NBF = 2
  224. END IF
  225. *
  226. * Determine the first row of the final block
  227. * and find out it is 1 by 1 or 2 by 2.
  228. *
  229. IF( ILST.GT.1 ) THEN
  230. IF( T( ILST, ILST-1 ).NE.ZERO )
  231. $ ILST = ILST - 1
  232. END IF
  233. NBL = 1
  234. IF( ILST.LT.N ) THEN
  235. IF( T( ILST+1, ILST ).NE.ZERO )
  236. $ NBL = 2
  237. END IF
  238. *
  239. IF( IFST.EQ.ILST )
  240. $ RETURN
  241. *
  242. IF( IFST.LT.ILST ) THEN
  243. *
  244. * Update ILST
  245. *
  246. IF( NBF.EQ.2 .AND. NBL.EQ.1 )
  247. $ ILST = ILST - 1
  248. IF( NBF.EQ.1 .AND. NBL.EQ.2 )
  249. $ ILST = ILST + 1
  250. *
  251. HERE = IFST
  252. *
  253. 10 CONTINUE
  254. *
  255. * Swap block with next one below
  256. *
  257. IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
  258. *
  259. * Current block either 1 by 1 or 2 by 2
  260. *
  261. NBNEXT = 1
  262. IF( HERE+NBF+1.LE.N ) THEN
  263. IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
  264. $ NBNEXT = 2
  265. END IF
  266. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
  267. $ WORK, INFO )
  268. IF( INFO.NE.0 ) THEN
  269. ILST = HERE
  270. RETURN
  271. END IF
  272. HERE = HERE + NBNEXT
  273. *
  274. * Test if 2 by 2 block breaks into two 1 by 1 blocks
  275. *
  276. IF( NBF.EQ.2 ) THEN
  277. IF( T( HERE+1, HERE ).EQ.ZERO )
  278. $ NBF = 3
  279. END IF
  280. *
  281. ELSE
  282. *
  283. * Current block consists of two 1 by 1 blocks each of which
  284. * must be swapped individually
  285. *
  286. NBNEXT = 1
  287. IF( HERE+3.LE.N ) THEN
  288. IF( T( HERE+3, HERE+2 ).NE.ZERO )
  289. $ NBNEXT = 2
  290. END IF
  291. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
  292. $ WORK, INFO )
  293. IF( INFO.NE.0 ) THEN
  294. ILST = HERE
  295. RETURN
  296. END IF
  297. IF( NBNEXT.EQ.1 ) THEN
  298. *
  299. * Swap two 1 by 1 blocks, no problems possible
  300. *
  301. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
  302. $ WORK, INFO )
  303. HERE = HERE + 1
  304. ELSE
  305. *
  306. * Recompute NBNEXT in case 2 by 2 split
  307. *
  308. IF( T( HERE+2, HERE+1 ).EQ.ZERO )
  309. $ NBNEXT = 1
  310. IF( NBNEXT.EQ.2 ) THEN
  311. *
  312. * 2 by 2 Block did not split
  313. *
  314. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
  315. $ NBNEXT, WORK, INFO )
  316. IF( INFO.NE.0 ) THEN
  317. ILST = HERE
  318. RETURN
  319. END IF
  320. HERE = HERE + 2
  321. ELSE
  322. *
  323. * 2 by 2 Block did split
  324. *
  325. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
  326. $ WORK, INFO )
  327. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
  328. $ WORK, INFO )
  329. HERE = HERE + 2
  330. END IF
  331. END IF
  332. END IF
  333. IF( HERE.LT.ILST )
  334. $ GO TO 10
  335. *
  336. ELSE
  337. *
  338. HERE = IFST
  339. 20 CONTINUE
  340. *
  341. * Swap block with next one above
  342. *
  343. IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
  344. *
  345. * Current block either 1 by 1 or 2 by 2
  346. *
  347. NBNEXT = 1
  348. IF( HERE.GE.3 ) THEN
  349. IF( T( HERE-1, HERE-2 ).NE.ZERO )
  350. $ NBNEXT = 2
  351. END IF
  352. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
  353. $ NBF, WORK, INFO )
  354. IF( INFO.NE.0 ) THEN
  355. ILST = HERE
  356. RETURN
  357. END IF
  358. HERE = HERE - NBNEXT
  359. *
  360. * Test if 2 by 2 block breaks into two 1 by 1 blocks
  361. *
  362. IF( NBF.EQ.2 ) THEN
  363. IF( T( HERE+1, HERE ).EQ.ZERO )
  364. $ NBF = 3
  365. END IF
  366. *
  367. ELSE
  368. *
  369. * Current block consists of two 1 by 1 blocks each of which
  370. * must be swapped individually
  371. *
  372. NBNEXT = 1
  373. IF( HERE.GE.3 ) THEN
  374. IF( T( HERE-1, HERE-2 ).NE.ZERO )
  375. $ NBNEXT = 2
  376. END IF
  377. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
  378. $ 1, WORK, INFO )
  379. IF( INFO.NE.0 ) THEN
  380. ILST = HERE
  381. RETURN
  382. END IF
  383. IF( NBNEXT.EQ.1 ) THEN
  384. *
  385. * Swap two 1 by 1 blocks, no problems possible
  386. *
  387. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
  388. $ WORK, INFO )
  389. HERE = HERE - 1
  390. ELSE
  391. *
  392. * Recompute NBNEXT in case 2 by 2 split
  393. *
  394. IF( T( HERE, HERE-1 ).EQ.ZERO )
  395. $ NBNEXT = 1
  396. IF( NBNEXT.EQ.2 ) THEN
  397. *
  398. * 2 by 2 Block did not split
  399. *
  400. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
  401. $ WORK, INFO )
  402. IF( INFO.NE.0 ) THEN
  403. ILST = HERE
  404. RETURN
  405. END IF
  406. HERE = HERE - 2
  407. ELSE
  408. *
  409. * 2 by 2 Block did split
  410. *
  411. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
  412. $ WORK, INFO )
  413. CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
  414. $ WORK, INFO )
  415. HERE = HERE - 2
  416. END IF
  417. END IF
  418. END IF
  419. IF( HERE.GT.ILST )
  420. $ GO TO 20
  421. END IF
  422. ILST = HERE
  423. *
  424. RETURN
  425. *
  426. * End of DTREXC
  427. *
  428. END
  429.  
  430.  
  431.  

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