1*> \brief \b DLAED9 used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLAED9 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed9.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed9.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed9.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, 22* S, LDS, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N 26* DOUBLE PRECISION RHO 27* .. 28* .. Array Arguments .. 29* DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), 30* $ W( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DLAED9 finds the roots of the secular equation, as defined by the 40*> values in D, Z, and RHO, between KSTART and KSTOP. It makes the 41*> appropriate calls to DLAED4 and then stores the new matrix of 42*> eigenvectors for use in calculating the next level of Z vectors. 43*> \endverbatim 44* 45* Arguments: 46* ========== 47* 48*> \param[in] K 49*> \verbatim 50*> K is INTEGER 51*> The number of terms in the rational function to be solved by 52*> DLAED4. K >= 0. 53*> \endverbatim 54*> 55*> \param[in] KSTART 56*> \verbatim 57*> KSTART is INTEGER 58*> \endverbatim 59*> 60*> \param[in] KSTOP 61*> \verbatim 62*> KSTOP is INTEGER 63*> The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP 64*> are to be computed. 1 <= KSTART <= KSTOP <= K. 65*> \endverbatim 66*> 67*> \param[in] N 68*> \verbatim 69*> N is INTEGER 70*> The number of rows and columns in the Q matrix. 71*> N >= K (delation may result in N > K). 72*> \endverbatim 73*> 74*> \param[out] D 75*> \verbatim 76*> D is DOUBLE PRECISION array, dimension (N) 77*> D(I) contains the updated eigenvalues 78*> for KSTART <= I <= KSTOP. 79*> \endverbatim 80*> 81*> \param[out] Q 82*> \verbatim 83*> Q is DOUBLE PRECISION array, dimension (LDQ,N) 84*> \endverbatim 85*> 86*> \param[in] LDQ 87*> \verbatim 88*> LDQ is INTEGER 89*> The leading dimension of the array Q. LDQ >= max( 1, N ). 90*> \endverbatim 91*> 92*> \param[in] RHO 93*> \verbatim 94*> RHO is DOUBLE PRECISION 95*> The value of the parameter in the rank one update equation. 96*> RHO >= 0 required. 97*> \endverbatim 98*> 99*> \param[in] DLAMDA 100*> \verbatim 101*> DLAMDA is DOUBLE PRECISION array, dimension (K) 102*> The first K elements of this array contain the old roots 103*> of the deflated updating problem. These are the poles 104*> of the secular equation. 105*> \endverbatim 106*> 107*> \param[in] W 108*> \verbatim 109*> W is DOUBLE PRECISION array, dimension (K) 110*> The first K elements of this array contain the components 111*> of the deflation-adjusted updating vector. 112*> \endverbatim 113*> 114*> \param[out] S 115*> \verbatim 116*> S is DOUBLE PRECISION array, dimension (LDS, K) 117*> Will contain the eigenvectors of the repaired matrix which 118*> will be stored for subsequent Z vector calculation and 119*> multiplied by the previously accumulated eigenvectors 120*> to update the system. 121*> \endverbatim 122*> 123*> \param[in] LDS 124*> \verbatim 125*> LDS is INTEGER 126*> The leading dimension of S. LDS >= max( 1, K ). 127*> \endverbatim 128*> 129*> \param[out] INFO 130*> \verbatim 131*> INFO is INTEGER 132*> = 0: successful exit. 133*> < 0: if INFO = -i, the i-th argument had an illegal value. 134*> > 0: if INFO = 1, an eigenvalue did not converge 135*> \endverbatim 136* 137* Authors: 138* ======== 139* 140*> \author Univ. of Tennessee 141*> \author Univ. of California Berkeley 142*> \author Univ. of Colorado Denver 143*> \author NAG Ltd. 144* 145*> \date December 2016 146* 147*> \ingroup auxOTHERcomputational 148* 149*> \par Contributors: 150* ================== 151*> 152*> Jeff Rutter, Computer Science Division, University of California 153*> at Berkeley, USA 154* 155* ===================================================================== 156 SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, 157 $ S, LDS, INFO ) 158* 159* -- LAPACK computational routine (version 3.7.0) -- 160* -- LAPACK is a software package provided by Univ. of Tennessee, -- 161* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 162* December 2016 163* 164* .. Scalar Arguments .. 165 INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N 166 DOUBLE PRECISION RHO 167* .. 168* .. Array Arguments .. 169 DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), 170 $ W( * ) 171* .. 172* 173* ===================================================================== 174* 175* .. Local Scalars .. 176 INTEGER I, J 177 DOUBLE PRECISION TEMP 178* .. 179* .. External Functions .. 180 DOUBLE PRECISION DLAMC3, DNRM2 181 EXTERNAL DLAMC3, DNRM2 182* .. 183* .. External Subroutines .. 184 EXTERNAL DCOPY, DLAED4, XERBLA 185* .. 186* .. Intrinsic Functions .. 187 INTRINSIC MAX, SIGN, SQRT 188* .. 189* .. Executable Statements .. 190* 191* Test the input parameters. 192* 193 INFO = 0 194* 195 IF( K.LT.0 ) THEN 196 INFO = -1 197 ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN 198 INFO = -2 199 ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) ) 200 $ THEN 201 INFO = -3 202 ELSE IF( N.LT.K ) THEN 203 INFO = -4 204 ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN 205 INFO = -7 206 ELSE IF( LDS.LT.MAX( 1, K ) ) THEN 207 INFO = -12 208 END IF 209 IF( INFO.NE.0 ) THEN 210 CALL XERBLA( 'DLAED9', -INFO ) 211 RETURN 212 END IF 213* 214* Quick return if possible 215* 216 IF( K.EQ.0 ) 217 $ RETURN 218* 219* Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can 220* be computed with high relative accuracy (barring over/underflow). 221* This is a problem on machines without a guard digit in 222* add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). 223* The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), 224* which on any of these machines zeros out the bottommost 225* bit of DLAMDA(I) if it is 1; this makes the subsequent 226* subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation 227* occurs. On binary machines with a guard digit (almost all 228* machines) it does not change DLAMDA(I) at all. On hexadecimal 229* and decimal machines with a guard digit, it slightly 230* changes the bottommost bits of DLAMDA(I). It does not account 231* for hexadecimal or decimal machines without guard digits 232* (we know of none). We use a subroutine call to compute 233* 2*DLAMBDA(I) to prevent optimizing compilers from eliminating 234* this code. 235* 236 DO 10 I = 1, N 237 DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) 238 10 CONTINUE 239* 240 DO 20 J = KSTART, KSTOP 241 CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) 242* 243* If the zero finder fails, the computation is terminated. 244* 245 IF( INFO.NE.0 ) 246 $ GO TO 120 247 20 CONTINUE 248* 249 IF( K.EQ.1 .OR. K.EQ.2 ) THEN 250 DO 40 I = 1, K 251 DO 30 J = 1, K 252 S( J, I ) = Q( J, I ) 253 30 CONTINUE 254 40 CONTINUE 255 GO TO 120 256 END IF 257* 258* Compute updated W. 259* 260 CALL DCOPY( K, W, 1, S, 1 ) 261* 262* Initialize W(I) = Q(I,I) 263* 264 CALL DCOPY( K, Q, LDQ+1, W, 1 ) 265 DO 70 J = 1, K 266 DO 50 I = 1, J - 1 267 W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 268 50 CONTINUE 269 DO 60 I = J + 1, K 270 W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 271 60 CONTINUE 272 70 CONTINUE 273 DO 80 I = 1, K 274 W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) ) 275 80 CONTINUE 276* 277* Compute eigenvectors of the modified rank-1 modification. 278* 279 DO 110 J = 1, K 280 DO 90 I = 1, K 281 Q( I, J ) = W( I ) / Q( I, J ) 282 90 CONTINUE 283 TEMP = DNRM2( K, Q( 1, J ), 1 ) 284 DO 100 I = 1, K 285 S( I, J ) = Q( I, J ) / TEMP 286 100 CONTINUE 287 110 CONTINUE 288* 289 120 CONTINUE 290 RETURN 291* 292* End of DLAED9 293* 294 END 295