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