LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zlaqr2()

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).

Download ZLAQR2 + dependencies [TGZ] [ZIP] [TXT]

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,ILOZ:IHIZ) 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
          The leading dimension of V just as declared in the
          calling subroutine.  NW .LE. LDV
[in]NH
          NH is INTEGER
          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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2017
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.7.1) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * June 2017
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 *
440 * ==== Return to Hessenberg form ====
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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
subroutine zlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
Definition: zlahqr.f:197
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
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 zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
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
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine ztrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
ZTREXC
Definition: ztrexc.f:128
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108
subroutine zunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMHR
Definition: zunmhr.f:180
Here is the call graph for this function:
Here is the caller graph for this function: