LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zlaqr3 ( 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 
)

ZLAQR3 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 ZLAQR3 + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
    Aggressive early deflation:

    ZLAQR3 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 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; ZLAQR3
          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 2016
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

Definition at line 269 of file zlaqr3.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: