LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dlaqr2()

subroutine dlaqr2 ( logical  WANTT,
logical  WANTZ,
integer  N,
integer  KTOP,
integer  KBOT,
integer  NW,
double precision, dimension( ldh, * )  H,
integer  LDH,
integer  ILOZ,
integer  IHIZ,
double precision, dimension( ldz, * )  Z,
integer  LDZ,
integer  NS,
integer  ND,
double precision, dimension( * )  SR,
double precision, dimension( * )  SI,
double precision, dimension( ldv, * )  V,
integer  LDV,
integer  NH,
double precision, dimension( ldt, * )  T,
integer  LDT,
integer  NV,
double precision, dimension( ldwv, * )  WV,
integer  LDWV,
double precision, dimension( * )  WORK,
integer  LWORK 
)

DLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).

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

Purpose:
    DLAQR2 is identical to DLAQR3 except that it avoids
    recursion by calling DLAHQR instead of DLAQR4.

    Aggressive early deflation:

    This subroutine accepts as input an upper Hessenberg matrix
    H and performs an orthogonal 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 orthogonal 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 quasi-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 orthogonal matrix Z is updated so
          so that the orthogonal 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 orthogonal 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 DOUBLE PRECISION 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 an orthogonal
          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 DOUBLE PRECISION array, dimension (LDZ,N)
          IF WANTZ is .TRUE., then on output, the orthogonal
          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]SR
          SR is DOUBLE PRECISION array, dimension (KBOT)
[out]SI
          SI is DOUBLE PRECISION array, dimension (KBOT)
          On output, the real and imaginary parts of approximate
          eigenvalues that may be used for shifts are stored in
          SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
          SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
          The real and imaginary parts of converged eigenvalues
          are stored in SR(KBOT-ND+1) through SR(KBOT) and
          SI(KBOT-ND+1) through SI(KBOT), respectively.
[out]V
          V is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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; DLAQR2
          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 280 of file dlaqr2.f.

