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