dlaqr2 (l) - Linux Manuals

NAME

SYNOPSIS

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 )

    

INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, LDZ, LWORK, N, ND, NH, NS, NV, NW

    

LOGICAL WANTT, WANTZ

    

DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ), V( LDV, * ), WORK( * ), WV( LDWV, * ), Z( LDZ, * )

    

DOUBLE PRECISION ZERO, ONE

    

PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )

    

DOUBLE PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S, SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP

    

INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL, KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, LWKOPT

    

LOGICAL BULGE, SORTED

    

DOUBLE PRECISION DLAMCH

    

EXTERNAL DLAMCH

    

EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR, DLANV2, DLARF, DLARFG, DLASET, DORMHR, DTREXC

    

INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT

    

JW = MIN( NW, KBOT-KTOP+1 )

    

IF( JW.LE.2 ) THEN

    

LWKOPT = 1

    

ELSE

    

CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )

    

LWK1 = INT( WORK( 1 ) )

    

CALL DORMHR( aqRaq, aqNaq, 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 ) = DBLE( LWKOPT )

    

RETURN

    

END IF

    

NS = 0

    

ND = 0

    

WORK( 1 ) = ONE

    

IF( KTOP.GT.KBOT ) RETURN

    

IF( NW.LT.1 ) RETURN

    

SAFMIN = DLAMCH( aqSAFE MINIMUMaq )

    

SAFMAX = ONE / SAFMIN

    

CALL DLABAD( SAFMIN, SAFMAX )

    

ULP = DLAMCH( aqPRECISIONaq )

    

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

    

SR( KWTOP ) = H( KWTOP, KWTOP )

    

SI( KWTOP ) = ZERO

    

NS = 1

    

ND = 0

    

IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( 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 DLACPY( aqUaq, JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )

    

CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )

    

CALL DLASET( aqAaq, JW, JW, ZERO, ONE, V, LDV )

    

CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), SI( KWTOP ), 1, JW, V, LDV, INFQR )

    

DO 10 J = 1, JW - 3

    

T( J+2, J ) = ZERO

    

T( J+3, J ) = ZERO

    

10 CONTINUE

    

IF( JW.GT.2 ) T( JW, JW-2 ) = ZERO

    

NS = JW

    

ILST = INFQR + 1

    

20 CONTINUE

    

IF( ILST.LE.NS ) THEN

    

IF( NS.EQ.1 ) THEN

    

BULGE = .FALSE.

    

ELSE

    

BULGE = T( NS, NS-1 ).NE.ZERO

    

END IF

    

IF(

    

FOO = ABS( T( NS, NS ) )

    

IF( FOO.EQ.ZERO ) FOO = ABS( S )

    

IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN

    

NS = NS - 1

    

ELSE

    

IFST = NS

    

CALL DTREXC( aqVaq, JW, T, LDT, V, LDV, IFST, ILST, WORK, INFO )

    

ILST = ILST + 1

    

END IF

    

ELSE

    

FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )* SQRT( ABS( T( NS-1, NS ) ) )

    

IF( FOO.EQ.ZERO ) FOO = ABS( S )

    

IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE. MAX( SMLNUM, ULP*FOO ) ) THEN

    

NS = NS - 2

    

ELSE

    

IFST = NS

    

CALL DTREXC( aqVaq, JW, T, LDT, V, LDV, IFST, ILST, WORK, INFO )

    

ILST = ILST + 2

    

END IF

    

END IF

    

GO TO 20

    

END IF

    

IF( NS.EQ.0 ) S = ZERO

    

IF( NS.LT.JW ) THEN

    

SORTED = .false.

    

I = NS + 1

    

30 CONTINUE

    

IF( SORTED ) GO TO 50

    

SORTED = .true.

    

KEND = I - 1

    

I = INFQR + 1

    

IF( I.EQ.NS ) THEN

    

K = I + 1

    

ELSE IF( T( I+1, I ).EQ.ZERO ) THEN

    

K = I + 1

    

ELSE

    

K = I + 2

    

END IF

    

40 CONTINUE

    

IF( K.LE.KEND ) THEN

    

IF( K.EQ.I+1 ) THEN

    

EVI = ABS( T( I, I ) )

    

ELSE

    

EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )* SQRT( ABS( T( I, I+1 ) ) )

    

END IF

    

IF( K.EQ.KEND ) THEN

    

EVK = ABS( T( K, K ) )

    

ELSE IF( T( K+1, K ).EQ.ZERO ) THEN

    

