LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 262 of file zlaqr3.f.

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