274 SUBROUTINE slaqr3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
275 $ ihiz, z, ldz, ns, nd, sr, si, v, ldv, nh, t,
276 $ ldt, nv, wv, ldwv, work, lwork )
284 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
285 $ ldz, lwork, n, nd, nh, ns, nv, nw
289 REAL H( ldh, * ), SI( * ), SR( * ), T( ldt, * ),
290 $ v( ldv, * ), work( * ), wv( ldwv, * ),
297 parameter ( zero = 0.0e0, one = 1.0e0 )
300 REAL 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, lwk3,
305 LOGICAL BULGE, SORTED
310 EXTERNAL slamch, ilaenv
318 INTRINSIC abs, int, max, min,
REAL, SQRT
324 jw = min( nw, kbot-ktop+1 )
331 CALL sgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
332 lwk1 = int( work( 1 ) )
336 CALL sormhr(
'R',
'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
338 lwk2 = int( work( 1 ) )
342 CALL slaqr4( .true., .true., jw, 1, jw, t, ldt, sr, si, 1, jw,
343 $ v, ldv, work, -1, infqr )
344 lwk3 = int( work( 1 ) )
348 lwkopt = max( jw+max( lwk1, lwk2 ), lwk3 )
353 IF( lwork.EQ.-1 )
THEN
354 work( 1 ) =
REAL( lwkopt )
371 safmin = slamch(
'SAFE MINIMUM' )
372 safmax = one / safmin
373 CALL slabad( safmin, safmax )
374 ulp = slamch(
'PRECISION' )
375 smlnum = safmin*(
REAL( N ) / ULP )
379 jw = min( nw, kbot-ktop+1 )
380 kwtop = kbot - jw + 1
381 IF( kwtop.EQ.ktop )
THEN
384 s = h( kwtop, kwtop-1 )
387 IF( kbot.EQ.kwtop )
THEN
391 sr( kwtop ) = h( kwtop, kwtop )
395 IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
400 $ h( kwtop, kwtop-1 ) = zero
412 CALL slacpy(
'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
413 CALL scopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
415 CALL slaset(
'A', jw, jw, zero, one, v, ldv )
416 nmin = ilaenv( 12,
'SLAQR3',
'SV', jw, 1, jw, lwork )
417 IF( jw.GT.nmin )
THEN
418 CALL slaqr4( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
419 $ si( kwtop ), 1, jw, v, ldv, work, lwork, infqr )
421 CALL slahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
422 $ si( kwtop ), 1, jw, v, ldv, infqr )
432 $ t( jw, jw-2 ) = zero
439 IF( ilst.LE.ns )
THEN
443 bulge = t( ns, ns-1 ).NE.zero
448 IF( .NOT. bulge )
THEN
452 foo = abs( t( ns, ns ) )
455 IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
THEN
466 CALL strexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, work,
474 foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
475 $ sqrt( abs( t( ns-1, ns ) ) )
478 IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
479 $ max( smlnum, ulp*foo ) )
THEN
491 CALL strexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, work,
524 ELSE IF( t( i+1, i ).EQ.zero )
THEN
532 evi = abs( t( i, i ) )
534 evi = abs( t( i, i ) ) + sqrt( abs( t( i+1, i ) ) )*
535 $ sqrt( abs( t( i, i+1 ) ) )
539 evk = abs( t( k, k ) )
540 ELSE IF( t( k+1, k ).EQ.zero )
THEN
541 evk = abs( t( k, k ) )
543 evk = abs( t( k, k ) ) + sqrt( abs( t( k+1, k ) ) )*
544 $ sqrt( abs( t( k, k+1 ) ) )
547 IF( evi.GE.evk )
THEN
553 CALL strexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, work,
563 ELSE IF( t( i+1, i ).EQ.zero )
THEN
578 IF( i.GE.infqr+1 )
THEN
579 IF( i.EQ.infqr+1 )
THEN
580 sr( kwtop+i-1 ) = t( i, i )
581 si( kwtop+i-1 ) = zero
583 ELSE IF( t( i, i-1 ).EQ.zero )
THEN
584 sr( kwtop+i-1 ) = t( i, i )
585 si( kwtop+i-1 ) = zero
592 CALL slanv2( aa, bb, cc, dd, sr( kwtop+i-2 ),
593 $ si( kwtop+i-2 ), sr( kwtop+i-1 ),
594 $ si( kwtop+i-1 ), cs, sn )
600 IF( ns.LT.jw .OR. s.EQ.zero )
THEN
601 IF( ns.GT.1 .AND. s.NE.zero )
THEN
605 CALL scopy( ns, v, ldv, work, 1 )
607 CALL slarfg( ns, beta, work( 2 ), 1, tau )
610 CALL slaset(
'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
612 CALL slarf(
'L', ns, jw, work, 1, tau, t, ldt,
614 CALL slarf(
'R', ns, ns, work, 1, tau, t, ldt,
616 CALL slarf(
'R', jw, ns, work, 1, tau, v, ldv,
619 CALL sgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
626 $ h( kwtop, kwtop-1 ) = s*v( 1, 1 )
627 CALL slacpy(
'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
628 CALL scopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
634 IF( ns.GT.1 .AND. s.NE.zero )
635 $
CALL sormhr(
'R',
'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
636 $ work( jw+1 ), lwork-jw, info )
645 DO 70 krow = ltop, kwtop - 1, nv
646 kln = min( nv, kwtop-krow )
647 CALL sgemm(
'N',
'N', kln, jw, jw, one, h( krow, kwtop ),
648 $ ldh, v, ldv, zero, wv, ldwv )
649 CALL slacpy(
'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
655 DO 80 kcol = kbot + 1, n, nh
656 kln = min( nh, n-kcol+1 )
657 CALL sgemm(
'C',
'N', jw, kln, jw, one, v, ldv,
658 $ h( kwtop, kcol ), ldh, zero, t, ldt )
659 CALL slacpy(
'A', jw, kln, t, ldt, h( kwtop, kcol ),
667 DO 90 krow = iloz, ihiz, nv
668 kln = min( nv, ihiz-krow+1 )
669 CALL sgemm(
'N',
'N', kln, jw, jw, one, z( krow, kwtop ),
670 $ ldz, v, ldv, zero, wv, ldwv )
671 CALL slacpy(
'A', kln, jw, wv, ldwv, z( krow, kwtop ),
691 work( 1 ) =
REAL( lwkopt )
subroutine slahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
subroutine slanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form...
subroutine slaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
SLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine strexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
STREXC
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
subroutine slaqr3(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)
SLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMHR
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.