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

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

Definition at line 275 of file dlaqr2.f.

278*
279* -- LAPACK auxiliary routine --
280* -- LAPACK is a software package provided by Univ. of Tennessee, --
281* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
282*
283* .. Scalar Arguments ..
284 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
285 $ LDZ, LWORK, N, ND, NH, NS, NV, NW
286 LOGICAL WANTT, WANTZ
287* ..
288* .. Array Arguments ..
289 DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
290 $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
291 $ Z( LDZ, * )
292* ..
293*
294* ================================================================
295* .. Parameters ..
296 DOUBLE PRECISION ZERO, ONE
297 parameter( zero = 0.0d0, one = 1.0d0 )
298* ..
299* .. Local Scalars ..
300 DOUBLE PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
301 $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
302 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
303 $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2,
304 $ LWKOPT
305 LOGICAL BULGE, SORTED
306* ..
307* .. External Functions ..
308 DOUBLE PRECISION DLAMCH
309 EXTERNAL dlamch
310* ..
311* .. External Subroutines ..
312 EXTERNAL dcopy, dgehrd, dgemm, dlacpy, dlahqr,
313 $ dlanv2, dlarf, dlarfg, dlaset, dormhr, dtrexc
314* ..
315* .. Intrinsic Functions ..
316 INTRINSIC abs, dble, int, max, min, sqrt
317* ..
318* .. Executable Statements ..
319*
320* ==== Estimate optimal workspace. ====
321*
322 jw = min( nw, kbot-ktop+1 )
323 IF( jw.LE.2 ) THEN
324 lwkopt = 1
325 ELSE
326*
327* ==== Workspace query call to DGEHRD ====
328*
329 CALL dgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
330 lwk1 = int( work( 1 ) )
331*
332* ==== Workspace query call to DORMHR ====
333*
334 CALL dormhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, 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 ) = dble( lwkopt )
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 = one / 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 sr( kwtop ) = h( kwtop, kwtop )
383 si( kwtop ) = zero
384 ns = 1
385 nd = 0
386 IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
387 $ THEN
388 ns = 0
389 nd = 1
390 IF( kwtop.GT.ktop )
391 $ h( kwtop, kwtop-1 ) = zero
392 END IF
393 work( 1 ) = one
394 RETURN
395 END IF
396*
397* ==== Convert to spike-triangular form. (In case of a
398* . rare QR failure, this routine continues to do
399* . aggressive early deflation using that part of
400* . the deflation window that converged using INFQR
401* . here and there to keep track.) ====
402*
403 CALL dlacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
404 CALL dcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
405*
406 CALL dlaset( 'A', jw, jw, zero, one, v, ldv )
407 CALL dlahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
408 $ si( kwtop ), 1, jw, v, ldv, infqr )
409*
410* ==== DTREXC needs a clean margin near the diagonal ====
411*
412 DO 10 j = 1, jw - 3
413 t( j+2, j ) = zero
414 t( j+3, j ) = zero
415 10 CONTINUE
416 IF( jw.GT.2 )
417 $ t( jw, jw-2 ) = zero
418*
419* ==== Deflation detection loop ====
420*
421 ns = jw
422 ilst = infqr + 1
423 20 CONTINUE
424 IF( ilst.LE.ns ) THEN
425 IF( ns.EQ.1 ) THEN
426 bulge = .false.
427 ELSE
428 bulge = t( ns, ns-1 ).NE.zero
429 END IF
430*
431* ==== Small spike tip test for deflation ====
432*
433 IF( .NOT.bulge ) THEN
434*
435* ==== Real eigenvalue ====
436*
437 foo = abs( t( ns, ns ) )
438 IF( foo.EQ.zero )
439 $ foo = abs( s )
440 IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) ) THEN
441*
442* ==== Deflatable ====
443*
444 ns = ns - 1
445 ELSE
446*
447* ==== Undeflatable. Move it up out of the way.
448* . (DTREXC can not fail in this case.) ====
449*
450 ifst = ns
451 CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
452 $ info )
453 ilst = ilst + 1
454 END IF
455 ELSE
456*
457* ==== Complex conjugate pair ====
458*
459 foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
460 $ sqrt( abs( t( ns-1, ns ) ) )
461 IF( foo.EQ.zero )
462 $ foo = abs( s )
463 IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
464 $ max( smlnum, ulp*foo ) ) THEN
465*
466* ==== Deflatable ====
467*
468 ns = ns - 2
469 ELSE
470*
471* ==== Undeflatable. Move them up out of the way.
472* . Fortunately, DTREXC does the right thing with
473* . ILST in case of a rare exchange failure. ====
474*
475 ifst = ns
476 CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
477 $ info )
478 ilst = ilst + 2
479 END IF
480 END IF
481*
482* ==== End deflation detection loop ====
483*
484 GO TO 20
485 END IF
486*
487* ==== Return to Hessenberg form ====
488*
489 IF( ns.EQ.0 )
490 $ s = zero
491*
492 IF( ns.LT.jw ) THEN
493*
494* ==== sorting diagonal blocks of T improves accuracy for
495* . graded matrices. Bubble sort deals well with
496* . exchange failures. ====
497*
498 sorted = .false.
499 i = ns + 1
500 30 CONTINUE
501 IF( sorted )
502 $ GO TO 50
503 sorted = .true.
504*
505 kend = i - 1
506 i = infqr + 1
507 IF( i.EQ.ns ) THEN
508 k = i + 1
509 ELSE IF( t( i+1, i ).EQ.zero ) THEN
510 k = i + 1
511 ELSE
512 k = i + 2
513 END IF
514 40 CONTINUE
515 IF( k.LE.kend ) THEN
516 IF( k.EQ.i+1 ) THEN
517 evi = abs( t( i, i ) )
518 ELSE
519 evi = abs( t( i, i ) ) + sqrt( abs( t( i+1, i ) ) )*
520 $ sqrt( abs( t( i, i+1 ) ) )
521 END IF
522*
523 IF( k.EQ.kend ) THEN
524 evk = abs( t( k, k ) )
525 ELSE IF( t( k+1, k ).EQ.zero ) THEN
526 evk = abs( t( k, k ) )
527 ELSE
528 evk = abs( t( k, k ) ) + sqrt( abs( t( k+1, k ) ) )*
529 $ sqrt( abs( t( k, k+1 ) ) )
530 END IF
531*
532 IF( evi.GE.evk ) THEN
533 i = k
534 ELSE
535 sorted = .false.
536 ifst = i
537 ilst = k
538 CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
539 $ info )
540 IF( info.EQ.0 ) THEN
541 i = ilst
542 ELSE
543 i = k
544 END IF
545 END IF
546 IF( i.EQ.kend ) THEN
547 k = i + 1
548 ELSE IF( t( i+1, i ).EQ.zero ) THEN
549 k = i + 1
550 ELSE
551 k = i + 2
552 END IF
553 GO TO 40
554 END IF
555 GO TO 30
556 50 CONTINUE
557 END IF
558*
559* ==== Restore shift/eigenvalue array from T ====
560*
561 i = jw
562 60 CONTINUE
563 IF( i.GE.infqr+1 ) THEN
564 IF( i.EQ.infqr+1 ) THEN
565 sr( kwtop+i-1 ) = t( i, i )
566 si( kwtop+i-1 ) = zero
567 i = i - 1
568 ELSE IF( t( i, i-1 ).EQ.zero ) THEN
569 sr( kwtop+i-1 ) = t( i, i )
570 si( kwtop+i-1 ) = zero
571 i = i - 1
572 ELSE
573 aa = t( i-1, i-1 )
574 cc = t( i, i-1 )
575 bb = t( i-1, i )
576 dd = t( i, i )
577 CALL dlanv2( aa, bb, cc, dd, sr( kwtop+i-2 ),
578 $ si( kwtop+i-2 ), sr( kwtop+i-1 ),
579 $ si( kwtop+i-1 ), cs, sn )
580 i = i - 2
581 END IF
582 GO TO 60
583 END IF
584*
585 IF( ns.LT.jw .OR. s.EQ.zero ) THEN
586 IF( ns.GT.1 .AND. s.NE.zero ) THEN
587*
588* ==== Reflect spike back into lower triangle ====
589*
590 CALL dcopy( ns, v, ldv, work, 1 )
591 beta = work( 1 )
592 CALL dlarfg( ns, beta, work( 2 ), 1, tau )
593 work( 1 ) = one
594*
595 CALL dlaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
596*
597 CALL dlarf( 'L', ns, jw, work, 1, tau, t, ldt,
598 $ work( jw+1 ) )
599 CALL dlarf( 'R', ns, ns, work, 1, tau, t, ldt,
600 $ work( jw+1 ) )
601 CALL dlarf( 'R', jw, ns, work, 1, tau, v, ldv,
602 $ work( jw+1 ) )
603*
604 CALL dgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
605 $ lwork-jw, info )
606 END IF
607*
608* ==== Copy updated reduced window into place ====
609*
610 IF( kwtop.GT.1 )
611 $ h( kwtop, kwtop-1 ) = s*v( 1, 1 )
612 CALL dlacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
613 CALL dcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
614 $ ldh+1 )
615*
616* ==== Accumulate orthogonal matrix in order update
617* . H and Z, if requested. ====
618*
619 IF( ns.GT.1 .AND. s.NE.zero )
620 $ CALL dormhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
621 $ work( jw+1 ), lwork-jw, info )
622*
623* ==== Update vertical slab in H ====
624*
625 IF( wantt ) THEN
626 ltop = 1
627 ELSE
628 ltop = ktop
629 END IF
630 DO 70 krow = ltop, kwtop - 1, nv
631 kln = min( nv, kwtop-krow )
632 CALL dgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
633 $ ldh, v, ldv, zero, wv, ldwv )
634 CALL dlacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
635 70 CONTINUE
636*
637* ==== Update horizontal slab in H ====
638*
639 IF( wantt ) THEN
640 DO 80 kcol = kbot + 1, n, nh
641 kln = min( nh, n-kcol+1 )
642 CALL dgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
643 $ h( kwtop, kcol ), ldh, zero, t, ldt )
644 CALL dlacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
645 $ ldh )
646 80 CONTINUE
647 END IF
648*
649* ==== Update vertical slab in Z ====
650*
651 IF( wantz ) THEN
652 DO 90 krow = iloz, ihiz, nv
653 kln = min( nv, ihiz-krow+1 )
654 CALL dgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
655 $ ldz, v, ldv, zero, wv, ldwv )
656 CALL dlacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
657 $ ldz )
658 90 CONTINUE
659 END IF
660 END IF
661*
662* ==== Return the number of deflations ... ====
663*
664 nd = jw - ns
665*
666* ==== ... and the number of shifts. (Subtracting
667* . INFQR from the spike length takes care
668* . of the case of a rare QR failure while
669* . calculating eigenvalues of the deflation
670* . window.) ====
671*
672 ns = ns - infqr
673*
674* ==== Return optimal workspace. ====
675*
676 work( 1 ) = dble( lwkopt )
677*
678* ==== End of DLAQR2 ====
679*
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgehrd
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:167
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlahqr
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,...
Definition dlahqr.f:207
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanv2
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:127
dlarf
subroutine dlarf(side, m, n, v, incv, tau, c, ldc, work)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition dlarf.f:124
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
dlaset
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:110
dtrexc
subroutine dtrexc(compq, n, t, ldt, q, ldq, ifst, ilst, work, info)
DTREXC
Definition dtrexc.f:148
dormhr
subroutine dormhr(side, trans, m, n, ilo, ihi, a, lda, tau, c, ldc, work, lwork, info)
DORMHR
Definition dormhr.f:178
Here is the call graph for this function:
Here is the caller graph for this function: