SUBROUTINE DLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER I DOUBLE PRECISION DSIGMA, RHO * .. * .. Array Arguments .. DOUBLE PRECISION D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 ) * .. * * Purpose * ======= * * This subroutine computes the square root of the I-th eigenvalue * of a positive symmetric rank-one modification of a 2-by-2 diagonal * matrix * * diag( D ) * diag( D ) + RHO * Z * transpose(Z) . * * The diagonal entries in the array D are assumed to satisfy * * 0 <= D(i) < D(j) for i < j . * * We also assume RHO > 0 and that the Euclidean norm of the vector * Z is one. * * Arguments * ========= * * I (input) INTEGER * The index of the eigenvalue to be computed. I = 1 or I = 2. * * D (input) DOUBLE PRECISION array, dimension ( 2 ) * The original eigenvalues. We assume 0 <= D(1) < D(2). * * Z (input) DOUBLE PRECISION array, dimension ( 2 ) * The components of the updating vector. * * DELTA (output) DOUBLE PRECISION array, dimension ( 2 ) * Contains (D(j) - sigma_I) in its j-th component. * The vector DELTA contains the information necessary * to construct the eigenvectors. * * RHO (input) DOUBLE PRECISION * The scalar in the symmetric updating formula. * * DSIGMA (output) DOUBLE PRECISION * The computed sigma_I, the I-th updated eigenvalue. * * WORK (workspace) DOUBLE PRECISION array, dimension ( 2 ) * WORK contains (D(j) + sigma_I) in its j-th component. * * Further Details * =============== * * Based on contributions by * Ren-Cang Li, Computer Science Division, University of California * at Berkeley, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, \$ THREE = 3.0D+0, FOUR = 4.0D+0 ) * .. * .. Local Scalars .. DOUBLE PRECISION B, C, DEL, DELSQ, TAU, W * .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. * .. Executable Statements .. * DEL = D( 2 ) - D( 1 ) DELSQ = DEL*( D( 2 )+D( 1 ) ) IF( I.EQ.1 ) THEN W = ONE + FOUR*RHO*( Z( 2 )*Z( 2 ) / ( D( 1 )+THREE*D( 2 ) )- \$ Z( 1 )*Z( 1 ) / ( THREE*D( 1 )+D( 2 ) ) ) / DEL IF( W.GT.ZERO ) THEN B = DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) C = RHO*Z( 1 )*Z( 1 )*DELSQ * * B > ZERO, always * * The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 ) * TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) ) * * The following TAU is DSIGMA - D( 1 ) * TAU = TAU / ( D( 1 )+SQRT( D( 1 )*D( 1 )+TAU ) ) DSIGMA = D( 1 ) + TAU DELTA( 1 ) = -TAU DELTA( 2 ) = DEL - TAU WORK( 1 ) = TWO*D( 1 ) + TAU WORK( 2 ) = ( D( 1 )+TAU ) + D( 2 ) * DELTA( 1 ) = -Z( 1 ) / TAU * DELTA( 2 ) = Z( 2 ) / ( DEL-TAU ) ELSE B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) C = RHO*Z( 2 )*Z( 2 )*DELSQ * * The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 ) * IF( B.GT.ZERO ) THEN TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) ) ELSE TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO END IF * * The following TAU is DSIGMA - D( 2 ) * TAU = TAU / ( D( 2 )+SQRT( ABS( D( 2 )*D( 2 )+TAU ) ) ) DSIGMA = D( 2 ) + TAU DELTA( 1 ) = -( DEL+TAU ) DELTA( 2 ) = -TAU WORK( 1 ) = D( 1 ) + TAU + D( 2 ) WORK( 2 ) = TWO*D( 2 ) + TAU * DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) * DELTA( 2 ) = -Z( 2 ) / TAU END IF * TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) * DELTA( 1 ) = DELTA( 1 ) / TEMP * DELTA( 2 ) = DELTA( 2 ) / TEMP ELSE * * Now I=2 * B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) C = RHO*Z( 2 )*Z( 2 )*DELSQ * * The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 ) * IF( B.GT.ZERO ) THEN TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO ELSE TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) ) END IF * * The following TAU is DSIGMA - D( 2 ) * TAU = TAU / ( D( 2 )+SQRT( D( 2 )*D( 2 )+TAU ) ) DSIGMA = D( 2 ) + TAU DELTA( 1 ) = -( DEL+TAU ) DELTA( 2 ) = -TAU WORK( 1 ) = D( 1 ) + TAU + D( 2 ) WORK( 2 ) = TWO*D( 2 ) + TAU * DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) * DELTA( 2 ) = -Z( 2 ) / TAU * TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) * DELTA( 1 ) = DELTA( 1 ) / TEMP * DELTA( 2 ) = DELTA( 2 ) / TEMP END IF RETURN * * End of DLASD5 * END