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