dlaqr5.3lapack

Langue: en

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

Section: 3 (Bibliothèques de fonctions)

NAME

SYNOPSIS

SUBROUTINE DLAQR5(
WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH )

    
INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV, LDWH, LDWV, LDZ, N, NH, NSHFTS, NV

    
LOGICAL WANTT, WANTZ

    
DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), U( LDU, * ), V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )

    
DOUBLE PRECISION ZERO, ONE

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

    
DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM, SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, TST1, TST2, ULP

    
INTEGER I, I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN, JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS, M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL, NS, NU

    
LOGICAL ACCUM, BLK22, BMP22

    
DOUBLE PRECISION DLAMCH

    
EXTERNAL DLAMCH

    
INTRINSIC ABS, DBLE, MAX, MIN, MOD

    
DOUBLE PRECISION VT( 3 )

    
EXTERNAL DGEMM, DLABAD, DLACPY, DLAQR1, DLARFG, DLASET, DTRMM

    
IF( NSHFTS.LT.2 ) RETURN

    
IF( KTOP.GE.KBOT ) RETURN

    
DO 10 I = 1, NSHFTS - 2, 2

    
IF( SI( I ).NE.-SI( I+1 ) ) THEN

    
SWAP = SR( I )

    
SR( I ) = SR( I+1 )

    
SR( I+1 ) = SR( I+2 )

    
SR( I+2 ) = SWAP

    
SWAP = SI( I )

    
SI( I ) = SI( I+1 )

    
SI( I+1 ) = SI( I+2 )

    
SI( I+2 ) = SWAP

    
END IF

    
10 CONTINUE

    
NS = NSHFTS - MOD( NSHFTS, 2 )

    
SAFMIN = DLAMCH( 'SAFE MINIMUM' )

    
SAFMAX = ONE / SAFMIN

    
CALL DLABAD( SAFMIN, SAFMAX )

    
ULP = DLAMCH( 'PRECISION' )

    
SMLNUM = SAFMIN*( DBLE( N ) / ULP )

    
ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )

    
BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )

    
IF( KTOP+2.LE.KBOT ) H( KTOP+2, KTOP ) = ZERO

    
NBMPS = NS / 2

    
KDU = 6*NBMPS - 3

    
DO 220 INCOL = 3*( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3*NBMPS - 2

    
NDCOL = INCOL + KDU

    
IF( ACCUM ) CALL DLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU )

    
DO 150 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )

    
MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )

    
MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )

    
M22 = MBOT + 1

    
BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ. ( KBOT-2 )

    
DO 20 M = MTOP, MBOT

    
K = KRCOL + 3*( M-1 )

    
IF( K.EQ.KTOP-1 ) THEN

    
CALL DLAQR1( 3, H( KTOP, KTOP ), LDH, SR( 2*M-1 ), SI( 2*M-1 ), SR( 2*M ), SI( 2*M ), V( 1, M ) )

    
ALPHA = V( 1, M )

    
CALL DLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )

    
ELSE

    
BETA = H( K+1, K )

    
V( 2, M ) = H( K+2, K )

    
V( 3, M ) = H( K+3, K )

    
CALL DLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) )

    
IF( H( K+3, K ).NE.ZERO .OR. H( K+3, K+1 ).NE. ZERO .OR. H( K+3, K+2 ).EQ.ZERO ) THEN

    
H( K+1, K ) = BETA

    
H( K+2, K ) = ZERO

    
H( K+3, K ) = ZERO

    
ELSE

    
CALL DLAQR1( 3, H( K+1, K+1 ), LDH, SR( 2*M-1 ), SI( 2*M-1 ), SR( 2*M ), SI( 2*M ), VT )

    
ALPHA = VT( 1 )

    
CALL DLARFG( 3, ALPHA, VT( 2 ), 1, VT( 1 ) )

    
REFSUM = VT( 1 )*( H( K+1, K )+VT( 2 )* H( K+2, K ) )

    
IF( ABS( H( K+2, K )-REFSUM*VT( 2 ) )+ ABS( REFSUM*VT( 3 ) ).GT.ULP* ( ABS( H( K, K ) )+ABS( H( K+1, K+1 ) )+ABS( H( K+2, K+2 ) ) ) ) THEN

    
H( K+1, K ) = BETA

    
H( K+2, K ) = ZERO

    
H( K+3, K ) = ZERO

    
ELSE

    
H( K+1, K ) = H( K+1, K ) - REFSUM

    
H( K+2, K ) = ZERO

    
H( K+3, K ) = ZERO

    
V( 1, M ) = VT( 1 )

    
V( 2, M ) = VT( 2 )

    
V( 3, M ) = VT( 3 )

    
END IF

    
END IF

    
END IF

    
20 CONTINUE

    
K = KRCOL + 3*( M22-1 )

    
IF( BMP22 ) THEN

    
IF( K.EQ.KTOP-1 ) THEN

    
CALL DLAQR1( 2, H( K+1, K+1 ), LDH, SR( 2*M22-1 ), SI( 2*M22-1 ), SR( 2*M22 ), SI( 2*M22 ), V( 1, M22 ) )

    
BETA = V( 1, M22 )

    
CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )

    
ELSE

    
BETA = H( K+1, K )

    
V( 2, M22 ) = H( K+2, K )

    
CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )

    
H( K+1, K ) = BETA

    
H( K+2, K ) = ZERO

    
END IF

    
END IF

    
IF( ACCUM ) THEN

    
JBOT = MIN( NDCOL, KBOT )

    
ELSE IF( WANTT ) THEN

    
JBOT = N

    
ELSE

    
JBOT = KBOT

    
END IF

    
DO 40 J = MAX( KTOP, KRCOL ), JBOT

    
MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )

    
DO 30 M = MTOP, MEND

    
K = KRCOL + 3*( M-1 )

    
REFSUM = V( 1, M )*( H( K+1, J )+V( 2, M )* H( K+2, J )+V( 3, M )*H( K+3, J ) )

    
H( K+1, J ) = H( K+1, J ) - REFSUM

    
H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M )

    
H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M )

    
30 CONTINUE

    
40 CONTINUE

    
IF( BMP22 ) THEN

    
K = KRCOL + 3*( M22-1 )

    
DO 50 J = MAX( K+1, KTOP ), JBOT

    
REFSUM = V( 1, M22 )*( H( K+1, J )+V( 2, M22 )* H( K+2, J ) )

    
H( K+1, J ) = H( K+1, J ) - REFSUM

    
H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 )

    
50 CONTINUE

    
END IF

    
IF( ACCUM ) THEN

    
JTOP = MAX( KTOP, INCOL )

    
ELSE IF( WANTT ) THEN

    
JTOP = 1

    
ELSE

    
JTOP = KTOP

    
END IF

    
DO 90 M = MTOP, MBOT

    
IF( V( 1, M ).NE.ZERO ) THEN

    
K = KRCOL + 3*( M-1 )

    
DO 60 J = JTOP, MIN( KBOT, K+3 )

    
REFSUM = V( 1, M )*( H( J, K+1 )+V( 2, M )* H( J, K+2 )+V( 3, M )*H( J, K+3 ) )

    
H( J, K+1 ) = H( J, K+1 ) - REFSUM

    
H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M )

    
H( J, K+3 ) = H( J, K+3 ) - REFSUM*V( 3, M )

    
60 CONTINUE

    
IF( ACCUM ) THEN

    
KMS = K - INCOL

    
DO 70 J = MAX( 1, KTOP-INCOL ), KDU

    
REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M )* U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) )

    
U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM

    
U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M )

    
U( J, KMS+3 ) = U( J, KMS+3 ) - REFSUM*V( 3, M )

    
70 CONTINUE

    
ELSE IF( WANTZ ) THEN

    
DO 80 J = ILOZ, IHIZ

    
REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )* Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) )

    
Z( J, K+1 ) = Z( J, K+1 ) - REFSUM

    
Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M )

    
Z( J, K+3 ) = Z( J, K+3 ) - REFSUM*V( 3, M )

    
80 CONTINUE

    
END IF

    
END IF

    
90 CONTINUE

    
K = KRCOL + 3*( M22-1 )

    
IF( BMP22 .AND. ( V( 1, M22 ).NE.ZERO ) ) THEN

    
DO 100 J = JTOP, MIN( KBOT, K+3 )

    
REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )* H( J, K+2 ) )

    
H( J, K+1 ) = H( J, K+1 ) - REFSUM

    
H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M22 )

    
100 CONTINUE

    
IF( ACCUM ) THEN

    
KMS = K - INCOL

    
DO 110 J = MAX( 1, KTOP-INCOL ), KDU

    
REFSUM = V( 1, M22 )*( U( J, KMS+1 )+V( 2, M22 )* U( J, KMS+2 ) )

    
U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM

    
U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M22 )

    
110 CONTINUE

    
ELSE IF( WANTZ ) THEN

    
DO 120 J = ILOZ, IHIZ

    
REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )* Z( J, K+2 ) )

    
Z( J, K+1 ) = Z( J, K+1 ) - REFSUM

    
Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M22 )

    
120 CONTINUE

    
END IF

    
END IF

    
MSTART = MTOP

    
IF( KRCOL+3*( MSTART-1 ).LT.KTOP ) MSTART = MSTART + 1

    
MEND = MBOT

    
IF( BMP22 ) MEND = MEND + 1

    
IF( KRCOL.EQ.KBOT-2 ) MEND = MEND + 1

    
DO 130 M = MSTART, MEND

    
K = MIN( KBOT-1, KRCOL+3*( M-1 ) )

    
IF( H( K+1, K ).NE.ZERO ) THEN

    
TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) )

    
IF( TST1.EQ.ZERO ) THEN

    
IF( K.GE.KTOP+1 ) TST1 = TST1 + ABS( H( K, K-1 ) )

    
IF( K.GE.KTOP+2 ) TST1 = TST1 + ABS( H( K, K-2 ) )

    
IF( K.GE.KTOP+3 ) TST1 = TST1 + ABS( H( K, K-3 ) )

    
IF( K.LE.KBOT-2 ) TST1 = TST1 + ABS( H( K+2, K+1 ) )

    
IF( K.LE.KBOT-3 ) TST1 = TST1 + ABS( H( K+3, K+1 ) )

    
IF( K.LE.KBOT-4 ) TST1 = TST1 + ABS( H( K+4, K+1 ) )

    
END IF

    
IF( ABS( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) ) THEN

    
H12 = MAX( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) )

    
H21 = MIN( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) )

    
H11 = MAX( ABS( H( K+1, K+1 ) ), ABS( H( K, K )-H( K+1, K+1 ) ) )

    
H22 = MIN( ABS( H( K+1, K+1 ) ), ABS( H( K, K )-H( K+1, K+1 ) ) )

    
SCL = H11 + H12

    
TST2 = H22*( H11 / SCL )

    
IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO

    
END IF

    
END IF

    
130 CONTINUE

    
MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )

    
DO 140 M = MTOP, MEND

    
K = KRCOL + 3*( M-1 )

    
REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )

    
H( K+4, K+1 ) = -REFSUM

    
H( K+4, K+2 ) = -REFSUM*V( 2, M )

    
H( K+4, K+3 ) = H( K+4, K+3 ) - REFSUM*V( 3, M )

    
140 CONTINUE

    
150 CONTINUE

    
IF( ACCUM ) THEN

    
IF( WANTT ) THEN

    
JTOP = 1

    
JBOT = N

    
ELSE

    
JTOP = KTOP

    
JBOT = KBOT

    
END IF

    
IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR. ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN

    
K1 = MAX( 1, KTOP-INCOL )

    
NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1

    
DO 160 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH

    
JLEN = MIN( NH, JBOT-JCOL+1 )

    
CALL DGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ), LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH, LDWH )

    
CALL DLACPY( 'ALL', NU, JLEN, WH, LDWH, H( INCOL+K1, JCOL ), LDH )

    
160 CONTINUE

    
DO 170 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV

    
JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )

    
CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE, H( JROW, INCOL+K1 ), LDH, U( K1, K1 ), LDU, ZERO, WV, LDWV )

    
CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV, H( JROW, INCOL+K1 ), LDH )

    
170 CONTINUE

    
IF( WANTZ ) THEN

    
DO 180 JROW = ILOZ, IHIZ, NV

    
JLEN = MIN( NV, IHIZ-JROW+1 )

    
CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE, Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ), LDU, ZERO, WV, LDWV )

    
CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV, Z( JROW, INCOL+K1 ), LDZ )

    
180 CONTINUE

    
END IF

    
ELSE

    
I2 = ( KDU+1 ) / 2

    
I4 = KDU

    
J2 = I4 - I2

    
J4 = KDU

    
KZS = ( J4-J2 ) - ( NS+1 )

    
KNZ = NS + 1

    
DO 190 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH

    
JLEN = MIN( NH, JBOT-JCOL+1 )

    
CALL DLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ), LDH, WH( KZS+1, 1 ), LDWH )

    
CALL DLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH )

    
CALL DTRMM( 'L', 'U', 'C', 'N', KNZ, JLEN, ONE, U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ), LDWH )

    
CALL DGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU, H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH )

    
CALL DLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH, WH( I2+1, 1 ), LDWH )

    
CALL DTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE, U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH )

    
CALL DGEMM( 'C', 'N', I4-I2, JLEN, J4-J2, ONE, U( J2+1, I2+1 ), LDU, H( INCOL+1+J2, JCOL ), LDH, ONE, WH( I2+1, 1 ), LDWH )

    
CALL DLACPY( 'ALL', KDU, JLEN, WH, LDWH, H( INCOL+1, JCOL ), LDH )

    
190 CONTINUE

    
DO 200 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV

    
JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW )

    
CALL DLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2 ), LDH, WV( 1, 1+KZS ), LDWV )

    
CALL DLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV )

    
CALL DTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE, U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ), LDWV )

    
CALL DGEMM( 'N', 'N', JLEN, I2, J2, ONE, H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV, LDWV )

    
CALL DLACPY( 'ALL', JLEN, J2, H( JROW, INCOL+1 ), LDH, WV( 1, 1+I2 ), LDWV )

    
CALL DTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE, U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )

    
CALL DGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE, H( JROW, INCOL+1+J2 ), LDH, U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ), LDWV )

    
CALL DLACPY( 'ALL', JLEN, KDU, WV, LDWV, H( JROW, INCOL+1 ), LDH )

    
200 CONTINUE

    
IF( WANTZ ) THEN

    
DO 210 JROW = ILOZ, IHIZ, NV

    
JLEN = MIN( NV, IHIZ-JROW+1 )

    
CALL DLACPY( 'ALL', JLEN, KNZ, Z( JROW, INCOL+1+J2 ), LDZ, WV( 1, 1+KZS ), LDWV )

    
CALL DLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV )

    
CALL DTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE, U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ), LDWV )

    
CALL DGEMM( 'N', 'N', JLEN, I2, J2, ONE, Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE, WV, LDWV )

    
CALL DLACPY( 'ALL', JLEN, J2, Z( JROW, INCOL+1 ), LDZ, WV( 1, 1+I2 ), LDWV )

    
CALL DTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE, U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )

    
CALL DGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE, Z( JROW, INCOL+1+J2 ), LDZ, U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ), LDWV )

    
CALL DLACPY( 'ALL', JLEN, KDU, WV, LDWV, Z( JROW, INCOL+1 ), LDZ )

    
210 CONTINUE

    
END IF

    
END IF

    
END IF

    
220 CONTINUE

    
END

PURPOSE