Télécharger dlaqr3.eso

Retour à la liste

Numérotation des lignes :

  1. C DLAQR3 SOURCE BP208322 20/09/18 21:16:01 10718
  2. *> \brief \b DLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DLAQR3 + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr3.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr3.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr3.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
  23. * IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
  24. * LDT, NV, WV, LDWV, WORK, LWORK )
  25. *
  26. * .. Scalar Arguments ..
  27. * INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
  28. * $ LDZ, LWORK, N, ND, NH, NS, NV, NW
  29. * LOGICAL WANTT, WANTZ
  30. * ..
  31. * .. Array Arguments ..
  32. * REAL*8 H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
  33. * $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
  34. * $ Z( LDZ, * )
  35. * ..
  36. *
  37. *
  38. *> \par Purpose:
  39. * =============
  40. *>
  41. *> \verbatim
  42. *>
  43. *> Aggressive early deflation:
  44. *>
  45. *> DLAQR3 accepts as input an upper Hessenberg matrix
  46. *> H and performs an orthogonal similarity transformation
  47. *> designed to detect and deflate fully converged eigenvalues from
  48. *> a trailing principal submatrix. On output H has been over-
  49. *> written by a new Hessenberg matrix that is a perturbation of
  50. *> an orthogonal similarity transformation of H. It is to be
  51. *> hoped that the final version of H has many zero subdiagonal
  52. *> entries.
  53. *> \endverbatim
  54. *
  55. * Arguments:
  56. * ==========
  57. *
  58. *> \param[in] WANTT
  59. *> \verbatim
  60. *> WANTT is LOGICAL
  61. *> If .TRUE., then the Hessenberg matrix H is fully updated
  62. *> so that the quasi-triangular Schur factor may be
  63. *> computed (in cooperation with the calling subroutine).
  64. *> If .FALSE., then only enough of H is updated to preserve
  65. *> the eigenvalues.
  66. *> \endverbatim
  67. *>
  68. *> \param[in] WANTZ
  69. *> \verbatim
  70. *> WANTZ is LOGICAL
  71. *> If .TRUE., then the orthogonal matrix Z is updated so
  72. *> so that the orthogonal Schur factor may be computed
  73. *> (in cooperation with the calling subroutine).
  74. *> If .FALSE., then Z is not referenced.
  75. *> \endverbatim
  76. *>
  77. *> \param[in] N
  78. *> \verbatim
  79. *> N is INTEGER
  80. *> The order of the matrix H and (if WANTZ is .TRUE.) the
  81. *> order of the orthogonal matrix Z.
  82. *> \endverbatim
  83. *>
  84. *> \param[in] KTOP
  85. *> \verbatim
  86. *> KTOP is INTEGER
  87. *> It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
  88. *> KBOT and KTOP together determine an isolated block
  89. *> along the diagonal of the Hessenberg matrix.
  90. *> \endverbatim
  91. *>
  92. *> \param[in] KBOT
  93. *> \verbatim
  94. *> KBOT is INTEGER
  95. *> It is assumed without a check that either
  96. *> KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together
  97. *> determine an isolated block along the diagonal of the
  98. *> Hessenberg matrix.
  99. *> \endverbatim
  100. *>
  101. *> \param[in] NW
  102. *> \verbatim
  103. *> NW is INTEGER
  104. *> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
  105. *> \endverbatim
  106. *>
  107. *> \param[in,out] H
  108. *> \verbatim
  109. *> H is REAL*8 array, dimension (LDH,N)
  110. *> On input the initial N-by-N section of H stores the
  111. *> Hessenberg matrix undergoing aggressive early deflation.
  112. *> On output H has been transformed by an orthogonal
  113. *> similarity transformation, perturbed, and the returned
  114. *> to Hessenberg form that (it is to be hoped) has some
  115. *> zero subdiagonal entries.
  116. *> \endverbatim
  117. *>
  118. *> \param[in] LDH
  119. *> \verbatim
  120. *> LDH is INTEGER
  121. *> Leading dimension of H just as declared in the calling
  122. *> subroutine. N .LE. LDH
  123. *> \endverbatim
  124. *>
  125. *> \param[in] ILOZ
  126. *> \verbatim
  127. *> ILOZ is INTEGER
  128. *> \endverbatim
  129. *>
  130. *> \param[in] IHIZ
  131. *> \verbatim
  132. *> IHIZ is INTEGER
  133. *> Specify the rows of Z to which transformations must be
  134. *> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
  135. *> \endverbatim
  136. *>
  137. *> \param[in,out] Z
  138. *> \verbatim
  139. *> Z is REAL*8 array, dimension (LDZ,N)
  140. *> IF WANTZ is .TRUE., then on output, the orthogonal
  141. *> similarity transformation mentioned above has been
  142. *> accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
  143. *> If WANTZ is .FALSE., then Z is unreferenced.
  144. *> \endverbatim
  145. *>
  146. *> \param[in] LDZ
  147. *> \verbatim
  148. *> LDZ is INTEGER
  149. *> The leading dimension of Z just as declared in the
  150. *> calling subroutine. 1 .LE. LDZ.
  151. *> \endverbatim
  152. *>
  153. *> \param[out] NS
  154. *> \verbatim
  155. *> NS is INTEGER
  156. *> The number of unconverged (ie approximate) eigenvalues
  157. *> returned in SR and SI that may be used as shifts by the
  158. *> calling subroutine.
  159. *> \endverbatim
  160. *>
  161. *> \param[out] ND
  162. *> \verbatim
  163. *> ND is INTEGER
  164. *> The number of converged eigenvalues uncovered by this
  165. *> subroutine.
  166. *> \endverbatim
  167. *>
  168. *> \param[out] SR
  169. *> \verbatim
  170. *> SR is REAL*8 array, dimension (KBOT)
  171. *> \endverbatim
  172. *>
  173. *> \param[out] SI
  174. *> \verbatim
  175. *> SI is REAL*8 array, dimension (KBOT)
  176. *> On output, the real and imaginary parts of approximate
  177. *> eigenvalues that may be used for shifts are stored in
  178. *> SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
  179. *> SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
  180. *> The real and imaginary parts of converged eigenvalues
  181. *> are stored in SR(KBOT-ND+1) through SR(KBOT) and
  182. *> SI(KBOT-ND+1) through SI(KBOT), respectively.
  183. *> \endverbatim
  184. *>
  185. *> \param[out] V
  186. *> \verbatim
  187. *> V is REAL*8 array, dimension (LDV,NW)
  188. *> An NW-by-NW work array.
  189. *> \endverbatim
  190. *>
  191. *> \param[in] LDV
  192. *> \verbatim
  193. *> LDV is INTEGER
  194. *> The leading dimension of V just as declared in the
  195. *> calling subroutine. NW .LE. LDV
  196. *> \endverbatim
  197. *>
  198. *> \param[in] NH
  199. *> \verbatim
  200. *> NH is INTEGER
  201. *> The number of columns of T. NH.GE.NW.
  202. *> \endverbatim
  203. *>
  204. *> \param[out] T
  205. *> \verbatim
  206. *> T is REAL*8 array, dimension (LDT,NW)
  207. *> \endverbatim
  208. *>
  209. *> \param[in] LDT
  210. *> \verbatim
  211. *> LDT is INTEGER
  212. *> The leading dimension of T just as declared in the
  213. *> calling subroutine. NW .LE. LDT
  214. *> \endverbatim
  215. *>
  216. *> \param[in] NV
  217. *> \verbatim
  218. *> NV is INTEGER
  219. *> The number of rows of work array WV available for
  220. *> workspace. NV.GE.NW.
  221. *> \endverbatim
  222. *>
  223. *> \param[out] WV
  224. *> \verbatim
  225. *> WV is REAL*8 array, dimension (LDWV,NW)
  226. *> \endverbatim
  227. *>
  228. *> \param[in] LDWV
  229. *> \verbatim
  230. *> LDWV is INTEGER
  231. *> The leading dimension of W just as declared in the
  232. *> calling subroutine. NW .LE. LDV
  233. *> \endverbatim
  234. *>
  235. *> \param[out] WORK
  236. *> \verbatim
  237. *> WORK is REAL*8 array, dimension (LWORK)
  238. *> On exit, WORK(1) is set to an estimate of the optimal value
  239. *> of LWORK for the given values of N, NW, KTOP and KBOT.
  240. *> \endverbatim
  241. *>
  242. *> \param[in] LWORK
  243. *> \verbatim
  244. *> LWORK is INTEGER
  245. *> The dimension of the work array WORK. LWORK = 2*NW
  246. *> suffices, but greater efficiency may result from larger
  247. *> values of LWORK.
  248. *>
  249. *> If LWORK = -1, then a workspace query is assumed; DLAQR3
  250. *> only estimates the optimal workspace size for the given
  251. *> values of N, NW, KTOP and KBOT. The estimate is returned
  252. *> in WORK(1). No error message related to LWORK is issued
  253. *> by XERBLA. Neither H nor Z are accessed.
  254. *> \endverbatim
  255. *
  256. * Authors:
  257. * ========
  258. *
  259. *> \author Univ. of Tennessee
  260. *> \author Univ. of California Berkeley
  261. *> \author Univ. of Colorado Denver
  262. *> \author NAG Ltd.
  263. *
  264. *> \date June 2016
  265. *
  266. *> \ingroup doubleOTHERauxiliary
  267. *
  268. *> \par Contributors:
  269. * ==================
  270. *>
  271. *> Karen Braman and Ralph Byers, Department of Mathematics,
  272. *> University of Kansas, USA
  273. *>
  274. * =====================================================================
  275. SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
  276. $ IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
  277. $ LDT, NV, WV, LDWV, WORK, LWORK )
  278. *
  279. * -- LAPACK auxiliary routine (version 3.7.1) --
  280. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  281. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  282. * June 2016
  283.  
  284. IMPLICIT INTEGER(I-N)
  285. IMPLICIT REAL*8(A-H,O-Z)
  286. *
  287. * .. Scalar Arguments ..
  288. INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
  289. $ LDZ, LWORK, N, ND, NH, NS, NV, NW
  290. LOGICAL WANTT, WANTZ
  291. * ..
  292. * .. Array Arguments ..
  293. REAL*8 H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
  294. $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
  295. $ Z( LDZ, * )
  296. * ..
  297. *
  298. * ================================================================
  299. * .. Parameters ..
  300. REAL*8 ZERO, ONE
  301. PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
  302. * ..
  303. * .. Local Scalars ..
  304. REAL*8 AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
  305. $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
  306. INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
  307. $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3,
  308. $ LWKOPT, NMIN
  309. LOGICAL BULGE, SORTED
  310. * ..
  311. * .. External Functions ..
  312. REAL*8 DLAMCH
  313. INTEGER ILAENV
  314. * EXTERNAL DLAMCH, ILAENV
  315. * ..
  316. * .. External Subroutines ..
  317. * EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR,
  318. * $ DLANV2, DLAQR4, DLARF, DLARFG, DLASET, DORMHR,
  319. * $ DTREXC
  320. * ..
  321. * .. Intrinsic Functions ..
  322. * INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT
  323. * ..
  324. * .. Executable Statements ..
  325. *
  326. * ==== Estimate optimal workspace. ====
  327. *
  328. JW = MIN( NW, KBOT-KTOP+1 )
  329. IF( JW.LE.2 ) THEN
  330. LWKOPT = 1
  331. ELSE
  332. *
  333. * ==== Workspace query call to DGEHRD ====
  334. *
  335. CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
  336. LWK1 = INT( WORK( 1 ) )
  337. *
  338. * ==== Workspace query call to DORMHR ====
  339. *
  340. CALL DORMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV,
  341. $ WORK, -1, INFO )
  342. LWK2 = INT( WORK( 1 ) )
  343. *
  344. * ==== Workspace query call to DLAQR4 ====
  345. *
  346. CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR, SI, 1, JW,
  347. $ V, LDV, WORK, -1, INFQR )
  348. LWK3 = INT( WORK( 1 ) )
  349. *
  350. * ==== Optimal workspace ====
  351. *
  352. LWKOPT = MAX( JW+MAX( LWK1, LWK2 ), LWK3 )
  353. END IF
  354. *
  355. * ==== Quick return in case of workspace query. ====
  356. *
  357. IF( LWORK.EQ.-1 ) THEN
  358. WORK( 1 ) = DBLE( LWKOPT )
  359. RETURN
  360. END IF
  361. *
  362. * ==== Nothing to do ...
  363. * ... for an empty active block ... ====
  364. NS = 0
  365. ND = 0
  366. WORK( 1 ) = ONE
  367. IF( KTOP.GT.KBOT )
  368. $ RETURN
  369. * ... nor for an empty deflation window. ====
  370. IF( NW.LT.1 )
  371. $ RETURN
  372. *
  373. * ==== Machine constants ====
  374. *
  375. SAFMIN = DLAMCH( 'SAFE MINIMUM' )
  376. SAFMAX = ONE / SAFMIN
  377. CALL DLABAD( SAFMIN, SAFMAX )
  378. ULP = DLAMCH( 'PRECISION' )
  379. SMLNUM = SAFMIN*( DBLE( N ) / ULP )
  380. *
  381. * ==== Setup deflation window ====
  382. *
  383. JW = MIN( NW, KBOT-KTOP+1 )
  384. KWTOP = KBOT - JW + 1
  385. IF( KWTOP.EQ.KTOP ) THEN
  386. S = ZERO
  387. ELSE
  388. S = H( KWTOP, KWTOP-1 )
  389. END IF
  390. *
  391. IF( KBOT.EQ.KWTOP ) THEN
  392. *
  393. * ==== 1-by-1 deflation window: not much to do ====
  394. *
  395. SR( KWTOP ) = H( KWTOP, KWTOP )
  396. SI( KWTOP ) = ZERO
  397. NS = 1
  398. ND = 0
  399. IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) )
  400. $ THEN
  401. NS = 0
  402. ND = 1
  403. IF( KWTOP.GT.KTOP )
  404. $ H( KWTOP, KWTOP-1 ) = ZERO
  405. END IF
  406. WORK( 1 ) = ONE
  407. RETURN
  408. END IF
  409. *
  410. * ==== Convert to spike-triangular form. (In case of a
  411. * . rare QR failure, this routine continues to do
  412. * . aggressive early deflation using that part of
  413. * . the deflation window that converged using INFQR
  414. * . here and there to keep track.) ====
  415. *
  416. CALL DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
  417. CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
  418. *
  419. CALL DLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
  420. NMIN = ILAENV( 12, 'DLAQR3', 'SV', JW, 1, JW, LWORK )
  421. IF( JW.GT.NMIN ) THEN
  422. CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
  423. $ SI( KWTOP ), 1, JW, V, LDV, WORK, LWORK, INFQR )
  424. ELSE
  425. CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
  426. $ SI( KWTOP ), 1, JW, V, LDV, INFQR )
  427. END IF
  428. *
  429. * ==== DTREXC needs a clean margin near the diagonal ====
  430. *
  431. DO 10 J = 1, JW - 3
  432. T( J+2, J ) = ZERO
  433. T( J+3, J ) = ZERO
  434. 10 CONTINUE
  435. IF( JW.GT.2 )
  436. $ T( JW, JW-2 ) = ZERO
  437. *
  438. * ==== Deflation detection loop ====
  439. *
  440. NS = JW
  441. ILST = INFQR + 1
  442. 20 CONTINUE
  443. IF( ILST.LE.NS ) THEN
  444. IF( NS.EQ.1 ) THEN
  445. BULGE = .FALSE.
  446. ELSE
  447. BULGE = T( NS, NS-1 ).NE.ZERO
  448. END IF
  449. *
  450. * ==== Small spike tip test for deflation ====
  451. *
  452. IF( .NOT. BULGE ) THEN
  453. *
  454. * ==== Real eigenvalue ====
  455. *
  456. FOO = ABS( T( NS, NS ) )
  457. IF( FOO.EQ.ZERO )
  458. $ FOO = ABS( S )
  459. IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN
  460. *
  461. * ==== Deflatable ====
  462. *
  463. NS = NS - 1
  464. ELSE
  465. *
  466. * ==== Undeflatable. Move it up out of the way.
  467. * . (DTREXC can not fail in this case.) ====
  468. *
  469. IFST = NS
  470. CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
  471. $ INFO )
  472. ILST = ILST + 1
  473. END IF
  474. ELSE
  475. *
  476. * ==== Complex conjugate pair ====
  477. *
  478. FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )*
  479. $ SQRT( ABS( T( NS-1, NS ) ) )
  480. IF( FOO.EQ.ZERO )
  481. $ FOO = ABS( S )
  482. IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE.
  483. $ MAX( SMLNUM, ULP*FOO ) ) THEN
  484. *
  485. * ==== Deflatable ====
  486. *
  487. NS = NS - 2
  488. ELSE
  489. *
  490. * ==== Undeflatable. Move them up out of the way.
  491. * . Fortunately, DTREXC does the right thing with
  492. * . ILST in case of a rare exchange failure. ====
  493. *
  494. IFST = NS
  495. CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
  496. $ INFO )
  497. ILST = ILST + 2
  498. END IF
  499. END IF
  500. *
  501. * ==== End deflation detection loop ====
  502. *
  503. GO TO 20
  504. END IF
  505. *
  506. * ==== Return to Hessenberg form ====
  507. *
  508. IF( NS.EQ.0 )
  509. $ S = ZERO
  510. *
  511. IF( NS.LT.JW ) THEN
  512. *
  513. * ==== sorting diagonal blocks of T improves accuracy for
  514. * . graded matrices. Bubble sort deals well with
  515. * . exchange failures. ====
  516. *
  517. SORTED = .false.
  518. I = NS + 1
  519. 30 CONTINUE
  520. IF( SORTED )
  521. $ GO TO 50
  522. SORTED = .true.
  523. *
  524. KEND = I - 1
  525. I = INFQR + 1
  526. IF( I.EQ.NS ) THEN
  527. K = I + 1
  528. ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
  529. K = I + 1
  530. ELSE
  531. K = I + 2
  532. END IF
  533. 40 CONTINUE
  534. IF( K.LE.KEND ) THEN
  535. IF( K.EQ.I+1 ) THEN
  536. EVI = ABS( T( I, I ) )
  537. ELSE
  538. EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )*
  539. $ SQRT( ABS( T( I, I+1 ) ) )
  540. END IF
  541. *
  542. IF( K.EQ.KEND ) THEN
  543. EVK = ABS( T( K, K ) )
  544. ELSE IF( T( K+1, K ).EQ.ZERO ) THEN
  545. EVK = ABS( T( K, K ) )
  546. ELSE
  547. EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )*
  548. $ SQRT( ABS( T( K, K+1 ) ) )
  549. END IF
  550. *
  551. IF( EVI.GE.EVK ) THEN
  552. I = K
  553. ELSE
  554. SORTED = .false.
  555. IFST = I
  556. ILST = K
  557. CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
  558. $ INFO )
  559. IF( INFO.EQ.0 ) THEN
  560. I = ILST
  561. ELSE
  562. I = K
  563. END IF
  564. END IF
  565. IF( I.EQ.KEND ) THEN
  566. K = I + 1
  567. ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
  568. K = I + 1
  569. ELSE
  570. K = I + 2
  571. END IF
  572. GO TO 40
  573. END IF
  574. GO TO 30
  575. 50 CONTINUE
  576. END IF
  577. *
  578. * ==== Restore shift/eigenvalue array from T ====
  579. *
  580. I = JW
  581. 60 CONTINUE
  582. IF( I.GE.INFQR+1 ) THEN
  583. IF( I.EQ.INFQR+1 ) THEN
  584. SR( KWTOP+I-1 ) = T( I, I )
  585. SI( KWTOP+I-1 ) = ZERO
  586. I = I - 1
  587. ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN
  588. SR( KWTOP+I-1 ) = T( I, I )
  589. SI( KWTOP+I-1 ) = ZERO
  590. I = I - 1
  591. ELSE
  592. AA = T( I-1, I-1 )
  593. CC = T( I, I-1 )
  594. BB = T( I-1, I )
  595. DD = T( I, I )
  596. CALL DLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ),
  597. $ SI( KWTOP+I-2 ), SR( KWTOP+I-1 ),
  598. $ SI( KWTOP+I-1 ), CS, SN )
  599. I = I - 2
  600. END IF
  601. GO TO 60
  602. END IF
  603. *
  604. IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
  605. IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
  606. *
  607. * ==== Reflect spike back into lower triangle ====
  608. *
  609. CALL DCOPY( NS, V, LDV, WORK, 1 )
  610. BETA = WORK( 1 )
  611. CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU )
  612. WORK( 1 ) = ONE
  613. *
  614. CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
  615. *
  616. CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT,
  617. $ WORK( JW+1 ) )
  618. CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT,
  619. $ WORK( JW+1 ) )
  620. CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV,
  621. $ WORK( JW+1 ) )
  622. *
  623. CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
  624. $ LWORK-JW, INFO )
  625. END IF
  626. *
  627. * ==== Copy updated reduced window into place ====
  628. *
  629. IF( KWTOP.GT.1 )
  630. $ H( KWTOP, KWTOP-1 ) = S*V( 1, 1 )
  631. CALL DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )
  632. CALL DCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
  633. $ LDH+1 )
  634. *
  635. * ==== Accumulate orthogonal matrix in order update
  636. * . H and Z, if requested. ====
  637. *
  638. IF( NS.GT.1 .AND. S.NE.ZERO )
  639. $ CALL DORMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV,
  640. $ WORK( JW+1 ), LWORK-JW, INFO )
  641. *
  642. * ==== Update vertical slab in H ====
  643. *
  644. IF( WANTT ) THEN
  645. LTOP = 1
  646. ELSE
  647. LTOP = KTOP
  648. END IF
  649. DO 70 KROW = LTOP, KWTOP - 1, NV
  650. KLN = MIN( NV, KWTOP-KROW )
  651. CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
  652. $ LDH, V, LDV, ZERO, WV, LDWV )
  653. CALL DLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )
  654. 70 CONTINUE
  655. *
  656. * ==== Update horizontal slab in H ====
  657. *
  658. IF( WANTT ) THEN
  659. DO 80 KCOL = KBOT + 1, N, NH
  660. KLN = MIN( NH, N-KCOL+1 )
  661. CALL DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV,
  662. $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT )
  663. CALL DLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ),
  664. $ LDH )
  665. 80 CONTINUE
  666. END IF
  667. *
  668. * ==== Update vertical slab in Z ====
  669. *
  670. IF( WANTZ ) THEN
  671. DO 90 KROW = ILOZ, IHIZ, NV
  672. KLN = MIN( NV, IHIZ-KROW+1 )
  673. CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
  674. $ LDZ, V, LDV, ZERO, WV, LDWV )
  675. CALL DLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
  676. $ LDZ )
  677. 90 CONTINUE
  678. END IF
  679. END IF
  680. *
  681. * ==== Return the number of deflations ... ====
  682. *
  683. ND = JW - NS
  684. *
  685. * ==== ... and the number of shifts. (Subtracting
  686. * . INFQR from the spike length takes care
  687. * . of the case of a rare QR failure while
  688. * . calculating eigenvalues of the deflation
  689. * . window.) ====
  690. *
  691. NS = NS - INFQR
  692. *
  693. * ==== Return optimal workspace. ====
  694. *
  695. WORK( 1 ) = DBLE( LWKOPT )
  696. *
  697. * ==== End of DLAQR3 ====
  698. *
  699. END
  700.  
  701.  
  702.  

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