277 SUBROUTINE dlaqr2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
278 $ ihiz, z, ldz, ns, nd, sr, si, v, ldv, nh, t,
279 $ ldt, nv, wv, ldwv, work, lwork )
287 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
288 $ ldz, lwork, n, nd, nh, ns, nv, nw
292 DOUBLE PRECISION H( ldh, * ), SI( * ), SR( * ), T( ldt, * ),
293 $ v( ldv, * ), work( * ), wv( ldwv, * ),
299 DOUBLE PRECISION ZERO, ONE
300 parameter ( zero = 0.0d0, one = 1.0d0 )
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,
308 LOGICAL BULGE, SORTED
311 DOUBLE PRECISION DLAMCH
319 INTRINSIC abs, dble, int, max, min, sqrt
325 jw = min( nw, kbot-ktop+1 )
332 CALL dgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
333 lwk1 = int( work( 1 ) )
337 CALL dormhr(
'R',
'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
339 lwk2 = int( work( 1 ) )
343 lwkopt = jw + max( lwk1, lwk2 )
348 IF( lwork.EQ.-1 )
THEN
349 work( 1 ) = dble( lwkopt )
366 safmin = dlamch(
'SAFE MINIMUM' )
367 safmax = one / safmin
368 CALL dlabad( safmin, safmax )
369 ulp = dlamch(
'PRECISION' )
370 smlnum = safmin*( dble( n ) / ulp )
374 jw = min( nw, kbot-ktop+1 )
375 kwtop = kbot - jw + 1
376 IF( kwtop.EQ.ktop )
THEN
379 s = h( kwtop, kwtop-1 )
382 IF( kbot.EQ.kwtop )
THEN
386 sr( kwtop ) = h( kwtop, kwtop )
390 IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
395 $ h( kwtop, kwtop-1 ) = zero
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 )
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 )
421 $ t( jw, jw-2 ) = zero
428 IF( ilst.LE.ns )
THEN
432 bulge = t( ns, ns-1 ).NE.zero
437 IF( .NOT.bulge )
THEN
441 foo = abs( t( ns, ns ) )
444 IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
THEN
455 CALL dtrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, work,
463 foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
464 $ sqrt( abs( t( ns-1, ns ) ) )
467 IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
468 $ max( smlnum, ulp*foo ) )
THEN
480 CALL dtrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, work,
513 ELSE IF( t( i+1, i ).EQ.zero )
THEN
521 evi = abs( t( i, i ) )
523 evi = abs( t( i, i ) ) + sqrt( abs( t( i+1, i ) ) )*
524 $ sqrt( abs( t( i, i+1 ) ) )
528 evk = abs( t( k, k ) )
529 ELSE IF( t( k+1, k ).EQ.zero )
THEN
530 evk = abs( t( k, k ) )
532 evk = abs( t( k, k ) ) + sqrt( abs( t( k+1, k ) ) )*
533 $ sqrt( abs( t( k, k+1 ) ) )
536 IF( evi.GE.evk )
THEN
542 CALL dtrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, work,
552 ELSE IF( t( i+1, i ).EQ.zero )
THEN
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
572 ELSE IF( t( i, i-1 ).EQ.zero )
THEN
573 sr( kwtop+i-1 ) = t( i, i )
574 si( kwtop+i-1 ) = zero
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 )
589 IF( ns.LT.jw .OR. s.EQ.zero )
THEN
590 IF( ns.GT.1 .AND. s.NE.zero )
THEN
594 CALL dcopy( ns, v, ldv, work, 1 )
596 CALL dlarfg( ns, beta, work( 2 ), 1, tau )
599 CALL dlaset(
'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
601 CALL dlarf(
'L', ns, jw, work, 1, tau, t, ldt,
603 CALL dlarf(
'R', ns, ns, work, 1, tau, t, ldt,
605 CALL dlarf(
'R', jw, ns, work, 1, tau, v, ldv,
608 CALL dgehrd( jw, 1, ns, t, ldt, work, work( jw+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 ),
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 )
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 )
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 ),
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 ),
680 work( 1 ) = dble( lwkopt )
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...
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dtrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
DTREXC
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
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.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
DLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
subroutine dormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMHR
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...