LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
slaqr2.f
Go to the documentation of this file.
1 *> \brief \b SLAQR2 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 SLAQR2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqr2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqr2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqr2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAQR2( 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 *> SLAQR2 is identical to SLAQR3 except that it avoids
43 *> recursion by calling SLAHQR instead of SLAQR4.
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 REAL 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 REAL 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 REAL array, dimension KBOT
173 *> \endverbatim
174 *>
175 *> \param[out] SI
176 *> \verbatim
177 *> SI is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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; SLAQR2
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 September 2012
267 *
268 *> \ingroup realOTHERauxiliary
269 *
270 *> \par Contributors:
271 * ==================
272 *>
273 *> Karen Braman and Ralph Byers, Department of Mathematics,
274 *> University of Kansas, USA
275 *>
276 * =====================================================================
277  SUBROUTINE slaqr2( 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 *
281 * -- LAPACK auxiliary routine (version 3.4.2) --
282 * -- LAPACK is a software package provided by Univ. of Tennessee, --
283 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
284 * September 2012
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  REAL H( ldh, * ), SI( * ), SR( * ), T( ldt, * ),
293  $ v( ldv, * ), work( * ), wv( ldwv, * ),
294  $ z( ldz, * )
295 * ..
296 *
297 * ================================================================
298 * .. Parameters ..
299  REAL ZERO, ONE
300  parameter ( zero = 0.0e0, one = 1.0e0 )
301 * ..
302 * .. Local Scalars ..
303  REAL 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  REAL SLAMCH
312  EXTERNAL slamch
313 * ..
314 * .. External Subroutines ..
315  EXTERNAL scopy, sgehrd, sgemm, slabad, slacpy, slahqr,
317 * ..
318 * .. Intrinsic Functions ..
319  INTRINSIC abs, int, max, min, REAL, 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 SGEHRD ====
331 *
332  CALL sgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
333  lwk1 = int( work( 1 ) )
334 *
335 * ==== Workspace query call to SORMHR ====
336 *
337  CALL sormhr( '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 ) = REAL( 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 = slamch( 'SAFE MINIMUM' )
367  safmax = one / safmin
368  CALL slabad( safmin, safmax )
369  ulp = slamch( 'PRECISION' )
370  smlnum = safmin*( REAL( 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 slacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
408  CALL scopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
409 *
410  CALL slaset( 'A', jw, jw, zero, one, v, ldv )
411  CALL slahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
412  $ si( kwtop ), 1, jw, v, ldv, infqr )
413 *
414 * ==== STREXC 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 * . (STREXC can not fail in this case.) ====
453 *
454  ifst = ns
455  CALL strexc( '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, STREXC does the right thing with
477 * . ILST in case of a rare exchange failure. ====
478 *
479  ifst = ns
480  CALL strexc( '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 strexc( '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 slanv2( 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 scopy( ns, v, ldv, work, 1 )
595  beta = work( 1 )
596  CALL slarfg( ns, beta, work( 2 ), 1, tau )
597  work( 1 ) = one
598 *
599  CALL slaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
600 *
601  CALL slarf( 'L', ns, jw, work, 1, tau, t, ldt,
602  $ work( jw+1 ) )
603  CALL slarf( 'R', ns, ns, work, 1, tau, t, ldt,
604  $ work( jw+1 ) )
605  CALL slarf( 'R', jw, ns, work, 1, tau, v, ldv,
606  $ work( jw+1 ) )
607 *
608  CALL sgehrd( 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 slacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
617  CALL scopy( 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 sormhr( '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 sgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
637  $ ldh, v, ldv, zero, wv, ldwv )
638  CALL slacpy( '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 sgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
647  $ h( kwtop, kcol ), ldh, zero, t, ldt )
648  CALL slacpy( '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 sgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
659  $ ldz, v, ldv, zero, wv, ldwv )
660  CALL slacpy( '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 ) = REAL( lwkopt )
681 *
682 * ==== End of SLAQR2 ====
683 *
684  END
subroutine slahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
Definition: slahqr.f:209
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108
subroutine slaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
SLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
Definition: slaqr2.f:280
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:169
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine strexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
STREXC
Definition: strexc.f:148
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine slanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form...
Definition: slanv2.f:129
subroutine sormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMHR
Definition: sormhr.f:181
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53