EVK = ABS( T( K, K ) )

    

ELSE

    

EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )* SQRT( ABS( T( K, K+1 ) ) )

    

END IF

    

IF( EVI.GE.EVK ) THEN

    

I = K

    

ELSE

    

SORTED = .false.

    

IFST = I

    

ILST = K

    

CALL DTREXC( aqVaq, JW, T, LDT, V, LDV, IFST, ILST, WORK, INFO )

    

IF( INFO.EQ.0 ) THEN

    

I = ILST

    

ELSE

    

I = K

    

END IF

    

END IF

    

IF( I.EQ.KEND ) THEN

    

K = I + 1

    

ELSE IF( T( I+1, I ).EQ.ZERO ) THEN

    

K = I + 1

    

ELSE

    

K = I + 2

    

END IF

    

GO TO 40

    

END IF

    

GO TO 30

    

50 CONTINUE

    

END IF

    

I = JW

    

60 CONTINUE

    

IF( I.GE.INFQR+1 ) THEN

    

IF( I.EQ.INFQR+1 ) THEN

    

SR( KWTOP+I-1 ) = T( I, I )

    

SI( KWTOP+I-1 ) = ZERO

    

I = I - 1

    

ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN

    

SR( KWTOP+I-1 ) = T( I, I )

    

SI( KWTOP+I-1 ) = ZERO

    

I = I - 1

    

ELSE

    

AA = T( I-1, I-1 )

    

CC = T( I, I-1 )

    

BB = T( I-1, I )

    

DD = T( I, I )

    

CALL DLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ), SI( KWTOP+I-2 ), SR( KWTOP+I-1 ), SI( KWTOP+I-1 ), CS, SN )

    

I = I - 2

    

END IF

    

GO TO 60

    

END IF

    

IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN

    

IF( NS.GT.1 .AND. S.NE.ZERO ) THEN

    

CALL DCOPY( NS, V, LDV, WORK, 1 )

    

BETA = WORK( 1 )

    

CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU )

    

WORK( 1 ) = ONE

    

CALL DLASET( aqLaq, JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )

    

CALL DLARF( aqLaq, NS, JW, WORK, 1, TAU, T, LDT, WORK( JW+1 ) )

    

CALL DLARF( aqRaq, NS, NS, WORK, 1, TAU, T, LDT, WORK( JW+1 ) )

    

CALL DLARF( aqRaq, JW, NS, WORK, 1, TAU, V, LDV, WORK( JW+1 ) )

    

CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), LWORK-JW, INFO )

    

END IF

    

IF( KWTOP.GT.1 ) H( KWTOP, KWTOP-1 ) = S*V( 1, 1 )

    

CALL DLACPY( aqUaq, JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )

    

CALL DCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ), LDH+1 )

    

IF( NS.GT.1 .AND. S.NE.ZERO ) CALL DORMHR( aqRaq, aqNaq, 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 70 KROW = LTOP, KWTOP - 1, NV

    

KLN = MIN( NV, KWTOP-KROW )

    

CALL DGEMM( aqNaq, aqNaq, KLN, JW, JW, ONE, H( KROW, KWTOP ), LDH, V, LDV, ZERO, WV, LDWV )

    

CALL DLACPY( aqAaq, KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )

    

70 CONTINUE

    

IF( WANTT ) THEN

    

DO 80 KCOL = KBOT + 1, N, NH

    

KLN = MIN( NH, N-KCOL+1 )

    

CALL DGEMM( aqCaq, aqNaq, JW, KLN, JW, ONE, V, LDV, H( KWTOP, KCOL ), LDH, ZERO, T, LDT )

    

CALL DLACPY( aqAaq, JW, KLN, T, LDT, H( KWTOP, KCOL ), LDH )

    

80 CONTINUE

    

END IF

    

IF( WANTZ ) THEN

    

DO 90 KROW = ILOZ, IHIZ, NV

    

KLN = MIN( NV, IHIZ-KROW+1 )

    

CALL DGEMM( aqNaq, aqNaq, KLN, JW, JW, ONE, Z( KROW, KWTOP ), LDZ, V, LDV, ZERO, WV, LDWV )

    

CALL DLACPY( aqAaq, KLN, JW, WV, LDWV, Z( KROW, KWTOP ), LDZ )

    

90 CONTINUE

    

END IF

    

END IF

    

ND = JW - NS

    

NS = NS - INFQR

    

WORK( 1 ) = DBLE( LWKOPT )

    

END

PURPOSE