Télécharger dlaexc.eso

Retour à la liste

Numérotation des lignes :

  1. C DLAEXC SOURCE BP208322 15/10/13 21:15:27 8670
  2. *> \brief \b DLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form, by an orthogonal similarity transformation.
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DLAEXC + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaexc.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaexc.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaexc.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
  23. * INFO )
  24. *
  25. * .. Scalar Arguments ..
  26. * LOGICAL WANTQ
  27. * INTEGER INFO, J1, LDQ, LDT, N, N1, N2
  28. * ..
  29. * .. Array Arguments ..
  30. * REAL*8 Q( LDQ, * ), T( LDT, * ), WORK( * )
  31. * ..
  32. *
  33. *
  34. *> \par Purpose:
  35. * =============
  36. *>
  37. *> \verbatim
  38. *>
  39. *> DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
  40. *> an upper quasi-triangular matrix T by an orthogonal similarity
  41. *> transformation.
  42. *>
  43. *> T must be in Schur canonical form, that is, block upper triangular
  44. *> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
  45. *> has its diagonal elemnts equal and its off-diagonal elements of
  46. *> opposite sign.
  47. *> \endverbatim
  48. *
  49. * Arguments:
  50. * ==========
  51. *
  52. *> \param[in] WANTQ
  53. *> \verbatim
  54. *> WANTQ is LOGICAL
  55. *> = .TRUE. : accumulate the transformation in the matrix Q;
  56. *> = .FALSE.: do not accumulate the transformation.
  57. *> \endverbatim
  58. *>
  59. *> \param[in] N
  60. *> \verbatim
  61. *> N is INTEGER
  62. *> The order of the matrix T. N >= 0.
  63. *> \endverbatim
  64. *>
  65. *> \param[in,out] T
  66. *> \verbatim
  67. *> T is DOUBLE PRECISION array, dimension (LDT,N)
  68. *> On entry, the upper quasi-triangular matrix T, in Schur
  69. *> canonical form.
  70. *> On exit, the updated matrix T, again in Schur canonical form.
  71. *> \endverbatim
  72. *>
  73. *> \param[in] LDT
  74. *> \verbatim
  75. *> LDT is INTEGER
  76. *> The leading dimension of the array T. LDT >= max(1,N).
  77. *> \endverbatim
  78. *>
  79. *> \param[in,out] Q
  80. *> \verbatim
  81. *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
  82. *> On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
  83. *> On exit, if WANTQ is .TRUE., the updated matrix Q.
  84. *> If WANTQ is .FALSE., Q is not referenced.
  85. *> \endverbatim
  86. *>
  87. *> \param[in] LDQ
  88. *> \verbatim
  89. *> LDQ is INTEGER
  90. *> The leading dimension of the array Q.
  91. *> LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
  92. *> \endverbatim
  93. *>
  94. *> \param[in] J1
  95. *> \verbatim
  96. *> J1 is INTEGER
  97. *> The index of the first row of the first block T11.
  98. *> \endverbatim
  99. *>
  100. *> \param[in] N1
  101. *> \verbatim
  102. *> N1 is INTEGER
  103. *> The order of the first block T11. N1 = 0, 1 or 2.
  104. *> \endverbatim
  105. *>
  106. *> \param[in] N2
  107. *> \verbatim
  108. *> N2 is INTEGER
  109. *> The order of the second block T22. N2 = 0, 1 or 2.
  110. *> \endverbatim
  111. *>
  112. *> \param[out] WORK
  113. *> \verbatim
  114. *> WORK is DOUBLE PRECISION array, dimension (N)
  115. *> \endverbatim
  116. *>
  117. *> \param[out] INFO
  118. *> \verbatim
  119. *> INFO is INTEGER
  120. *> = 0: successful exit
  121. *> = 1: the transformed matrix T would be too far from Schur
  122. *> form; the blocks are not swapped and T and Q are
  123. *> unchanged.
  124. *> \endverbatim
  125. *
  126. * Authors:
  127. * ========
  128. *
  129. *> \author Univ. of Tennessee
  130. *> \author Univ. of California Berkeley
  131. *> \author Univ. of Colorado Denver
  132. *> \author NAG Ltd.
  133. *
  134. *> \date September 2012
  135. *
  136. *> \ingroup doubleOTHERauxiliary
  137. *
  138. * =====================================================================
  139. SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
  140. $ INFO )
  141. *
  142. * -- LAPACK auxiliary routine (version 3.4.2) --
  143. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  144. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  145. * September 2012
  146. *
  147. * .. Scalar Arguments ..
  148. LOGICAL WANTQ
  149. INTEGER INFO, J1, LDQ, LDT, N, N1, N2
  150. * ..
  151. * .. Array Arguments ..
  152. REAL*8 Q( LDQ, * ), T( LDT, * ), WORK( * )
  153. * ..
  154. *
  155. * =====================================================================
  156. *
  157. * .. Parameters ..
  158. REAL*8 ZERO, ONE
  159. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  160. REAL*8 TEN
  161. PARAMETER ( TEN = 1.0D+1 )
  162. INTEGER LDD, LDX
  163. PARAMETER ( LDD = 4, LDX = 2 )
  164. * ..
  165. * .. Local Scalars ..
  166. INTEGER IERR, J2, J3, J4, K, ND
  167. REAL*8 CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
  168. $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
  169. $ WR1, WR2, XNORM
  170. * ..
  171. * .. Local Arrays ..
  172. REAL*8 D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
  173. $ X( LDX, 2 )
  174. * ..
  175. * .. External Functions ..
  176. REAL*8 DLAMCH, DLANGE
  177. EXTERNAL DLAMCH, DLANGE
  178. * ..
  179. * .. External Subroutines ..
  180. $ DROT
  181. * ..
  182. ** .. Intrinsic Functions ..
  183. * INTRINSIC ABS, MAX
  184. ** ..
  185. ** .. Executable Statements ..
  186. *
  187. INFO = 0
  188. *
  189. * Quick return if possible
  190. *
  191. IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
  192. $ RETURN
  193. IF( J1+N1.GT.N )
  194. $ RETURN
  195. *
  196. J2 = J1 + 1
  197. J3 = J1 + 2
  198. J4 = J1 + 3
  199. *
  200. IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
  201. *
  202. * Swap two 1-by-1 blocks.
  203. *
  204. T11 = T( J1, J1 )
  205. T22 = T( J2, J2 )
  206. *
  207. * Determine the transformation to perform the interchange.
  208. *
  209. CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
  210. *
  211. * Apply transformation to the matrix T.
  212. *
  213. IF( J3.LE.N )
  214. $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
  215. $ SN )
  216. CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
  217. *
  218. T( J1, J1 ) = T22
  219. T( J2, J2 ) = T11
  220. *
  221. IF( WANTQ ) THEN
  222. *
  223. * Accumulate transformation in the matrix Q.
  224. *
  225. CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
  226. END IF
  227. *
  228. ELSE
  229. *
  230. * Swapping involves at least one 2-by-2 block.
  231. *
  232. * Copy the diagonal block of order N1+N2 to the local array D
  233. * and compute its norm.
  234. *
  235. ND = N1 + N2
  236. CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
  237. DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
  238. *
  239. * Compute machine-dependent threshold for test for accepting
  240. * swap.
  241. *
  242. EPS = DLAMCH( 'P' )
  243. SMLNUM = DLAMCH( 'S' ) / EPS
  244. THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
  245. *
  246. * Solve T11*X - X*T22 = scale*T12 for X.
  247. *
  248. CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
  249. $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
  250. $ LDX, XNORM, IERR )
  251. *
  252. * Swap the adjacent diagonal blocks.
  253. *
  254. K = N1 + N1 + N2 - 3
  255. GO TO ( 10, 20, 30 )K
  256. *
  257. 10 CONTINUE
  258. *
  259. * N1 = 1, N2 = 2: generate elementary reflector H so that:
  260. *
  261. * ( scale, X11, X12 ) H = ( 0, 0, * )
  262. *
  263. U( 1 ) = SCALE
  264. U( 2 ) = X( 1, 1 )
  265. U( 3 ) = X( 1, 2 )
  266. CALL DLARFG( 3, U( 3 ), U, 1, TAU )
  267. U( 3 ) = ONE
  268. T11 = T( J1, J1 )
  269. *
  270. * Perform swap provisionally on diagonal block in D.
  271. *
  272. CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
  273. CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
  274. *
  275. * Test whether to reject swap.
  276. *
  277. IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
  278. $ 3 )-T11 ) ).GT.THRESH )GO TO 50
  279. *
  280. * Accept swap: apply transformation to the entire matrix T.
  281. *
  282. CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
  283. CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
  284. *
  285. T( J3, J1 ) = ZERO
  286. T( J3, J2 ) = ZERO
  287. T( J3, J3 ) = T11
  288. *
  289. IF( WANTQ ) THEN
  290. *
  291. * Accumulate transformation in the matrix Q.
  292. *
  293. CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
  294. END IF
  295. GO TO 40
  296. *
  297. 20 CONTINUE
  298. *
  299. * N1 = 2, N2 = 1: generate elementary reflector H so that:
  300. *
  301. * H ( -X11 ) = ( * )
  302. * ( -X21 ) = ( 0 )
  303. * ( scale ) = ( 0 )
  304. *
  305. U( 1 ) = -X( 1, 1 )
  306. U( 2 ) = -X( 2, 1 )
  307. U( 3 ) = SCALE
  308. CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
  309. U( 1 ) = ONE
  310. T33 = T( J3, J3 )
  311. *
  312. * Perform swap provisionally on diagonal block in D.
  313. *
  314. CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
  315. CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
  316. *
  317. * Test whether to reject swap.
  318. *
  319. IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
  320. $ 1 )-T33 ) ).GT.THRESH )GO TO 50
  321. *
  322. * Accept swap: apply transformation to the entire matrix T.
  323. *
  324. CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
  325. CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
  326. *
  327. T( J1, J1 ) = T33
  328. T( J2, J1 ) = ZERO
  329. T( J3, J1 ) = ZERO
  330. *
  331. IF( WANTQ ) THEN
  332. *
  333. * Accumulate transformation in the matrix Q.
  334. *
  335. CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
  336. END IF
  337. GO TO 40
  338. *
  339. 30 CONTINUE
  340. *
  341. * N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
  342. * that:
  343. *
  344. * H(2) H(1) ( -X11 -X12 ) = ( * * )
  345. * ( -X21 -X22 ) ( 0 * )
  346. * ( scale 0 ) ( 0 0 )
  347. * ( 0 scale ) ( 0 0 )
  348. *
  349. U1( 1 ) = -X( 1, 1 )
  350. U1( 2 ) = -X( 2, 1 )
  351. U1( 3 ) = SCALE
  352. CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
  353. U1( 1 ) = ONE
  354. *
  355. TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
  356. U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
  357. U2( 2 ) = -TEMP*U1( 3 )
  358. U2( 3 ) = SCALE
  359. CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
  360. U2( 1 ) = ONE
  361. *
  362. * Perform swap provisionally on diagonal block in D.
  363. *
  364. CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
  365. CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
  366. CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
  367. CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
  368. *
  369. * Test whether to reject swap.
  370. *
  371. IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
  372. $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
  373. *
  374. * Accept swap: apply transformation to the entire matrix T.
  375. *
  376. CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
  377. CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
  378. CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
  379. CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
  380. *
  381. T( J3, J1 ) = ZERO
  382. T( J3, J2 ) = ZERO
  383. T( J4, J1 ) = ZERO
  384. T( J4, J2 ) = ZERO
  385. *
  386. IF( WANTQ ) THEN
  387. *
  388. * Accumulate transformation in the matrix Q.
  389. *
  390. CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
  391. CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
  392. END IF
  393. *
  394. 40 CONTINUE
  395. *
  396. IF( N2.EQ.2 ) THEN
  397. *
  398. * Standardize new 2-by-2 block T11
  399. *
  400. CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
  401. $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
  402. CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
  403. $ CS, SN )
  404. CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
  405. IF( WANTQ )
  406. $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
  407. END IF
  408. *
  409. IF( N1.EQ.2 ) THEN
  410. *
  411. * Standardize new 2-by-2 block T22
  412. *
  413. J3 = J1 + N2
  414. J4 = J3 + 1
  415. CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
  416. $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
  417. IF( J3+2.LE.N )
  418. $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
  419. $ LDT, CS, SN )
  420. CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
  421. IF( WANTQ )
  422. $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
  423. END IF
  424. *
  425. END IF
  426. RETURN
  427. *
  428. * Exit with INFO = 1 if swap was rejected.
  429. *
  430. 50 CONTINUE
  431. INFO = 1
  432. RETURN
  433. *
  434. * End of DLAEXC
  435. *
  436. END
  437.  
  438.  

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