LAPACK 3.12.1
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 273 of file dlaqr2.f.

277*
278* -- LAPACK auxiliary routine --
279* -- LAPACK is a software package provided by Univ. of Tennessee, --
280* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
281*
282* .. Scalar Arguments ..
283 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
284 $ LDZ, LWORK, N, ND, NH, NS, NV, NW
285 LOGICAL WANTT, WANTZ
286* ..
287* .. Array Arguments ..
288 DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
289 $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
290 $ Z( LDZ, * )
291* ..
292*
293* ================================================================
294* .. Parameters ..
295 DOUBLE PRECISION ZERO, ONE
296 parameter( zero = 0.0d0, one = 1.0d0 )
297* ..
298* .. Local Scalars ..
299 DOUBLE PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
300 $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
301 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
302 $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2,
303 $ LWKOPT
304 LOGICAL BULGE, SORTED
305* ..
306* .. External Functions ..
307 DOUBLE PRECISION DLAMCH
308 EXTERNAL dlamch
309* ..
310* .. External Subroutines ..
311 EXTERNAL dcopy, dgehrd, dgemm, dlacpy,
312 $ dlahqr,
313 $ dlanv2, dlarf1f, 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,
335 $ ldv,
336 $ work, -1, info )
337 lwk2 = int( work( 1 ) )
338*
339* ==== Optimal workspace ====
340*
341 lwkopt = jw + max( lwk1, lwk2 )
342 END IF
343*
344* ==== Quick return in case of workspace query. ====
345*
346 IF( lwork.EQ.-1 ) THEN
347 work( 1 ) = dble( lwkopt )
348 RETURN
349 END IF
350*
351* ==== Nothing to do ...
352* ... for an empty active block ... ====
353 ns = 0
354 nd = 0
355 work( 1 ) = one
356 IF( ktop.GT.kbot )
357 $ RETURN
358* ... nor for an empty deflation window. ====
359 IF( nw.LT.1 )
360 $ RETURN
361*
362* ==== Machine constants ====
363*
364 safmin = dlamch( 'SAFE MINIMUM' )
365 safmax = one / safmin
366 ulp = dlamch( 'PRECISION' )
367 smlnum = safmin*( dble( n ) / ulp )
368*
369* ==== Setup deflation window ====
370*
371 jw = min( nw, kbot-ktop+1 )
372 kwtop = kbot - jw + 1
373 IF( kwtop.EQ.ktop ) THEN
374 s = zero
375 ELSE
376 s = h( kwtop, kwtop-1 )
377 END IF
378*
379 IF( kbot.EQ.kwtop ) THEN
380*
381* ==== 1-by-1 deflation window: not much to do ====
382*
383 sr( kwtop ) = h( kwtop, kwtop )
384 si( kwtop ) = zero
385 ns = 1
386 nd = 0
387 IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
388 $ THEN
389 ns = 0
390 nd = 1
391 IF( kwtop.GT.ktop )
392 $ h( kwtop, kwtop-1 ) = zero
393 END IF
394 work( 1 ) = one
395 RETURN
396 END IF
397*
398* ==== Convert to spike-triangular form. (In case of a
399* . rare QR failure, this routine continues to do
400* . aggressive early deflation using that part of
401* . the deflation window that converged using INFQR
402* . here and there to keep track.) ====
403*
404 CALL dlacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
405 CALL dcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ),
406 $ ldt+1 )
407*
408 CALL dlaset( 'A', jw, jw, zero, one, v, ldv )
409 CALL dlahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
410 $ si( kwtop ), 1, jw, v, ldv, infqr )
411*
412* ==== DTREXC needs a clean margin near the diagonal ====
413*
414 DO 10 j = 1, jw - 3
415 t( j+2, j ) = zero
416 t( j+3, j ) = zero
417 10 CONTINUE
418 IF( jw.GT.2 )
419 $ t( jw, jw-2 ) = zero
420*
421* ==== Deflation detection loop ====
422*
423 ns = jw
424 ilst = infqr + 1
425 20 CONTINUE
426 IF( ilst.LE.ns ) THEN
427 IF( ns.EQ.1 ) THEN
428 bulge = .false.
429 ELSE
430 bulge = t( ns, ns-1 ).NE.zero
431 END IF
432*
433* ==== Small spike tip test for deflation ====
434*
435 IF( .NOT.bulge ) THEN
436*
437* ==== Real eigenvalue ====
438*
439 foo = abs( t( ns, ns ) )
440 IF( foo.EQ.zero )
441 $ foo = abs( s )
442 IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) ) THEN
443*
444* ==== Deflatable ====
445*
446 ns = ns - 1
447 ELSE
448*
449* ==== Undeflatable. Move it up out of the way.
450* . (DTREXC can not fail in this case.) ====
451*
452 ifst = ns
453 CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst,
454 $ work,
455 $ info )
456 ilst = ilst + 1
457 END IF
458 ELSE
459*
460* ==== Complex conjugate pair ====
461*
462 foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
463 $ sqrt( abs( t( ns-1, ns ) ) )
464 IF( foo.EQ.zero )
465 $ foo = abs( s )
466 IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
467 $ max( smlnum, ulp*foo ) ) THEN
468*
469* ==== Deflatable ====
470*
471 ns = ns - 2
472 ELSE
473*
474* ==== Undeflatable. Move them up out of the way.
475* . Fortunately, DTREXC does the right thing with
476* . ILST in case of a rare exchange failure. ====
477*
478 ifst = ns
479 CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst,
480 $ 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,
543 $ work,
544 $ info )
545 IF( info.EQ.0 ) THEN
546 i = ilst
547 ELSE
548 i = k
549 END IF
550 END IF
551 IF( i.EQ.kend ) THEN
552 k = i + 1
553 ELSE IF( t( i+1, i ).EQ.zero ) THEN
554 k = i + 1
555 ELSE
556 k = i + 2
557 END IF
558 GO TO 40
559 END IF
560 GO TO 30
561 50 CONTINUE
562 END IF
563*
564* ==== Restore shift/eigenvalue array from T ====
565*
566 i = jw
567 60 CONTINUE
568 IF( i.GE.infqr+1 ) THEN
569 IF( i.EQ.infqr+1 ) THEN
570 sr( kwtop+i-1 ) = t( i, i )
571 si( kwtop+i-1 ) = zero
572 i = i - 1
573 ELSE IF( t( i, i-1 ).EQ.zero ) THEN
574 sr( kwtop+i-1 ) = t( i, i )
575 si( kwtop+i-1 ) = zero
576 i = i - 1
577 ELSE
578 aa = t( i-1, i-1 )
579 cc = t( i, i-1 )
580 bb = t( i-1, i )
581 dd = t( i, i )
582 CALL dlanv2( aa, bb, cc, dd, sr( kwtop+i-2 ),
583 $ si( kwtop+i-2 ), sr( kwtop+i-1 ),
584 $ si( kwtop+i-1 ), cs, sn )
585 i = i - 2
586 END IF
587 GO TO 60
588 END IF
589*
590 IF( ns.LT.jw .OR. s.EQ.zero ) THEN
591 IF( ns.GT.1 .AND. s.NE.zero ) THEN
592*
593* ==== Reflect spike back into lower triangle ====
594*
595 CALL dcopy( ns, v, ldv, work, 1 )
596 beta = work( 1 )
597 CALL dlarfg( ns, beta, work( 2 ), 1, tau )
598*
599 CALL dlaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ),
600 $ ldt )
601*
602 CALL dlarf1f( 'L', ns, jw, work, 1, tau, t, ldt,
603 $ work( jw+1 ) )
604 CALL dlarf1f( 'R', ns, ns, work, 1, tau, t, ldt,
605 $ work( jw+1 ) )
606 CALL dlarf1f( 'R', jw, ns, work, 1, tau, v, ldv,
607 $ work( jw+1 ) )
608*
609 CALL dgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
610 $ lwork-jw, info )
611 END IF
612*
613* ==== Copy updated reduced window into place ====
614*
615 IF( kwtop.GT.1 )
616 $ h( kwtop, kwtop-1 ) = s*v( 1, 1 )
617 CALL dlacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
618 CALL dcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
619 $ ldh+1 )
620*
621* ==== Accumulate orthogonal matrix in order update
622* . H and Z, if requested. ====
623*
624 IF( ns.GT.1 .AND. s.NE.zero )
625 $ CALL dormhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v,
626 $ ldv,
627 $ work( jw+1 ), lwork-jw, info )
628*
629* ==== Update vertical slab in H ====
630*
631 IF( wantt ) THEN
632 ltop = 1
633 ELSE
634 ltop = ktop
635 END IF
636 DO 70 krow = ltop, kwtop - 1, nv
637 kln = min( nv, kwtop-krow )
638 CALL dgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
639 $ ldh, v, ldv, zero, wv, ldwv )
640 CALL dlacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ),
641 $ ldh )
642 70 CONTINUE
643*
644* ==== Update horizontal slab in H ====
645*
646 IF( wantt ) THEN
647 DO 80 kcol = kbot + 1, n, nh
648 kln = min( nh, n-kcol+1 )
649 CALL dgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
650 $ h( kwtop, kcol ), ldh, zero, t, ldt )
651 CALL dlacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
652 $ ldh )
653 80 CONTINUE
654 END IF
655*
656* ==== Update vertical slab in Z ====
657*
658 IF( wantz ) THEN
659 DO 90 krow = iloz, ihiz, nv
660 kln = min( nv, ihiz-krow+1 )
661 CALL dgemm( 'N', 'N', kln, jw, jw, one, z( krow,
662 $ kwtop ),
663 $ ldz, v, ldv, zero, wv, ldwv )
664 CALL dlacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
665 $ ldz )
666 90 CONTINUE
667 END IF
668 END IF
669*
670* ==== Return the number of deflations ... ====
671*
672 nd = jw - ns
673*
674* ==== ... and the number of shifts. (Subtracting
675* . INFQR from the spike length takes care
676* . of the case of a rare QR failure while
677* . calculating eigenvalues of the deflation
678* . window.) ====
679*
680 ns = ns - infqr
681*
682* ==== Return optimal workspace. ====
683*
684 work( 1 ) = dble( lwkopt )
685*
686* ==== End of DLAQR2 ====
687*
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:166
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:101
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:205
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:125
dlarf1f
subroutine dlarf1f(side, m, n, v, incv, tau, c, ldc, work)
DLARF1F applies an elementary reflector to a general rectangular
Definition dlarf1f.f:157
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:104
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:108
dtrexc
subroutine dtrexc(compq, n, t, ldt, q, ldq, ifst, ilst, work, info)
DTREXC
Definition dtrexc.f:146
dormhr
subroutine dormhr(side, trans, m, n, ilo, ihi, a, lda, tau, c, ldc, work, lwork, info)
DORMHR
Definition dormhr.f:176
Here is the call graph for this function:
Here is the caller graph for this function: