LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine zlaqr2 ( logical WANTT, logical WANTZ, integer N, integer KTOP, integer KBOT, integer NW, complex*16, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer IHIZ, complex*16, dimension( ldz, * ) Z, integer LDZ, integer NS, integer ND, complex*16, dimension( * ) SH, complex*16, dimension( ldv, * ) V, integer LDV, integer NH, complex*16, dimension( ldt, * ) T, integer LDT, integer NV, complex*16, dimension( ldwv, * ) WV, integer LDWV, complex*16, dimension( * ) WORK, integer LWORK )

ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).

Purpose:
```    ZLAQR2 is identical to ZLAQR3 except that it avoids
recursion by calling ZLAHQR instead of ZLAQR4.

Aggressive early deflation:

ZLAQR2 accepts as input an upper Hessenberg matrix
H and performs an unitary similarity transformation
designed to detect and deflate fully converged eigenvalues from
a trailing principal submatrix.  On output H has been over-
written by a new Hessenberg matrix that is a perturbation of
an unitary similarity transformation of H.  It is to be
hoped that the final version of H has many zero subdiagonal
entries.```
Parameters
 [in] WANTT ``` WANTT is LOGICAL If .TRUE., then the Hessenberg matrix H is fully updated so that the triangular Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then only enough of H is updated to preserve the eigenvalues.``` [in] WANTZ ``` WANTZ is LOGICAL If .TRUE., then the unitary matrix Z is updated so so that the unitary Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then Z is not referenced.``` [in] N ``` N is INTEGER The order of the matrix H and (if WANTZ is .TRUE.) the order of the unitary matrix Z.``` [in] KTOP ``` KTOP is INTEGER It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix.``` [in] KBOT ``` KBOT is INTEGER It is assumed without a check that either KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix.``` [in] NW ``` NW is INTEGER Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).``` [in,out] H ``` H is COMPLEX*16 array, dimension (LDH,N) On input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. On output H has been transformed by a unitary similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries.``` [in] LDH ``` LDH is integer Leading dimension of H just as declared in the calling subroutine. N .LE. LDH``` [in] ILOZ ` ILOZ is INTEGER` [in] IHIZ ``` IHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.``` [in,out] Z ``` Z is COMPLEX*16 array, dimension (LDZ,N) IF WANTZ is .TRUE., then on output, the unitary similarity transformation mentioned above has been accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right. If WANTZ is .FALSE., then Z is unreferenced.``` [in] LDZ ``` LDZ is integer The leading dimension of Z just as declared in the calling subroutine. 1 .LE. LDZ.``` [out] NS ``` NS is integer The number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine.``` [out] ND ``` ND is integer The number of converged eigenvalues uncovered by this subroutine.``` [out] SH ``` SH is COMPLEX*16 array, dimension KBOT On output, approximate eigenvalues that may be used for shifts are stored in SH(KBOT-ND-NS+1) through SR(KBOT-ND). Converged eigenvalues are stored in SH(KBOT-ND+1) through SH(KBOT).``` [out] V ``` V is COMPLEX*16 array, dimension (LDV,NW) An NW-by-NW work array.``` [in] LDV ``` LDV is integer scalar The leading dimension of V just as declared in the calling subroutine. NW .LE. LDV``` [in] NH ``` NH is integer scalar The number of columns of T. NH.GE.NW.``` [out] T ` T is COMPLEX*16 array, dimension (LDT,NW)` [in] LDT ``` LDT is integer The leading dimension of T just as declared in the calling subroutine. NW .LE. LDT``` [in] NV ``` NV is integer The number of rows of work array WV available for workspace. NV.GE.NW.``` [out] WV ` WV is COMPLEX*16 array, dimension (LDWV,NW)` [in] LDWV ``` LDWV is integer The leading dimension of W just as declared in the calling subroutine. NW .LE. LDV``` [out] WORK ``` WORK is COMPLEX*16 array, dimension LWORK. On exit, WORK(1) is set to an estimate of the optimal value of LWORK for the given values of N, NW, KTOP and KBOT.``` [in] LWORK ``` LWORK is integer The dimension of the work array WORK. LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK. If LWORK = -1, then a workspace query is assumed; ZLAQR2 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed.```
Date
September 2012
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

Definition at line 272 of file zlaqr2.f.

272 *
273 * -- LAPACK auxiliary routine (version 3.4.2) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * September 2012
277 *
278 * .. Scalar Arguments ..
279  INTEGER ihiz, iloz, kbot, ktop, ldh, ldt, ldv, ldwv,
280  \$ ldz, lwork, n, nd, nh, ns, nv, nw
281  LOGICAL wantt, wantz
282 * ..
283 * .. Array Arguments ..
284  COMPLEX*16 h( ldh, * ), sh( * ), t( ldt, * ), v( ldv, * ),
285  \$ work( * ), wv( ldwv, * ), z( ldz, * )
286 * ..
287 *
288 * ================================================================
289 *
290 * .. Parameters ..
291  COMPLEX*16 zero, one
292  parameter ( zero = ( 0.0d0, 0.0d0 ),
293  \$ one = ( 1.0d0, 0.0d0 ) )
294  DOUBLE PRECISION rzero, rone
295  parameter ( rzero = 0.0d0, rone = 1.0d0 )
296 * ..
297 * .. Local Scalars ..
298  COMPLEX*16 beta, cdum, s, tau
299  DOUBLE PRECISION foo, safmax, safmin, smlnum, ulp
300  INTEGER i, ifst, ilst, info, infqr, j, jw, kcol, kln,
301  \$ knt, krow, kwtop, ltop, lwk1, lwk2, lwkopt
302 * ..
303 * .. External Functions ..
304  DOUBLE PRECISION dlamch
305  EXTERNAL dlamch
306 * ..
307 * .. External Subroutines ..
308  EXTERNAL dlabad, zcopy, zgehrd, zgemm, zlacpy, zlahqr,
310 * ..
311 * .. Intrinsic Functions ..
312  INTRINSIC abs, dble, dcmplx, dconjg, dimag, int, max, min
313 * ..
314 * .. Statement Functions ..
315  DOUBLE PRECISION cabs1
316 * ..
317 * .. Statement Function definitions ..
318  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
319 * ..
320 * .. Executable Statements ..
321 *
322 * ==== Estimate optimal workspace. ====
323 *
324  jw = min( nw, kbot-ktop+1 )
325  IF( jw.LE.2 ) THEN
326  lwkopt = 1
327  ELSE
328 *
329 * ==== Workspace query call to ZGEHRD ====
330 *
331  CALL zgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
332  lwk1 = int( work( 1 ) )
333 *
334 * ==== Workspace query call to ZUNMHR ====
335 *
336  CALL zunmhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
337  \$ work, -1, info )
338  lwk2 = int( work( 1 ) )
339 *
340 * ==== Optimal workspace ====
341 *
342  lwkopt = jw + max( lwk1, lwk2 )
343  END IF
344 *
345 * ==== Quick return in case of workspace query. ====
346 *
347  IF( lwork.EQ.-1 ) THEN
348  work( 1 ) = dcmplx( lwkopt, 0 )
349  RETURN
350  END IF
351 *
352 * ==== Nothing to do ...
353 * ... for an empty active block ... ====
354  ns = 0
355  nd = 0
356  work( 1 ) = one
357  IF( ktop.GT.kbot )
358  \$ RETURN
359 * ... nor for an empty deflation window. ====
360  IF( nw.LT.1 )
361  \$ RETURN
362 *
363 * ==== Machine constants ====
364 *
365  safmin = dlamch( 'SAFE MINIMUM' )
366  safmax = rone / safmin
367  CALL dlabad( safmin, safmax )
368  ulp = dlamch( 'PRECISION' )
369  smlnum = safmin*( dble( n ) / ulp )
370 *
371 * ==== Setup deflation window ====
372 *
373  jw = min( nw, kbot-ktop+1 )
374  kwtop = kbot - jw + 1
375  IF( kwtop.EQ.ktop ) THEN
376  s = zero
377  ELSE
378  s = h( kwtop, kwtop-1 )
379  END IF
380 *
381  IF( kbot.EQ.kwtop ) THEN
382 *
383 * ==== 1-by-1 deflation window: not much to do ====
384 *
385  sh( kwtop ) = h( kwtop, kwtop )
386  ns = 1
387  nd = 0
388  IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
389  \$ kwtop ) ) ) ) THEN
390  ns = 0
391  nd = 1
392  IF( kwtop.GT.ktop )
393  \$ h( kwtop, kwtop-1 ) = zero
394  END IF
395  work( 1 ) = one
396  RETURN
397  END IF
398 *
399 * ==== Convert to spike-triangular form. (In case of a
400 * . rare QR failure, this routine continues to do
401 * . aggressive early deflation using that part of
402 * . the deflation window that converged using INFQR
403 * . here and there to keep track.) ====
404 *
405  CALL zlacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
406  CALL zcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
407 *
408  CALL zlaset( 'A', jw, jw, zero, one, v, ldv )
409  CALL zlahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
410  \$ jw, v, ldv, infqr )
411 *
412 * ==== Deflation detection loop ====
413 *
414  ns = jw
415  ilst = infqr + 1
416  DO 10 knt = infqr + 1, jw
417 *
418 * ==== Small spike tip deflation test ====
419 *
420  foo = cabs1( t( ns, ns ) )
421  IF( foo.EQ.rzero )
422  \$ foo = cabs1( s )
423  IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
424  \$ THEN
425 *
426 * ==== One more converged eigenvalue ====
427 *
428  ns = ns - 1
429  ELSE
430 *
431 * ==== One undeflatable eigenvalue. Move it up out of the
432 * . way. (ZTREXC can not fail in this case.) ====
433 *
434  ifst = ns
435  CALL ztrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
436  ilst = ilst + 1
437  END IF
438  10 CONTINUE
439 *
441 *
442  IF( ns.EQ.0 )
443  \$ s = zero
444 *
445  IF( ns.LT.jw ) THEN
446 *
447 * ==== sorting the diagonal of T improves accuracy for
448 * . graded matrices. ====
449 *
450  DO 30 i = infqr + 1, ns
451  ifst = i
452  DO 20 j = i + 1, ns
453  IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
454  \$ ifst = j
455  20 CONTINUE
456  ilst = i
457  IF( ifst.NE.ilst )
458  \$ CALL ztrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
459  30 CONTINUE
460  END IF
461 *
462 * ==== Restore shift/eigenvalue array from T ====
463 *
464  DO 40 i = infqr + 1, jw
465  sh( kwtop+i-1 ) = t( i, i )
466  40 CONTINUE
467 *
468 *
469  IF( ns.LT.jw .OR. s.EQ.zero ) THEN
470  IF( ns.GT.1 .AND. s.NE.zero ) THEN
471 *
472 * ==== Reflect spike back into lower triangle ====
473 *
474  CALL zcopy( ns, v, ldv, work, 1 )
475  DO 50 i = 1, ns
476  work( i ) = dconjg( work( i ) )
477  50 CONTINUE
478  beta = work( 1 )
479  CALL zlarfg( ns, beta, work( 2 ), 1, tau )
480  work( 1 ) = one
481 *
482  CALL zlaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
483 *
484  CALL zlarf( 'L', ns, jw, work, 1, dconjg( tau ), t, ldt,
485  \$ work( jw+1 ) )
486  CALL zlarf( 'R', ns, ns, work, 1, tau, t, ldt,
487  \$ work( jw+1 ) )
488  CALL zlarf( 'R', jw, ns, work, 1, tau, v, ldv,
489  \$ work( jw+1 ) )
490 *
491  CALL zgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
492  \$ lwork-jw, info )
493  END IF
494 *
495 * ==== Copy updated reduced window into place ====
496 *
497  IF( kwtop.GT.1 )
498  \$ h( kwtop, kwtop-1 ) = s*dconjg( v( 1, 1 ) )
499  CALL zlacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
500  CALL zcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
501  \$ ldh+1 )
502 *
503 * ==== Accumulate orthogonal matrix in order update
504 * . H and Z, if requested. ====
505 *
506  IF( ns.GT.1 .AND. s.NE.zero )
507  \$ CALL zunmhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
508  \$ work( jw+1 ), lwork-jw, info )
509 *
510 * ==== Update vertical slab in H ====
511 *
512  IF( wantt ) THEN
513  ltop = 1
514  ELSE
515  ltop = ktop
516  END IF
517  DO 60 krow = ltop, kwtop - 1, nv
518  kln = min( nv, kwtop-krow )
519  CALL zgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
520  \$ ldh, v, ldv, zero, wv, ldwv )
521  CALL zlacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
522  60 CONTINUE
523 *
524 * ==== Update horizontal slab in H ====
525 *
526  IF( wantt ) THEN
527  DO 70 kcol = kbot + 1, n, nh
528  kln = min( nh, n-kcol+1 )
529  CALL zgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
530  \$ h( kwtop, kcol ), ldh, zero, t, ldt )
531  CALL zlacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
532  \$ ldh )
533  70 CONTINUE
534  END IF
535 *
536 * ==== Update vertical slab in Z ====
537 *
538  IF( wantz ) THEN
539  DO 80 krow = iloz, ihiz, nv
540  kln = min( nv, ihiz-krow+1 )
541  CALL zgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
542  \$ ldz, v, ldv, zero, wv, ldwv )
543  CALL zlacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
544  \$ ldz )
545  80 CONTINUE
546  END IF
547  END IF
548 *
549 * ==== Return the number of deflations ... ====
550 *
551  nd = jw - ns
552 *
553 * ==== ... and the number of shifts. (Subtracting
554 * . INFQR from the spike length takes care
555 * . of the case of a rare QR failure while
556 * . calculating eigenvalues of the deflation
557 * . window.) ====
558 *
559  ns = ns - infqr
560 *
561 * ==== Return optimal workspace. ====
562 *
563  work( 1 ) = dcmplx( lwkopt, 0 )
564 *
565 * ==== End of ZLAQR2 ====
566 *
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMHR
Definition: zunmhr.f:180
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108