Télécharger dlahqr.eso

Retour à la liste

Numérotation des lignes :

  1. C DLAHQR SOURCE BP208322 15/10/13 21:15:27 8670
  2. SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
  3. $ ILOZ, IHIZ, Z, LDZ, INFO )
  4. *
  5. * -- LAPACK auxiliary routine (version 2.0) --
  6. * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  7. * Courant Institute, Argonne National Lab, and Rice University
  8. * October 31, 1992
  9. *
  10. * .. Scalar Arguments ..
  11. LOGICAL WANTT, WANTZ
  12. INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
  13. * ..
  14. * .. Array Arguments ..
  15. REAL*8 H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
  16. * ..
  17. *
  18. * Purpose
  19. * =======
  20. *
  21. * DLAHQR is an auxiliary routine called by DHSEQR to update the
  22. * eigenvalues and Schur decomposition already computed by DHSEQR, by
  23. * dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
  24. *
  25. * Arguments
  26. * =========
  27. *
  28. * WANTT (input) LOGICAL
  29. * = .TRUE. : the full Schur form T is required;
  30. * = .FALSE.: only eigenvalues are required.
  31. *
  32. * WANTZ (input) LOGICAL
  33. * = .TRUE. : the matrix of Schur vectors Z is required;
  34. * = .FALSE.: Schur vectors are not required.
  35. *
  36. * N (input) INTEGER
  37. * The order of the matrix H. N >= 0.
  38. *
  39. * ILO (input) INTEGER
  40. * IHI (input) INTEGER
  41. * It is assumed that H is already upper quasi-triangular in
  42. * rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
  43. * ILO = 1). DLAHQR works primarily with the Hessenberg
  44. * submatrix in rows and columns ILO to IHI, but applies
  45. * transformations to all of H if WANTT is .TRUE..
  46. * 1 <= ILO <= max(1,IHI); IHI <= N.
  47. *
  48. * H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
  49. * On entry, the upper Hessenberg matrix H.
  50. * On exit, if WANTT is .TRUE., H is upper quasi-triangular in
  51. * rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in
  52. * standard form. If WANTT is .FALSE., the contents of H are
  53. * unspecified on exit.
  54. *
  55. * LDH (input) INTEGER
  56. * The leading dimension of the array H. LDH >= max(1,N).
  57. *
  58. * WR (output) DOUBLE PRECISION array, dimension (N)
  59. * WI (output) DOUBLE PRECISION array, dimension (N)
  60. * The real and imaginary parts, respectively, of the computed
  61. * eigenvalues ILO to IHI are stored in the corresponding
  62. * elements of WR and WI. If two eigenvalues are computed as a
  63. * complex conjugate pair, they are stored in consecutive
  64. * elements of WR and WI, say the i-th and (i+1)th, with
  65. * WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
  66. * eigenvalues are stored in the same order as on the diagonal
  67. * of the Schur form returned in H, with WR(i) = H(i,i), and, if
  68. * H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
  69. * WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
  70. *
  71. * ILOZ (input) INTEGER
  72. * IHIZ (input) INTEGER
  73. * Specify the rows of Z to which transformations must be
  74. * applied if WANTZ is .TRUE..
  75. * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
  76. *
  77. * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
  78. * If WANTZ is .TRUE., on entry Z must contain the current
  79. * matrix Z of transformations accumulated by DHSEQR, and on
  80. * exit Z has been updated; transformations are applied only to
  81. * the submatrix Z(ILOZ:IHIZ,ILO:IHI).
  82. * If WANTZ is .FALSE., Z is not referenced.
  83. *
  84. * LDZ (input) INTEGER
  85. * The leading dimension of the array Z. LDZ >= max(1,N).
  86. *
  87. * INFO (output) INTEGER
  88. * = 0: successful exit
  89. * > 0: DLAHQR failed to compute all the eigenvalues ILO to IHI
  90. * in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
  91. * elements i+1:ihi of WR and WI contain those eigenvalues
  92. * which have been successfully computed.
  93. *
  94. * =====================================================================
  95. *
  96. * .. Parameters ..
  97. REAL*8 ZERO, ONE
  98. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  99. REAL*8 DAT1, DAT2
  100. PARAMETER ( DAT1 = 0.75D+0, DAT2 = -0.4375D+0 )
  101. * ..
  102. * .. Local Scalars ..
  103. INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
  104. REAL*8 CS, H00, H10, H11, H12, H21, H22, H33, H33S,
  105. $ H43H34, H44, H44S, OVFL, S, SMLNUM, SN, SUM,
  106. $ T1, T2, T3, TST1, ULP, UNFL, V1, V2, V3
  107. * ..
  108. * .. Local Arrays ..
  109. REAL*8 V( 3 ), WORK( 1 )
  110. * ..
  111. * .. External Functions ..
  112. REAL*8 DLAMCH, DLANHS
  113. EXTERNAL DLAMCH, DLANHS
  114. * ..
  115. * .. External Subroutines ..
  116. EXTERNAL DCOPY, DLABAD, DLANV2, DLARFG, DROT
  117. * ..
  118. ** .. Intrinsic Functions ..
  119. * INTRINSIC ABS, MAX, MIN
  120. ** ..
  121. ** .. Executable Statements ..
  122. *
  123. INFO = 0
  124. *
  125. * Quick return if possible
  126. *
  127. IF( N.EQ.0 )
  128. $ RETURN
  129. IF( ILO.EQ.IHI ) THEN
  130. WR( ILO ) = H( ILO, ILO )
  131. WI( ILO ) = ZERO
  132. RETURN
  133. END IF
  134. *
  135. NH = IHI - ILO + 1
  136. NZ = IHIZ - ILOZ + 1
  137. *
  138. * Set machine-dependent constants for the stopping criterion.
  139. * If norm(H) <= sqrt(OVFL), overflow should not occur.
  140. *
  141. UNFL = DLAMCH( 'Safe minimum' )
  142. OVFL = ONE / UNFL
  143. CALL DLABAD( UNFL, OVFL )
  144. ULP = DLAMCH( 'Precision' )
  145. SMLNUM = UNFL*( NH / ULP )
  146. *
  147. * I1 and I2 are the indices of the first row and last column of H
  148. * to which transformations must be applied. If eigenvalues only are
  149. * being computed, I1 and I2 are set inside the main loop.
  150. *
  151. IF( WANTT ) THEN
  152. I1 = 1
  153. I2 = N
  154. END IF
  155. *
  156. * ITN is the total number of QR iterations allowed.
  157. *
  158. ITN = 30*NH
  159. *
  160. * The main loop begins here. I is the loop index and decreases from
  161. * IHI to ILO in steps of 1 or 2. Each iteration of the loop works
  162. * with the active submatrix in rows and columns L to I.
  163. * Eigenvalues I+1 to IHI have already converged. Either L = ILO or
  164. * H(L,L-1) is negligible so that the matrix splits.
  165. *
  166. I = IHI
  167. 10 CONTINUE
  168. L = ILO
  169. IF( I.LT.ILO )
  170. $ GO TO 150
  171. *
  172. * Perform QR iterations on rows and columns ILO to I until a
  173. * submatrix of order 1 or 2 splits off at the bottom because a
  174. * subdiagonal element has become negligible.
  175. *
  176. DO 130 ITS = 0, ITN
  177. *
  178. * Look for a single small subdiagonal element.
  179. *
  180. DO 20 K = I, L + 1, -1
  181. TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
  182. IF( TST1.EQ.ZERO )
  183. $ TST1 = DLANHS( '1', I-L+1, H( L, L ), LDH, WORK )
  184. IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) )
  185. $ GO TO 30
  186. 20 CONTINUE
  187. 30 CONTINUE
  188. L = K
  189. IF( L.GT.ILO ) THEN
  190. *
  191. * H(L,L-1) is negligible
  192. *
  193. H( L, L-1 ) = ZERO
  194. END IF
  195. *
  196. * Exit from loop if a submatrix of order 1 or 2 has split off.
  197. *
  198. IF( L.GE.I-1 )
  199. $ GO TO 140
  200. *
  201. * Now the active submatrix is in rows and columns L to I. If
  202. * eigenvalues only are being computed, only the active submatrix
  203. * need be transformed.
  204. *
  205. IF( .NOT.WANTT ) THEN
  206. I1 = L
  207. I2 = I
  208. END IF
  209. *
  210. IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN
  211. *
  212. * Exceptional shift.
  213. *
  214. S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
  215. H44 = DAT1*S
  216. H33 = H44
  217. H43H34 = DAT2*S*S
  218. ELSE
  219. *
  220. * Prepare to use Wilkinson's double shift
  221. *
  222. H44 = H( I, I )
  223. H33 = H( I-1, I-1 )
  224. H43H34 = H( I, I-1 )*H( I-1, I )
  225. END IF
  226. *
  227. * Look for two consecutive small subdiagonal elements.
  228. *
  229. DO 40 M = I - 2, L, -1
  230. *
  231. * Determine the effect of starting the double-shift QR
  232. * iteration at row M, and see if this would make H(M,M-1)
  233. * negligible.
  234. *
  235. H11 = H( M, M )
  236. H22 = H( M+1, M+1 )
  237. H21 = H( M+1, M )
  238. H12 = H( M, M+1 )
  239. H44S = H44 - H11
  240. H33S = H33 - H11
  241. V1 = ( H33S*H44S-H43H34 ) / H21 + H12
  242. V2 = H22 - H11 - H33S - H44S
  243. V3 = H( M+2, M+1 )
  244. S = ABS( V1 ) + ABS( V2 ) + ABS( V3 )
  245. V1 = V1 / S
  246. V2 = V2 / S
  247. V3 = V3 / S
  248. V( 1 ) = V1
  249. V( 2 ) = V2
  250. V( 3 ) = V3
  251. IF( M.EQ.L )
  252. $ GO TO 50
  253. H00 = H( M-1, M-1 )
  254. H10 = H( M, M-1 )
  255. TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) )
  256. IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) ).LE.ULP*TST1 )
  257. $ GO TO 50
  258. 40 CONTINUE
  259. 50 CONTINUE
  260. *
  261. * Double-shift QR step
  262. *
  263. DO 120 K = M, I - 1
  264. *
  265. * The first iteration of this loop determines a reflection G
  266. * from the vector V and applies it from left and right to H,
  267. * thus creating a nonzero bulge below the subdiagonal.
  268. *
  269. * Each subsequent iteration determines a reflection G to
  270. * restore the Hessenberg form in the (K-1)th column, and thus
  271. * chases the bulge one step toward the bottom of the active
  272. * submatrix. NR is the order of G.
  273. *
  274. NR = MIN( 3, I-K+1 )
  275. IF( K.GT.M )
  276. $ CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 )
  277. CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 )
  278. IF( K.GT.M ) THEN
  279. H( K, K-1 ) = V( 1 )
  280. H( K+1, K-1 ) = ZERO
  281. IF( K.LT.I-1 )
  282. $ H( K+2, K-1 ) = ZERO
  283. ELSE IF( M.GT.L ) THEN
  284. H( K, K-1 ) = -H( K, K-1 )
  285. END IF
  286. V2 = V( 2 )
  287. T2 = T1*V2
  288. IF( NR.EQ.3 ) THEN
  289. V3 = V( 3 )
  290. T3 = T1*V3
  291. *
  292. * Apply G from the left to transform the rows of the matrix
  293. * in columns K to I2.
  294. *
  295. DO 60 J = K, I2
  296. SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )
  297. H( K, J ) = H( K, J ) - SUM*T1
  298. H( K+1, J ) = H( K+1, J ) - SUM*T2
  299. H( K+2, J ) = H( K+2, J ) - SUM*T3
  300. 60 CONTINUE
  301. *
  302. * Apply G from the right to transform the columns of the
  303. * matrix in rows I1 to min(K+3,I).
  304. *
  305. DO 70 J = I1, MIN( K+3, I )
  306. SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )
  307. H( J, K ) = H( J, K ) - SUM*T1
  308. H( J, K+1 ) = H( J, K+1 ) - SUM*T2
  309. H( J, K+2 ) = H( J, K+2 ) - SUM*T3
  310. 70 CONTINUE
  311. *
  312. IF( WANTZ ) THEN
  313. *
  314. * Accumulate transformations in the matrix Z
  315. *
  316. DO 80 J = ILOZ, IHIZ
  317. SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )
  318. Z( J, K ) = Z( J, K ) - SUM*T1
  319. Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
  320. Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3
  321. 80 CONTINUE
  322. END IF
  323. ELSE IF( NR.EQ.2 ) THEN
  324. *
  325. * Apply G from the left to transform the rows of the matrix
  326. * in columns K to I2.
  327. *
  328. DO 90 J = K, I2
  329. SUM = H( K, J ) + V2*H( K+1, J )
  330. H( K, J ) = H( K, J ) - SUM*T1
  331. H( K+1, J ) = H( K+1, J ) - SUM*T2
  332. 90 CONTINUE
  333. *
  334. * Apply G from the right to transform the columns of the
  335. * matrix in rows I1 to min(K+3,I).
  336. *
  337. DO 100 J = I1, I
  338. SUM = H( J, K ) + V2*H( J, K+1 )
  339. H( J, K ) = H( J, K ) - SUM*T1
  340. H( J, K+1 ) = H( J, K+1 ) - SUM*T2
  341. 100 CONTINUE
  342. *
  343. IF( WANTZ ) THEN
  344. *
  345. * Accumulate transformations in the matrix Z
  346. *
  347. DO 110 J = ILOZ, IHIZ
  348. SUM = Z( J, K ) + V2*Z( J, K+1 )
  349. Z( J, K ) = Z( J, K ) - SUM*T1
  350. Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
  351. 110 CONTINUE
  352. END IF
  353. END IF
  354. 120 CONTINUE
  355. *
  356. 130 CONTINUE
  357. *
  358. * Failure to converge in remaining number of iterations
  359. *
  360. INFO = I
  361. RETURN
  362. *
  363. 140 CONTINUE
  364. *
  365. IF( L.EQ.I ) THEN
  366. *
  367. * H(I,I-1) is negligible: one eigenvalue has converged.
  368. *
  369. WR( I ) = H( I, I )
  370. WI( I ) = ZERO
  371. ELSE IF( L.EQ.I-1 ) THEN
  372. *
  373. * H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
  374. *
  375. * Transform the 2-by-2 submatrix to standard Schur form,
  376. * and compute and store the eigenvalues.
  377. *
  378. CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),
  379. $ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),
  380. $ CS, SN )
  381. *
  382. IF( WANTT ) THEN
  383. *
  384. * Apply the transformation to the rest of H.
  385. *
  386. IF( I2.GT.I )
  387. $ CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,
  388. $ CS, SN )
  389. CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
  390. END IF
  391. IF( WANTZ ) THEN
  392. *
  393. * Apply the transformation to Z.
  394. *
  395. CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
  396. END IF
  397. END IF
  398. *
  399. * Decrement number of remaining iterations, and return to start of
  400. * the main loop with new value of I.
  401. *
  402. ITN = ITN - ITS
  403. I = L - 1
  404. GO TO 10
  405. *
  406. 150 CONTINUE
  407. RETURN
  408. *
  409. * End of DLAHQR
  410. *
  411. END
  412.  
  413.  

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