280 *
281 * -- LAPACK auxiliary routine (version 3.7.1) --
282 * -- LAPACK is a software package provided by Univ. of Tennessee, --
283 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
284 * June 2017
285 *
286 * .. Scalar Arguments ..
287  INTEGER ihiz, iloz, kbot, ktop, ldh, ldt, ldv, ldwv,
288  $ ldz, lwork, n, nd, nh, ns, nv, nw
289  LOGICAL wantt, wantz
290 * ..
291 * .. Array Arguments ..
292  DOUBLE PRECISION h( ldh, * ), si( * ), sr( * ), t( ldt, * ),
293  $ v( ldv, * ), work( * ), wv( ldwv, * ),
294  $ z( ldz, * )
295 * ..
296 *
297 * ================================================================
298 * .. Parameters ..
299  DOUBLE PRECISION zero, one
300  parameter( zero = 0.0d0, one = 1.0d0 )
301 * ..
302 * .. Local Scalars ..
303  DOUBLE PRECISION aa, bb, beta, cc, cs, dd, evi, evk, foo, s,
304  $ safmax, safmin, smlnum, sn, tau, ulp
305  INTEGER i, ifst, ilst, info, infqr, j, jw, k, kcol,
306  $ kend, kln, krow, kwtop, ltop, lwk1, lwk2,
307  $ lwkopt
308  LOGICAL bulge, sorted
309 * ..
310 * .. External Functions ..
311  DOUBLE PRECISION dlamch
312  EXTERNAL dlamch
313 * ..
314 * .. External Subroutines ..
315  EXTERNAL dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr,
317 * ..
318 * .. Intrinsic Functions ..
319  INTRINSIC abs, dble, int, max, min, sqrt
320 * ..
321 * .. Executable Statements ..
322 *
323 * ==== Estimate optimal workspace. ====
324 *
325  jw = min( nw, kbot-ktop+1 )
326  IF( jw.LE.2 ) THEN
327  lwkopt = 1
328  ELSE
329 *
330 * ==== Workspace query call to DGEHRD ====
331 *
332  CALL dgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
333  lwk1 = int( work( 1 ) )
334 *
335 * ==== Workspace query call to DORMHR ====
336 *
337  CALL dormhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
338  $ work, -1, info )
339  lwk2 = int( work( 1 ) )
340 *
341 * ==== Optimal workspace ====
342 *
343  lwkopt = jw + max( lwk1, lwk2 )
344  END IF
345 *
346 * ==== Quick return in case of workspace query. ====
347 *
348  IF( lwork.EQ.-1 ) THEN
349  work( 1 ) = dble( lwkopt )
350  RETURN
351  END IF
352 *
353 * ==== Nothing to do ...
354 * ... for an empty active block ... ====
355  ns = 0
356  nd = 0
357  work( 1 ) = one
358  IF( ktop.GT.kbot )
359  $ RETURN
360 * ... nor for an empty deflation window. ====
361  IF( nw.LT.1 )
362  $ RETURN
363 *
364 * ==== Machine constants ====
365 *
366  safmin = dlamch( 'SAFE MINIMUM' )
367  safmax = one / safmin
368  CALL dlabad( safmin, safmax )
369  ulp = dlamch( 'PRECISION' )
370  smlnum = safmin*( dble( n ) / ulp )
371 *
372 * ==== Setup deflation window ====
373 *
374  jw = min( nw, kbot-ktop+1 )
375  kwtop = kbot - jw + 1
376  IF( kwtop.EQ.ktop ) THEN
377  s = zero
378  ELSE
379  s = h( kwtop, kwtop-1 )
380  END IF
381 *
382  IF( kbot.EQ.kwtop ) THEN
383 *
384 * ==== 1-by-1 deflation window: not much to do ====
385 *
386  sr( kwtop ) = h( kwtop, kwtop )
387  si( kwtop ) = zero
388  ns = 1
389  nd = 0
390  IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
391  $ 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 dlacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
408  CALL dcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
409 *
410  CALL dlaset( 'A', jw, jw, zero, one, v, ldv )
411  CALL dlahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
412  $ si( kwtop ), 1, jw, v, ldv, infqr )
413 *
414 * ==== DTREXC needs a clean margin near the diagonal ====
415 *
416  DO 10 j = 1, jw - 3
417  t( j+2, j ) = zero
418  t( j+3, j ) = zero
419  10 CONTINUE
420  IF( jw.GT.2 )
421  $ t( jw, jw-2 ) = zero
422 *
423 * ==== Deflation detection loop ====
424 *
425  ns = jw
426  ilst = infqr + 1
427  20 CONTINUE
428  IF( ilst.LE.ns ) THEN
429  IF( ns.EQ.1 ) THEN
430  bulge = .false.
431  ELSE
432  bulge = t( ns, ns-1 ).NE.zero
433  END IF
434 *
435 * ==== Small spike tip test for deflation ====
436 *
437  IF( .NOT.bulge ) THEN
438 *
439 * ==== Real eigenvalue ====
440 *
441  foo = abs( t( ns, ns ) )
442  IF( foo.EQ.zero )
443  $ foo = abs( s )
444  IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) ) THEN
445 *
446 * ==== Deflatable ====
447 *
448  ns = ns - 1
449  ELSE
450 *
451 * ==== Undeflatable. Move it up out of the way.
452 * . (DTREXC can not fail in this case.) ====
453 *
454  ifst = ns
455  CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
456  $ info )
457  ilst = ilst + 1
458  END IF
459  ELSE
460 *
461 * ==== Complex conjugate pair ====
462 *
463  foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
464  $ sqrt( abs( t( ns-1, ns ) ) )
465  IF( foo.EQ.zero )
466  $ foo = abs( s )
467  IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
468  $ max( smlnum, ulp*foo ) ) THEN
469 *
470 * ==== Deflatable ====
471 *
472  ns = ns - 2
473  ELSE
474 *
475 * ==== Undeflatable. Move them up out of the way.
476 * . Fortunately, DTREXC does the right thing with
477 * . ILST in case of a rare exchange failure. ====
478 *
479  ifst = ns
480  CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
481  $ info )
482  ilst = ilst + 2
483  END IF
484  END IF
485 *
486 * ==== End deflation detection loop ====
487 *
488  GO TO 20
489  END IF
490 *
491 * ==== Return to Hessenberg form ====
492 *
493  IF( ns.EQ.0 )
494  $ s = zero
495 *
496  IF( ns.LT.jw ) THEN
497 *
498 * ==== sorting diagonal blocks of T improves accuracy for
499 * . graded matrices. Bubble sort deals well with
500 * . exchange failures. ====
501 *
502  sorted = .false.
503  i = ns + 1
504  30 CONTINUE
505  IF( sorted )
506  $ GO TO 50
507  sorted = .true.
508 *
509  kend = i - 1
510  i = infqr + 1
511  IF( i.EQ.ns ) THEN
512  k = i + 1
513  ELSE IF( t( i+1, i ).EQ.zero ) THEN
514  k = i + 1
515  ELSE
516  k = i + 2
517  END IF
518  40 CONTINUE
519  IF( k.LE.kend ) THEN
520  IF( k.EQ.i+1 ) THEN
521  evi = abs( t( i, i ) )
522  ELSE
523  evi = abs( t( i, i ) ) + sqrt( abs( t( i+1, i ) ) )*
524  $ sqrt( abs( t( i, i+1 ) ) )
525  END IF
526 *
527  IF( k.EQ.kend ) THEN
528  evk = abs( t( k, k ) )
529  ELSE IF( t( k+1, k ).EQ.zero ) THEN
530  evk = abs( t( k, k ) )
531  ELSE
532  evk = abs( t( k, k ) ) + sqrt( abs( t( k+1, k ) ) )*
533  $ sqrt( abs( t( k, k+1 ) ) )
534  END IF
535 *
536  IF( evi.GE.evk ) THEN
537  i = k
538  ELSE
539  sorted = .false.
540  ifst = i
541  ilst = k
542  CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
543  $ info )
544  IF( info.EQ.0 ) THEN
545  i = ilst
546  ELSE
547  i = k
548  END IF
549  END IF
550  IF( i.EQ.kend ) THEN
551  k = i + 1
552  ELSE IF( t( i+1, i ).EQ.zero ) THEN
553  k = i + 1
554  ELSE
555  k = i + 2
556  END IF
557  GO TO 40
558  END IF
559  GO TO 30
560  50 CONTINUE
561  END IF
562 *
563 * ==== Restore shift/eigenvalue array from T ====
564 *
565  i = jw
566  60 CONTINUE
567  IF( i.GE.infqr+1 ) THEN
568  IF( i.EQ.infqr+1 ) THEN
569  sr( kwtop+i-1 ) = t( i, i )
570  si( kwtop+i-1 ) = zero
571  i = i - 1
572  ELSE IF( t( i, i-1 ).EQ.zero ) THEN
573  sr( kwtop+i-1 ) = t( i, i )
574  si( kwtop+i-1 ) = zero
575  i = i - 1
576  ELSE
577  aa = t( i-1, i-1 )
578  cc = t( i, i-1 )
579  bb = t( i-1, i )
580  dd = t( i, i )
581  CALL dlanv2( aa, bb, cc, dd, sr( kwtop+i-2 ),
582  $ si( kwtop+i-2 ), sr( kwtop+i-1 ),
583  $ si( kwtop+i-1 ), cs, sn )
584  i = i - 2
585  END IF
586  GO TO 60
587  END IF
588 *
589  IF( ns.LT.jw .OR. s.EQ.zero ) THEN
590  IF( ns.GT.1 .AND. s.NE.zero ) THEN
591 *
592 * ==== Reflect spike back into lower triangle ====
593 *
594  CALL dcopy( ns, v, ldv, work, 1 )
595  beta = work( 1 )
596  CALL dlarfg( ns, beta, work( 2 ), 1, tau )
597  work( 1 ) = one
598 *
599  CALL dlaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
600 *
601  CALL dlarf( 'L', ns, jw, work, 1, tau, t, ldt,
602  $ work( jw+1 ) )
603  CALL dlarf( 'R', ns, ns, work, 1, tau, t, ldt,
604  $ work( jw+1 ) )
605  CALL dlarf( 'R', jw, ns, work, 1, tau, v, ldv,
606  $ work( jw+1 ) )
607 *
608  CALL dgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
609  $ lwork-jw, info )
610  END IF
611 *
612 * ==== Copy updated reduced window into place ====
613 *
614  IF( kwtop.GT.1 )
615  $ h( kwtop, kwtop-1 ) = s*v( 1, 1 )
616  CALL dlacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
617  CALL dcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
618  $ ldh+1 )
619 *
620 * ==== Accumulate orthogonal matrix in order update
621 * . H and Z, if requested. ====
622 *
623  IF( ns.GT.1 .AND. s.NE.zero )
624  $ CALL dormhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
625  $ work( jw+1 ), lwork-jw, info )
626 *
627 * ==== Update vertical slab in H ====
628 *
629  IF( wantt ) THEN
630  ltop = 1
631  ELSE
632  ltop = ktop
633  END IF
634  DO 70 krow = ltop, kwtop - 1, nv
635  kln = min( nv, kwtop-krow )
636  CALL dgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
637  $ ldh, v, ldv, zero, wv, ldwv )
638  CALL dlacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
639  70 CONTINUE
640 *
641 * ==== Update horizontal slab in H ====
642 *
643  IF( wantt ) THEN
644  DO 80 kcol = kbot + 1, n, nh
645  kln = min( nh, n-kcol+1 )
646  CALL dgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
647  $ h( kwtop, kcol ), ldh, zero, t, ldt )
648  CALL dlacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
649  $ ldh )
650  80 CONTINUE
651  END IF
652 *
653 * ==== Update vertical slab in Z ====
654 *
655  IF( wantz ) THEN
656  DO 90 krow = iloz, ihiz, nv
657  kln = min( nv, ihiz-krow+1 )
658  CALL dgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
659  $ ldz, v, ldv, zero, wv, ldwv )
660  CALL dlacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
661  $ ldz )
662  90 CONTINUE
663  END IF
664  END IF
665 *
666 * ==== Return the number of deflations ... ====
667 *
668  nd = jw - ns
669 *
670 * ==== ... and the number of shifts. (Subtracting
671 * . INFQR from the spike length takes care
672 * . of the case of a rare QR failure while
673 * . calculating eigenvalues of the deflation
674 * . window.) ====
675 *
676  ns = ns - infqr
677 *
678 * ==== Return optimal workspace. ====
679 *
680  work( 1 ) = dble( lwkopt )
681 *
682 * ==== End of DLAQR2 ====
683 *
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:169
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
subroutine dlanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form...
Definition: dlanv2.f:129
subroutine dlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
Definition: dlahqr.f:209
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dtrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
DTREXC
Definition: dtrexc.f:150
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMHR
Definition: dormhr.f:180
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126
Here is the call graph for this function:
Here is the caller graph for this function: