SYNOPSIS¶
- SUBROUTINE
ZLAQR2(
- WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH,
V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK )
INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, LDZ, LWORK,
N, ND, NH, NS, NV, NW LOGICAL WANTT, WANTZ COMPLEX*16 H( LDH, * ), SH( * ),
T( LDT, * ), V( LDV, * ), WORK( * ), WV( LDWV, * ), Z( LDZ, * ) COMPLEX*16
ZERO, ONE PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ), ONE = ( 1.0d0, 0.0d0 ) )
DOUBLE PRECISION RZERO, RONE PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0 )
COMPLEX*16 BETA, CDUM, S, TAU DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM,
ULP INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN, KNT, KROW, KWTOP,
LTOP, LWK1, LWK2, LWKOPT DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH EXTERNAL
DLABAD, ZCOPY, ZGEHRD, ZGEMM, ZLACPY, ZLAHQR, ZLARF, ZLARFG, ZLASET, ZTREXC,
ZUNMHR INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, MAX, MIN DOUBLE
PRECISION CABS1 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
JW = MIN( NW, KBOT-KTOP+1 ) IF( JW.LE.2 ) THEN LWKOPT = 1 ELSE CALL ZGEHRD(
JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO ) LWK1 = INT( WORK( 1 ) ) CALL
ZUNMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV, WORK, -1, INFO )
LWK2 = INT( WORK( 1 ) ) LWKOPT = JW + MAX( LWK1, LWK2 ) END IF IF(
LWORK.EQ.-1 ) THEN WORK( 1 ) = DCMPLX( LWKOPT, 0 ) RETURN END IF NS = 0 ND =
0 WORK( 1 ) = ONE IF( KTOP.GT.KBOT ) RETURN IF( NW.LT.1 ) RETURN SAFMIN =
DLAMCH( 'SAFE MINIMUM' ) SAFMAX = RONE / SAFMIN CALL DLABAD( SAFMIN, SAFMAX
) ULP = DLAMCH( 'PRECISION' ) SMLNUM = SAFMIN*( DBLE( N ) / ULP ) JW = MIN(
NW, KBOT-KTOP+1 ) KWTOP = KBOT - JW + 1 IF( KWTOP.EQ.KTOP ) THEN S = ZERO
ELSE S = H( KWTOP, KWTOP-1 ) END IF IF( KBOT.EQ.KWTOP ) THEN SH( KWTOP ) =
H( KWTOP, KWTOP ) NS = 1 ND = 0 IF( CABS1( S ).LE.MAX( SMLNUM, ULP*CABS1( H(
KWTOP, KWTOP ) ) ) ) THEN NS = 0 ND = 1 IF( KWTOP.GT.KTOP ) H( KWTOP,
KWTOP-1 ) = ZERO END IF WORK( 1 ) = ONE RETURN END IF CALL ZLACPY( 'U', JW,
JW, H( KWTOP, KWTOP ), LDH, T, LDT ) CALL ZCOPY( JW-1, H( KWTOP+1, KWTOP ),
LDH+1, T( 2, 1 ), LDT+1 ) CALL ZLASET( 'A', JW, JW, ZERO, ONE, V, LDV ) CALL
ZLAHQR( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP ), 1, JW, V, LDV, INFQR
) NS = JW ILST = INFQR + 1 DO 10 KNT = INFQR + 1, JW FOO = CABS1( T( NS, NS
) ) IF( FOO.EQ.RZERO ) FOO = CABS1( S ) IF( CABS1( S )*CABS1( V( 1, NS )
).LE.MAX( SMLNUM, ULP*FOO ) ) THEN NS = NS - 1 ELSE IFST = NS CALL ZTREXC(
'V', JW, T, LDT, V, LDV, IFST, ILST, INFO ) ILST = ILST + 1 END IF 10
CONTINUE IF( NS.EQ.0 ) S = ZERO IF( NS.LT.JW ) THEN DO 30 I = INFQR + 1, NS
IFST = I DO 20 J = I + 1, NS IF( CABS1( T( J, J ) ).GT.CABS1( T( IFST, IFST
) ) ) IFST = J 20 CONTINUE ILST = I IF( IFST.NE.ILST ) CALL ZTREXC( 'V', JW,
T, LDT, V, LDV, IFST, ILST, INFO ) 30 CONTINUE END IF DO 40 I = INFQR + 1,
JW SH( KWTOP+I-1 ) = T( I, I ) 40 CONTINUE IF( NS.LT.JW .OR. S.EQ.ZERO )
THEN IF( NS.GT.1 .AND. S.NE.ZERO ) THEN CALL ZCOPY( NS, V, LDV, WORK, 1 ) DO
50 I = 1, NS WORK( I ) = DCONJG( WORK( I ) ) 50 CONTINUE BETA = WORK( 1 )
CALL ZLARFG( NS, BETA, WORK( 2 ), 1, TAU ) WORK( 1 ) = ONE CALL ZLASET( 'L',
JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT ) CALL ZLARF( 'L', NS, JW, WORK, 1,
DCONJG( TAU ), T, LDT, WORK( JW+1 ) ) CALL ZLARF( 'R', NS, NS, WORK, 1, TAU,
T, LDT, WORK( JW+1 ) ) CALL ZLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, WORK(
JW+1 ) ) CALL ZGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), LWORK-JW, INFO
) END IF IF( KWTOP.GT.1 ) H( KWTOP, KWTOP-1 ) = S*DCONJG( V( 1, 1 ) ) CALL
ZLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH ) CALL ZCOPY( JW-1, T(
2, 1 ), LDT+1, H( KWTOP+1, KWTOP ), LDH+1 ) IF( NS.GT.1 .AND. S.NE.ZERO )
CALL ZUNMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV, WORK( JW+1 ),
LWORK-JW, INFO ) IF( WANTT ) THEN LTOP = 1 ELSE LTOP = KTOP END IF DO 60
KROW = LTOP, KWTOP - 1, NV KLN = MIN( NV, KWTOP-KROW ) CALL ZGEMM( 'N', 'N',
KLN, JW, JW, ONE, H( KROW, KWTOP ), LDH, V, LDV, ZERO, WV, LDWV ) CALL
ZLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH ) 60 CONTINUE IF(
WANTT ) THEN DO 70 KCOL = KBOT + 1, N, NH KLN = MIN( NH, N-KCOL+1 ) CALL
ZGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, H( KWTOP, KCOL ), LDH, ZERO, T,
LDT ) CALL ZLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ), LDH ) 70 CONTINUE
END IF IF( WANTZ ) THEN DO 80 KROW = ILOZ, IHIZ, NV KLN = MIN( NV,
IHIZ-KROW+1 ) CALL ZGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ), LDZ,
V, LDV, ZERO, WV, LDWV ) CALL ZLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP
), LDZ ) 80 CONTINUE END IF END IF ND = JW - NS NS = NS - INFQR WORK( 1 ) =
DCMPLX( LWKOPT, 0 ) END