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