LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
claqr2.f
Go to the documentation of this file.
1 *> \brief \b CLAQR2 performs the unitary 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 CLAQR2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claqr2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claqr2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claqr2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
22 * IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT,
23 * 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 * COMPLEX H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
32 * $ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> CLAQR2 is identical to CLAQR3 except that it avoids
42 *> recursion by calling CLAHQR instead of CLAQR4.
43 *>
44 *> Aggressive early deflation:
45 *>
46 *> This subroutine accepts as input an upper Hessenberg matrix
47 *> H and performs an unitary similarity transformation
48 *> designed to detect and deflate fully converged eigenvalues from
49 *> a trailing principal submatrix. On output H has been over-
50 *> written by a new Hessenberg matrix that is a perturbation of
51 *> an unitary similarity transformation of H. It is to be
52 *> hoped that the final version of H has many zero subdiagonal
53 *> entries.
54 *> \endverbatim
55 *
56 * Arguments:
57 * ==========
58 *
59 *> \param[in] WANTT
60 *> \verbatim
61 *> WANTT is LOGICAL
62 *> If .TRUE., then the Hessenberg matrix H is fully updated
63 *> so that the triangular Schur factor may be
64 *> computed (in cooperation with the calling subroutine).
65 *> If .FALSE., then only enough of H is updated to preserve
66 *> the eigenvalues.
67 *> \endverbatim
68 *>
69 *> \param[in] WANTZ
70 *> \verbatim
71 *> WANTZ is LOGICAL
72 *> If .TRUE., then the unitary matrix Z is updated so
73 *> so that the unitary Schur factor may be computed
74 *> (in cooperation with the calling subroutine).
75 *> If .FALSE., then Z is not referenced.
76 *> \endverbatim
77 *>
78 *> \param[in] N
79 *> \verbatim
80 *> N is INTEGER
81 *> The order of the matrix H and (if WANTZ is .TRUE.) the
82 *> order of the unitary matrix Z.
83 *> \endverbatim
84 *>
85 *> \param[in] KTOP
86 *> \verbatim
87 *> KTOP is INTEGER
88 *> It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
89 *> KBOT and KTOP together determine an isolated block
90 *> along the diagonal of the Hessenberg matrix.
91 *> \endverbatim
92 *>
93 *> \param[in] KBOT
94 *> \verbatim
95 *> KBOT is INTEGER
96 *> It is assumed without a check that either
97 *> KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together
98 *> determine an isolated block along the diagonal of the
99 *> Hessenberg matrix.
100 *> \endverbatim
101 *>
102 *> \param[in] NW
103 *> \verbatim
104 *> NW is INTEGER
105 *> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
106 *> \endverbatim
107 *>
108 *> \param[in,out] H
109 *> \verbatim
110 *> H is COMPLEX array, dimension (LDH,N)
111 *> On input the initial N-by-N section of H stores the
112 *> Hessenberg matrix undergoing aggressive early deflation.
113 *> On output H has been transformed by a unitary
114 *> similarity transformation, perturbed, and the returned
115 *> to Hessenberg form that (it is to be hoped) has some
116 *> zero subdiagonal entries.
117 *> \endverbatim
118 *>
119 *> \param[in] LDH
120 *> \verbatim
121 *> LDH is integer
122 *> Leading dimension of H just as declared in the calling
123 *> subroutine. N .LE. LDH
124 *> \endverbatim
125 *>
126 *> \param[in] ILOZ
127 *> \verbatim
128 *> ILOZ is INTEGER
129 *> \endverbatim
130 *>
131 *> \param[in] IHIZ
132 *> \verbatim
133 *> IHIZ is INTEGER
134 *> Specify the rows of Z to which transformations must be
135 *> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
136 *> \endverbatim
137 *>
138 *> \param[in,out] Z
139 *> \verbatim
140 *> Z is COMPLEX array, dimension (LDZ,N)
141 *> IF WANTZ is .TRUE., then on output, the unitary
142 *> similarity transformation mentioned above has been
143 *> accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
144 *> If WANTZ is .FALSE., then Z is unreferenced.
145 *> \endverbatim
146 *>
147 *> \param[in] LDZ
148 *> \verbatim
149 *> LDZ is integer
150 *> The leading dimension of Z just as declared in the
151 *> calling subroutine. 1 .LE. LDZ.
152 *> \endverbatim
153 *>
154 *> \param[out] NS
155 *> \verbatim
156 *> NS is integer
157 *> The number of unconverged (ie approximate) eigenvalues
158 *> returned in SR and SI that may be used as shifts by the
159 *> calling subroutine.
160 *> \endverbatim
161 *>
162 *> \param[out] ND
163 *> \verbatim
164 *> ND is integer
165 *> The number of converged eigenvalues uncovered by this
166 *> subroutine.
167 *> \endverbatim
168 *>
169 *> \param[out] SH
170 *> \verbatim
171 *> SH is COMPLEX array, dimension KBOT
172 *> On output, approximate eigenvalues that may
173 *> be used for shifts are stored in SH(KBOT-ND-NS+1)
174 *> through SR(KBOT-ND). Converged eigenvalues are
175 *> stored in SH(KBOT-ND+1) through SH(KBOT).
176 *> \endverbatim
177 *>
178 *> \param[out] V
179 *> \verbatim
180 *> V is COMPLEX array, dimension (LDV,NW)
181 *> An NW-by-NW work array.
182 *> \endverbatim
183 *>
184 *> \param[in] LDV
185 *> \verbatim
186 *> LDV is integer scalar
187 *> The leading dimension of V just as declared in the
188 *> calling subroutine. NW .LE. LDV
189 *> \endverbatim
190 *>
191 *> \param[in] NH
192 *> \verbatim
193 *> NH is integer scalar
194 *> The number of columns of T. NH.GE.NW.
195 *> \endverbatim
196 *>
197 *> \param[out] T
198 *> \verbatim
199 *> T is COMPLEX array, dimension (LDT,NW)
200 *> \endverbatim
201 *>
202 *> \param[in] LDT
203 *> \verbatim
204 *> LDT is integer
205 *> The leading dimension of T just as declared in the
206 *> calling subroutine. NW .LE. LDT
207 *> \endverbatim
208 *>
209 *> \param[in] NV
210 *> \verbatim
211 *> NV is integer
212 *> The number of rows of work array WV available for
213 *> workspace. NV.GE.NW.
214 *> \endverbatim
215 *>
216 *> \param[out] WV
217 *> \verbatim
218 *> WV is COMPLEX array, dimension (LDWV,NW)
219 *> \endverbatim
220 *>
221 *> \param[in] LDWV
222 *> \verbatim
223 *> LDWV is integer
224 *> The leading dimension of W just as declared in the
225 *> calling subroutine. NW .LE. LDV
226 *> \endverbatim
227 *>
228 *> \param[out] WORK
229 *> \verbatim
230 *> WORK is COMPLEX array, dimension LWORK.
231 *> On exit, WORK(1) is set to an estimate of the optimal value
232 *> of LWORK for the given values of N, NW, KTOP and KBOT.
233 *> \endverbatim
234 *>
235 *> \param[in] LWORK
236 *> \verbatim
237 *> LWORK is integer
238 *> The dimension of the work array WORK. LWORK = 2*NW
239 *> suffices, but greater efficiency may result from larger
240 *> values of LWORK.
241 *>
242 *> If LWORK = -1, then a workspace query is assumed; CLAQR2
243 *> only estimates the optimal workspace size for the given
244 *> values of N, NW, KTOP and KBOT. The estimate is returned
245 *> in WORK(1). No error message related to LWORK is issued
246 *> by XERBLA. Neither H nor Z are accessed.
247 *> \endverbatim
248 *
249 * Authors:
250 * ========
251 *
252 *> \author Univ. of Tennessee
253 *> \author Univ. of California Berkeley
254 *> \author Univ. of Colorado Denver
255 *> \author NAG Ltd.
256 *
257 *> \date September 2012
258 *
259 *> \ingroup complexOTHERauxiliary
260 *
261 *> \par Contributors:
262 * ==================
263 *>
264 *> Karen Braman and Ralph Byers, Department of Mathematics,
265 *> University of Kansas, USA
266 *>
267 * =====================================================================
268  SUBROUTINE claqr2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
269  $ ihiz, z, ldz, ns, nd, sh, v, ldv, nh, t, ldt,
270  $ nv, wv, ldwv, work, lwork )
271 *
272 * -- LAPACK auxiliary routine (version 3.4.2) --
273 * -- LAPACK is a software package provided by Univ. of Tennessee, --
274 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
275 * September 2012
276 *
277 * .. Scalar Arguments ..
278  INTEGER ihiz, iloz, kbot, ktop, ldh, ldt, ldv, ldwv,
279  $ ldz, lwork, n, nd, nh, ns, nv, nw
280  LOGICAL wantt, wantz
281 * ..
282 * .. Array Arguments ..
283  COMPLEX h( ldh, * ), sh( * ), t( ldt, * ), v( ldv, * ),
284  $ work( * ), wv( ldwv, * ), z( ldz, * )
285 * ..
286 *
287 * ================================================================
288 *
289 * .. Parameters ..
290  COMPLEX zero, one
291  parameter( zero = ( 0.0e0, 0.0e0 ),
292  $ one = ( 1.0e0, 0.0e0 ) )
293  REAL rzero, rone
294  parameter( rzero = 0.0e0, rone = 1.0e0 )
295 * ..
296 * .. Local Scalars ..
297  COMPLEX beta, cdum, s, tau
298  REAL foo, safmax, safmin, smlnum, ulp
299  INTEGER i, ifst, ilst, info, infqr, j, jw, kcol, kln,
300  $ knt, krow, kwtop, ltop, lwk1, lwk2, lwkopt
301 * ..
302 * .. External Functions ..
303  REAL slamch
304  EXTERNAL slamch
305 * ..
306 * .. External Subroutines ..
307  EXTERNAL ccopy, cgehrd, cgemm, clacpy, clahqr, clarf,
309 * ..
310 * .. Intrinsic Functions ..
311  INTRINSIC abs, aimag, cmplx, conjg, int, max, min, real
312 * ..
313 * .. Statement Functions ..
314  REAL cabs1
315 * ..
316 * .. Statement Function definitions ..
317  cabs1( cdum ) = abs( REAL( CDUM ) ) + abs( aimag( cdum ) )
318 * ..
319 * .. Executable Statements ..
320 *
321 * ==== Estimate optimal workspace. ====
322 *
323  jw = min( nw, kbot-ktop+1 )
324  IF( jw.LE.2 ) THEN
325  lwkopt = 1
326  ELSE
327 *
328 * ==== Workspace query call to CGEHRD ====
329 *
330  CALL cgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
331  lwk1 = int( work( 1 ) )
332 *
333 * ==== Workspace query call to CUNMHR ====
334 *
335  CALL cunmhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
336  $ work, -1, info )
337  lwk2 = int( work( 1 ) )
338 *
339 * ==== Optimal workspace ====
340 *
341  lwkopt = jw + max( lwk1, lwk2 )
342  END IF
343 *
344 * ==== Quick return in case of workspace query. ====
345 *
346  IF( lwork.EQ.-1 ) THEN
347  work( 1 ) = cmplx( lwkopt, 0 )
348  return
349  END IF
350 *
351 * ==== Nothing to do ...
352 * ... for an empty active block ... ====
353  ns = 0
354  nd = 0
355  work( 1 ) = one
356  IF( ktop.GT.kbot )
357  $ return
358 * ... nor for an empty deflation window. ====
359  IF( nw.LT.1 )
360  $ return
361 *
362 * ==== Machine constants ====
363 *
364  safmin = slamch( 'SAFE MINIMUM' )
365  safmax = rone / safmin
366  CALL slabad( safmin, safmax )
367  ulp = slamch( 'PRECISION' )
368  smlnum = safmin*( REAL( N ) / ulp )
369 *
370 * ==== Setup deflation window ====
371 *
372  jw = min( nw, kbot-ktop+1 )
373  kwtop = kbot - jw + 1
374  IF( kwtop.EQ.ktop ) THEN
375  s = zero
376  ELSE
377  s = h( kwtop, kwtop-1 )
378  END IF
379 *
380  IF( kbot.EQ.kwtop ) THEN
381 *
382 * ==== 1-by-1 deflation window: not much to do ====
383 *
384  sh( kwtop ) = h( kwtop, kwtop )
385  ns = 1
386  nd = 0
387  IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
388  $ kwtop ) ) ) ) THEN
389  ns = 0
390  nd = 1
391  IF( kwtop.GT.ktop )
392  $ h( kwtop, kwtop-1 ) = zero
393  END IF
394  work( 1 ) = one
395  return
396  END IF
397 *
398 * ==== Convert to spike-triangular form. (In case of a
399 * . rare QR failure, this routine continues to do
400 * . aggressive early deflation using that part of
401 * . the deflation window that converged using INFQR
402 * . here and there to keep track.) ====
403 *
404  CALL clacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
405  CALL ccopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
406 *
407  CALL claset( 'A', jw, jw, zero, one, v, ldv )
408  CALL clahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
409  $ jw, v, ldv, infqr )
410 *
411 * ==== Deflation detection loop ====
412 *
413  ns = jw
414  ilst = infqr + 1
415  DO 10 knt = infqr + 1, jw
416 *
417 * ==== Small spike tip deflation test ====
418 *
419  foo = cabs1( t( ns, ns ) )
420  IF( foo.EQ.rzero )
421  $ foo = cabs1( s )
422  IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
423  $ THEN
424 *
425 * ==== One more converged eigenvalue ====
426 *
427  ns = ns - 1
428  ELSE
429 *
430 * ==== One undeflatable eigenvalue. Move it up out of the
431 * . way. (CTREXC can not fail in this case.) ====
432 *
433  ifst = ns
434  CALL ctrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
435  ilst = ilst + 1
436  END IF
437  10 continue
438 *
439 * ==== Return to Hessenberg form ====
440 *
441  IF( ns.EQ.0 )
442  $ s = zero
443 *
444  IF( ns.LT.jw ) THEN
445 *
446 * ==== sorting the diagonal of T improves accuracy for
447 * . graded matrices. ====
448 *
449  DO 30 i = infqr + 1, ns
450  ifst = i
451  DO 20 j = i + 1, ns
452  IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
453  $ ifst = j
454  20 continue
455  ilst = i
456  IF( ifst.NE.ilst )
457  $ CALL ctrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
458  30 continue
459  END IF
460 *
461 * ==== Restore shift/eigenvalue array from T ====
462 *
463  DO 40 i = infqr + 1, jw
464  sh( kwtop+i-1 ) = t( i, i )
465  40 continue
466 *
467 *
468  IF( ns.LT.jw .OR. s.EQ.zero ) THEN
469  IF( ns.GT.1 .AND. s.NE.zero ) THEN
470 *
471 * ==== Reflect spike back into lower triangle ====
472 *
473  CALL ccopy( ns, v, ldv, work, 1 )
474  DO 50 i = 1, ns
475  work( i ) = conjg( work( i ) )
476  50 continue
477  beta = work( 1 )
478  CALL clarfg( ns, beta, work( 2 ), 1, tau )
479  work( 1 ) = one
480 *
481  CALL claset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
482 *
483  CALL clarf( 'L', ns, jw, work, 1, conjg( tau ), t, ldt,
484  $ work( jw+1 ) )
485  CALL clarf( 'R', ns, ns, work, 1, tau, t, ldt,
486  $ work( jw+1 ) )
487  CALL clarf( 'R', jw, ns, work, 1, tau, v, ldv,
488  $ work( jw+1 ) )
489 *
490  CALL cgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
491  $ lwork-jw, info )
492  END IF
493 *
494 * ==== Copy updated reduced window into place ====
495 *
496  IF( kwtop.GT.1 )
497  $ h( kwtop, kwtop-1 ) = s*conjg( v( 1, 1 ) )
498  CALL clacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
499  CALL ccopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
500  $ ldh+1 )
501 *
502 * ==== Accumulate orthogonal matrix in order update
503 * . H and Z, if requested. ====
504 *
505  IF( ns.GT.1 .AND. s.NE.zero )
506  $ CALL cunmhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
507  $ work( jw+1 ), lwork-jw, info )
508 *
509 * ==== Update vertical slab in H ====
510 *
511  IF( wantt ) THEN
512  ltop = 1
513  ELSE
514  ltop = ktop
515  END IF
516  DO 60 krow = ltop, kwtop - 1, nv
517  kln = min( nv, kwtop-krow )
518  CALL cgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
519  $ ldh, v, ldv, zero, wv, ldwv )
520  CALL clacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
521  60 continue
522 *
523 * ==== Update horizontal slab in H ====
524 *
525  IF( wantt ) THEN
526  DO 70 kcol = kbot + 1, n, nh
527  kln = min( nh, n-kcol+1 )
528  CALL cgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
529  $ h( kwtop, kcol ), ldh, zero, t, ldt )
530  CALL clacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
531  $ ldh )
532  70 continue
533  END IF
534 *
535 * ==== Update vertical slab in Z ====
536 *
537  IF( wantz ) THEN
538  DO 80 krow = iloz, ihiz, nv
539  kln = min( nv, ihiz-krow+1 )
540  CALL cgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
541  $ ldz, v, ldv, zero, wv, ldwv )
542  CALL clacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
543  $ ldz )
544  80 continue
545  END IF
546  END IF
547 *
548 * ==== Return the number of deflations ... ====
549 *
550  nd = jw - ns
551 *
552 * ==== ... and the number of shifts. (Subtracting
553 * . INFQR from the spike length takes care
554 * . of the case of a rare QR failure while
555 * . calculating eigenvalues of the deflation
556 * . window.) ====
557 *
558  ns = ns - infqr
559 *
560 * ==== Return optimal workspace. ====
561 *
562  work( 1 ) = cmplx( lwkopt, 0 )
563 *
564 * ==== End of CLAQR2 ====
565 *
566  END