|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER I 00010 DOUBLE PRECISION DSIGMA, RHO 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * This subroutine computes the square root of the I-th eigenvalue 00020 * of a positive symmetric rank-one modification of a 2-by-2 diagonal 00021 * matrix 00022 * 00023 * diag( D ) * diag( D ) + RHO * Z * transpose(Z) . 00024 * 00025 * The diagonal entries in the array D are assumed to satisfy 00026 * 00027 * 0 <= D(i) < D(j) for i < j . 00028 * 00029 * We also assume RHO > 0 and that the Euclidean norm of the vector 00030 * Z is one. 00031 * 00032 * Arguments 00033 * ========= 00034 * 00035 * I (input) INTEGER 00036 * The index of the eigenvalue to be computed. I = 1 or I = 2. 00037 * 00038 * D (input) DOUBLE PRECISION array, dimension ( 2 ) 00039 * The original eigenvalues. We assume 0 <= D(1) < D(2). 00040 * 00041 * Z (input) DOUBLE PRECISION array, dimension ( 2 ) 00042 * The components of the updating vector. 00043 * 00044 * DELTA (output) DOUBLE PRECISION array, dimension ( 2 ) 00045 * Contains (D(j) - sigma_I) in its j-th component. 00046 * The vector DELTA contains the information necessary 00047 * to construct the eigenvectors. 00048 * 00049 * RHO (input) DOUBLE PRECISION 00050 * The scalar in the symmetric updating formula. 00051 * 00052 * DSIGMA (output) DOUBLE PRECISION 00053 * The computed sigma_I, the I-th updated eigenvalue. 00054 * 00055 * WORK (workspace) DOUBLE PRECISION array, dimension ( 2 ) 00056 * WORK contains (D(j) + sigma_I) in its j-th component. 00057 * 00058 * Further Details 00059 * =============== 00060 * 00061 * Based on contributions by 00062 * Ren-Cang Li, Computer Science Division, University of California 00063 * at Berkeley, USA 00064 * 00065 * ===================================================================== 00066 * 00067 * .. Parameters .. 00068 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR 00069 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, 00070 $ THREE = 3.0D+0, FOUR = 4.0D+0 ) 00071 * .. 00072 * .. Local Scalars .. 00073 DOUBLE PRECISION B, C, DEL, DELSQ, TAU, W 00074 * .. 00075 * .. Intrinsic Functions .. 00076 INTRINSIC ABS, SQRT 00077 * .. 00078 * .. Executable Statements .. 00079 * 00080 DEL = D( 2 ) - D( 1 ) 00081 DELSQ = DEL*( D( 2 )+D( 1 ) ) 00082 IF( I.EQ.1 ) THEN 00083 W = ONE + FOUR*RHO*( Z( 2 )*Z( 2 ) / ( D( 1 )+THREE*D( 2 ) )- 00084 $ Z( 1 )*Z( 1 ) / ( THREE*D( 1 )+D( 2 ) ) ) / DEL 00085 IF( W.GT.ZERO ) THEN 00086 B = DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) 00087 C = RHO*Z( 1 )*Z( 1 )*DELSQ 00088 * 00089 * B > ZERO, always 00090 * 00091 * The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 ) 00092 * 00093 TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) ) 00094 * 00095 * The following TAU is DSIGMA - D( 1 ) 00096 * 00097 TAU = TAU / ( D( 1 )+SQRT( D( 1 )*D( 1 )+TAU ) ) 00098 DSIGMA = D( 1 ) + TAU 00099 DELTA( 1 ) = -TAU 00100 DELTA( 2 ) = DEL - TAU 00101 WORK( 1 ) = TWO*D( 1 ) + TAU 00102 WORK( 2 ) = ( D( 1 )+TAU ) + D( 2 ) 00103 * DELTA( 1 ) = -Z( 1 ) / TAU 00104 * DELTA( 2 ) = Z( 2 ) / ( DEL-TAU ) 00105 ELSE 00106 B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) 00107 C = RHO*Z( 2 )*Z( 2 )*DELSQ 00108 * 00109 * The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 ) 00110 * 00111 IF( B.GT.ZERO ) THEN 00112 TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) ) 00113 ELSE 00114 TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO 00115 END IF 00116 * 00117 * The following TAU is DSIGMA - D( 2 ) 00118 * 00119 TAU = TAU / ( D( 2 )+SQRT( ABS( D( 2 )*D( 2 )+TAU ) ) ) 00120 DSIGMA = D( 2 ) + TAU 00121 DELTA( 1 ) = -( DEL+TAU ) 00122 DELTA( 2 ) = -TAU 00123 WORK( 1 ) = D( 1 ) + TAU + D( 2 ) 00124 WORK( 2 ) = TWO*D( 2 ) + TAU 00125 * DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) 00126 * DELTA( 2 ) = -Z( 2 ) / TAU 00127 END IF 00128 * TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) 00129 * DELTA( 1 ) = DELTA( 1 ) / TEMP 00130 * DELTA( 2 ) = DELTA( 2 ) / TEMP 00131 ELSE 00132 * 00133 * Now I=2 00134 * 00135 B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) 00136 C = RHO*Z( 2 )*Z( 2 )*DELSQ 00137 * 00138 * The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 ) 00139 * 00140 IF( B.GT.ZERO ) THEN 00141 TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO 00142 ELSE 00143 TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) ) 00144 END IF 00145 * 00146 * The following TAU is DSIGMA - D( 2 ) 00147 * 00148 TAU = TAU / ( D( 2 )+SQRT( D( 2 )*D( 2 )+TAU ) ) 00149 DSIGMA = D( 2 ) + TAU 00150 DELTA( 1 ) = -( DEL+TAU ) 00151 DELTA( 2 ) = -TAU 00152 WORK( 1 ) = D( 1 ) + TAU + D( 2 ) 00153 WORK( 2 ) = TWO*D( 2 ) + TAU 00154 * DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) 00155 * DELTA( 2 ) = -Z( 2 ) / TAU 00156 * TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) 00157 * DELTA( 1 ) = DELTA( 1 ) / TEMP 00158 * DELTA( 2 ) = DELTA( 2 ) / TEMP 00159 END IF 00160 RETURN 00161 * 00162 * End of DLASD5 00163 * 00164 END
1.7.3 
