Télécharger dlaqr2.eso

Retour à la liste

Numérotation des lignes :

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

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