SYNOPSIS¶
- SUBROUTINE
ZLAQR3(
- 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, LWK3, LWKOPT, NMIN DOUBLE PRECISION DLAMCH INTEGER ILAENV
EXTERNAL DLAMCH, ILAENV EXTERNAL DLABAD, ZCOPY, ZGEHRD, ZGEMM, ZLACPY,
ZLAHQR, ZLAQR4, 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 ) ) CALL ZLAQR4(
.true., .true., JW, 1, JW, T, LDT, SH, 1, JW, V, LDV, WORK, -1, INFQR ) LWK3
= INT( WORK( 1 ) ) LWKOPT = MAX( JW+MAX( LWK1, LWK2 ), LWK3 ) 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 ) NMIN
= ILAENV( 12, 'ZLAQR3', 'SV', JW, 1, JW, LWORK ) IF( JW.GT.NMIN ) THEN CALL
ZLAQR4( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP ), 1, JW, V, LDV, WORK,
LWORK, INFQR ) ELSE CALL ZLAHQR( .true., .true., JW, 1, JW, T, LDT, SH(
KWTOP ), 1, JW, V, LDV, INFQR ) END IF 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