dlaqr3.3lapack

Langue: en

Version: 299298 (debian - 07/07/09)

Section: 3 (Bibliothèques de fonctions)

NAME

SYNOPSIS

SUBROUTINE DLAQR3(
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, LWK3, LWKOPT, NMIN

    
LOGICAL BULGE, SORTED

    
DOUBLE PRECISION DLAMCH

    
INTEGER ILAENV

    
EXTERNAL DLAMCH, ILAENV

    
EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR, DLANV2, DLAQR4, 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( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV, WORK, -1, INFO )

    
LWK2 = INT( WORK( 1 ) )

    
CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR, SI, 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 ) = DBLE( LWKOPT )

    
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 = ONE / 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

    
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( 'U', 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( 'A', JW, JW, ZERO, ONE, V, LDV )

    
NMIN = ILAENV( 12, 'DLAQR3', 'SV', JW, 1, JW, LWORK )

    
IF( JW.GT.NMIN ) THEN

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

    
ELSE

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

    
END IF

    
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( 'V', 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( 'V', 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( 'V', 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( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )

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

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

    
CALL DLARF( 'R', 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( 'U', 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( '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 70 KROW = LTOP, KWTOP - 1, NV

    
KLN = MIN( NV, KWTOP-KROW )

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

    
CALL DLACPY( 'A', 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( 'C', 'N', JW, KLN, JW, ONE, V, LDV, H( KWTOP, KCOL ), LDH, ZERO, T, LDT )

    
CALL DLACPY( 'A', 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( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ), LDZ, V, LDV, ZERO, WV, LDWV )

    
CALL DLACPY( 'A', 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