|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 |
SUBROUTINE DLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
* * -- LAPACK auxiliary routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver 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 |