1 SUBROUTINE DLAED5( I, D, Z, DELTA, RHO, DLAM ) 2* 3* -- LAPACK routine (instrumented to count operations, version 3.0) -- 4* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, 5* Courant Institute, NAG Ltd., and Rice University 6* September 30, 1994 7* 8* .. Scalar Arguments .. 9 INTEGER I 10 DOUBLE PRECISION DLAM, RHO 11* .. 12* .. Array Arguments .. 13 DOUBLE PRECISION D( 2 ), DELTA( 2 ), Z( 2 ) 14* .. 15* Common block to return operation count and iteration count 16* ITCNT is unchanged, OPS is only incremented 17* .. Common blocks .. 18 COMMON / LATIME / OPS, ITCNT 19* .. 20* .. Scalars in Common .. 21 DOUBLE PRECISION ITCNT, OPS 22* .. 23* 24* Purpose 25* ======= 26* 27* This subroutine computes the I-th eigenvalue of a symmetric rank-one 28* modification of a 2-by-2 diagonal matrix 29* 30* diag( D ) + RHO * Z * transpose(Z) . 31* 32* The diagonal elements in the array D are assumed to satisfy 33* 34* D(i) < D(j) for i < j . 35* 36* We also assume RHO > 0 and that the Euclidean norm of the vector 37* Z is one. 38* 39* Arguments 40* ========= 41* 42* I (input) INTEGER 43* The index of the eigenvalue to be computed. I = 1 or I = 2. 44* 45* D (input) DOUBLE PRECISION array, dimension (2) 46* The original eigenvalues. We assume D(1) < D(2). 47* 48* Z (input) DOUBLE PRECISION array, dimension (2) 49* The components of the updating vector. 50* 51* DELTA (output) DOUBLE PRECISION array, dimension (2) 52* The vector DELTA contains the information necessary 53* to construct the eigenvectors. 54* 55* RHO (input) DOUBLE PRECISION 56* The scalar in the symmetric updating formula. 57* 58* DLAM (output) DOUBLE PRECISION 59* The computed lambda_I, the I-th updated eigenvalue. 60* 61* Further Details 62* =============== 63* 64* Based on contributions by 65* Ren-Cang Li, Computer Science Division, University of California 66* at Berkeley, USA 67* 68* ===================================================================== 69* 70* .. Parameters .. 71 DOUBLE PRECISION ZERO, ONE, TWO, FOUR 72 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 73 $ FOUR = 4.0D0 ) 74* .. 75* .. Local Scalars .. 76 DOUBLE PRECISION B, C, DEL, TAU, TEMP, W 77* .. 78* .. Intrinsic Functions .. 79 INTRINSIC ABS, SQRT 80* .. 81* .. Executable Statements .. 82* 83 DEL = D( 2 ) - D( 1 ) 84 IF( I.EQ.1 ) THEN 85 W = ONE + TWO*RHO*( Z( 2 )*Z( 2 )-Z( 1 )*Z( 1 ) ) / DEL 86 IF( W.GT.ZERO ) THEN 87 OPS = OPS + 33 88 B = DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) 89 C = RHO*Z( 1 )*Z( 1 )*DEL 90* 91* B > ZERO, always 92* 93 TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) ) 94 DLAM = D( 1 ) + TAU 95 DELTA( 1 ) = -Z( 1 ) / TAU 96 DELTA( 2 ) = Z( 2 ) / ( DEL-TAU ) 97 ELSE 98 OPS = OPS + 31 99 B = -DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) 100 C = RHO*Z( 2 )*Z( 2 )*DEL 101 IF( B.GT.ZERO ) THEN 102 TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) ) 103 ELSE 104 TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO 105 END IF 106 DLAM = D( 2 ) + TAU 107 DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) 108 DELTA( 2 ) = -Z( 2 ) / TAU 109 END IF 110 TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) 111 DELTA( 1 ) = DELTA( 1 ) / TEMP 112 DELTA( 2 ) = DELTA( 2 ) / TEMP 113 ELSE 114* 115* Now I=2 116* 117 OPS = OPS + 24 118 B = -DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) 119 C = RHO*Z( 2 )*Z( 2 )*DEL 120 IF( B.GT.ZERO ) THEN 121 TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO 122 ELSE 123 TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) ) 124 END IF 125 DLAM = D( 2 ) + TAU 126 DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) 127 DELTA( 2 ) = -Z( 2 ) / TAU 128 TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) 129 DELTA( 1 ) = DELTA( 1 ) / TEMP 130 DELTA( 2 ) = DELTA( 2 ) / TEMP 131 END IF 132 RETURN 133* 134* End OF DLAED5 135* 136 END 137