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