LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages

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

Definition at line 265 of file zlaqr2.f.

269*
270* -- LAPACK auxiliary routine --
271* -- LAPACK is a software package provided by Univ. of Tennessee, --
272* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
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*16 H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
281 $ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
282* ..
283*
284* ================================================================
285*
286* .. Parameters ..
287 COMPLEX*16 ZERO, ONE
288 parameter( zero = ( 0.0d0, 0.0d0 ),
289 $ one = ( 1.0d0, 0.0d0 ) )
290 DOUBLE PRECISION RZERO, RONE
291 parameter( rzero = 0.0d0, rone = 1.0d0 )
292* ..
293* .. Local Scalars ..
294 COMPLEX*16 CDUM, S, TAU
295 DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM, ULP
296 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
297 $ KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWKOPT
298* ..
299* .. External Functions ..
300 DOUBLE PRECISION DLAMCH
301 EXTERNAL dlamch
302* ..
303* .. External Subroutines ..
304 EXTERNAL zcopy, zgehrd, zgemm, zlacpy,
305 $ zlahqr,
306 $ zlarf1f, zlarfg, zlaset, ztrexc, zunmhr
307* ..
308* .. Intrinsic Functions ..
309 INTRINSIC abs, dble, dcmplx, dconjg, dimag, int, max, min
310* ..
311* .. Statement Functions ..
312 DOUBLE PRECISION CABS1
313* ..
314* .. Statement Function definitions ..
315 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
316* ..
317* .. Executable Statements ..
318*
319* ==== Estimate optimal workspace. ====
320*
321 jw = min( nw, kbot-ktop+1 )
322 IF( jw.LE.2 ) THEN
323 lwkopt = 1
324 ELSE
325*
326* ==== Workspace query call to ZGEHRD ====
327*
328 CALL zgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
329 lwk1 = int( work( 1 ) )
330*
331* ==== Workspace query call to ZUNMHR ====
332*
333 CALL zunmhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v,
334 $ ldv,
335 $ work, -1, info )
336 lwk2 = int( work( 1 ) )
337*
338* ==== Optimal workspace ====
339*
340 lwkopt = jw + max( lwk1, lwk2 )
341 END IF
342*
343* ==== Quick return in case of workspace query. ====
344*
345 IF( lwork.EQ.-1 ) THEN
346 work( 1 ) = dcmplx( lwkopt, 0 )
347 RETURN
348 END IF
349*
350* ==== Nothing to do ...
351* ... for an empty active block ... ====
352 ns = 0
353 nd = 0
354 work( 1 ) = one
355 IF( ktop.GT.kbot )
356 $ RETURN
357* ... nor for an empty deflation window. ====
358 IF( nw.LT.1 )
359 $ RETURN
360*
361* ==== Machine constants ====
362*
363 safmin = dlamch( 'SAFE MINIMUM' )
364 safmax = rone / safmin
365 ulp = dlamch( 'PRECISION' )
366 smlnum = safmin*( dble( n ) / ulp )
367*
368* ==== Setup deflation window ====
369*
370 jw = min( nw, kbot-ktop+1 )
371 kwtop = kbot - jw + 1
372 IF( kwtop.EQ.ktop ) THEN
373 s = zero
374 ELSE
375 s = h( kwtop, kwtop-1 )
376 END IF
377*
378 IF( kbot.EQ.kwtop ) THEN
379*
380* ==== 1-by-1 deflation window: not much to do ====
381*
382 sh( kwtop ) = h( kwtop, kwtop )
383 ns = 1
384 nd = 0
385 IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
386 $ kwtop ) ) ) ) THEN
387 ns = 0
388 nd = 1
389 IF( kwtop.GT.ktop )
390 $ h( kwtop, kwtop-1 ) = zero
391 END IF
392 work( 1 ) = one
393 RETURN
394 END IF
395*
396* ==== Convert to spike-triangular form. (In case of a
397* . rare QR failure, this routine continues to do
398* . aggressive early deflation using that part of
399* . the deflation window that converged using INFQR
400* . here and there to keep track.) ====
401*
402 CALL zlacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
403 CALL zcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ),
404 $ ldt+1 )
405*
406 CALL zlaset( 'A', jw, jw, zero, one, v, ldv )
407 CALL zlahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
408 $ jw, v, ldv, infqr )
409*
410* ==== Deflation detection loop ====
411*
412 ns = jw
413 ilst = infqr + 1
414 DO 10 knt = infqr + 1, jw
415*
416* ==== Small spike tip deflation test ====
417*
418 foo = cabs1( t( ns, ns ) )
419 IF( foo.EQ.rzero )
420 $ foo = cabs1( s )
421 IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
422 $ THEN
423*
424* ==== One more converged eigenvalue ====
425*
426 ns = ns - 1
427 ELSE
428*
429* ==== One undeflatable eigenvalue. Move it up out of the
430* . way. (ZTREXC can not fail in this case.) ====
431*
432 ifst = ns
433 CALL ztrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
434 ilst = ilst + 1
435 END IF
436 10 CONTINUE
437*
438* ==== Return to Hessenberg form ====
439*
440 IF( ns.EQ.0 )
441 $ s = zero
442*
443 IF( ns.LT.jw ) THEN
444*
445* ==== sorting the diagonal of T improves accuracy for
446* . graded matrices. ====
447*
448 DO 30 i = infqr + 1, ns
449 ifst = i
450 DO 20 j = i + 1, ns
451 IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
452 $ ifst = j
453 20 CONTINUE
454 ilst = i
455 IF( ifst.NE.ilst )
456 $ CALL ztrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst,
457 $ info )
458 30 CONTINUE
459 END IF
460*
461* ==== Restore shift/eigenvalue array from T ====
462*
463 DO 40 i = infqr + 1, jw
464 sh( kwtop+i-1 ) = t( i, i )
465 40 CONTINUE
466*
467*
468 IF( ns.LT.jw .OR. s.EQ.zero ) THEN
469 IF( ns.GT.1 .AND. s.NE.zero ) THEN
470*
471* ==== Reflect spike back into lower triangle ====
472*
473 CALL zcopy( ns, v, ldv, work, 1 )
474 DO 50 i = 1, ns
475 work( i ) = dconjg( work( i ) )
476 50 CONTINUE
477 CALL zlarfg( ns, work( 1 ), work( 2 ), 1, tau )
478*
479 CALL zlaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ),
480 $ ldt )
481*
482 CALL zlarf1f( 'L', ns, jw, work, 1, conjg( tau ), t, ldt,
483 $ work( jw+1 ) )
484 CALL zlarf1f( 'R', ns, ns, work, 1, tau, t, ldt,
485 $ work( jw+1 ) )
486 CALL zlarf1f( 'R', jw, ns, work, 1, tau, v, ldv,
487 $ work( jw+1 ) )
488*
489 CALL zgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
490 $ lwork-jw, info )
491 END IF
492*
493* ==== Copy updated reduced window into place ====
494*
495 IF( kwtop.GT.1 )
496 $ h( kwtop, kwtop-1 ) = s*dconjg( v( 1, 1 ) )
497 CALL zlacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
498 CALL zcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
499 $ ldh+1 )
500*
501* ==== Accumulate orthogonal matrix in order update
502* . H and Z, if requested. ====
503*
504 IF( ns.GT.1 .AND. s.NE.zero )
505 $ CALL zunmhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v,
506 $ ldv,
507 $ work( jw+1 ), lwork-jw, info )
508*
509* ==== Update vertical slab in H ====
510*
511 IF( wantt ) THEN
512 ltop = 1
513 ELSE
514 ltop = ktop
515 END IF
516 DO 60 krow = ltop, kwtop - 1, nv
517 kln = min( nv, kwtop-krow )
518 CALL zgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
519 $ ldh, v, ldv, zero, wv, ldwv )
520 CALL zlacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ),
521 $ 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,
542 $ kwtop ),
543 $ ldz, v, ldv, zero, wv, ldwv )
544 CALL zlacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
545 $ ldz )
546 80 CONTINUE
547 END IF
548 END IF
549*
550* ==== Return the number of deflations ... ====
551*
552 nd = jw - ns
553*
554* ==== ... and the number of shifts. (Subtracting
555* . INFQR from the spike length takes care
556* . of the case of a rare QR failure while
557* . calculating eigenvalues of the deflation
558* . window.) ====
559*
560 ns = ns - infqr
561*
562* ==== Return optimal workspace. ====
563*
564 work( 1 ) = dcmplx( lwkopt, 0 )
565*
566* ==== End of ZLAQR2 ====
567*
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgehrd
subroutine zgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZGEHRD
Definition zgehrd.f:166
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
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:101
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,...
Definition zlahqr.f:193
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlarf1f
subroutine zlarf1f(side, m, n, v, incv, tau, c, ldc, work)
ZLARF1F applies an elementary reflector to a general rectangular
Definition zlarf1f.f:157
zlarfg
subroutine zlarfg(n, alpha, x, incx, tau)
ZLARFG generates an elementary reflector (Householder matrix).
Definition zlarfg.f:104
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:104
ztrexc
subroutine ztrexc(compq, n, t, ldt, q, ldq, ifst, ilst, info)
ZTREXC
Definition ztrexc.f:124
zunmhr
subroutine zunmhr(side, trans, m, n, ilo, ihi, a, lda, tau, c, ldc, work, lwork, info)
ZUNMHR
Definition zunmhr.f:176
Here is the call graph for this function:
Here is the caller graph for this function: