LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ zlaqr3()

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 <= NW <= (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 <= 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 <= ILOZ <= IHIZ <= 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 <= 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 <= LDV
[in]NH
          NH is INTEGER
          The number of columns of T.  NH >= 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 <= LDT
[in]NV
          NV is INTEGER
          The number of rows of work array WV available for
          workspace.  NV >= 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 <= 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.
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

Definition at line 264 of file zlaqr3.f.

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