1*> \brief \b DLAED8 used by DSTEDC. Merges eigenvalues and deflates secular equation. 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 DLAED8 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed8.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed8.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed8.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, 22* CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, 23* GIVCOL, GIVNUM, INDXP, INDX, INFO ) 24* 25* .. Scalar Arguments .. 26* INTEGER CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N, 27* $ QSIZ 28* DOUBLE PRECISION RHO 29* .. 30* .. Array Arguments .. 31* INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), 32* $ INDXQ( * ), PERM( * ) 33* DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), 34* $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) 35* .. 36* 37* 38*> \par Purpose: 39* ============= 40*> 41*> \verbatim 42*> 43*> DLAED8 merges the two sets of eigenvalues together into a single 44*> sorted set. Then it tries to deflate the size of the problem. 45*> There are two ways in which deflation can occur: when two or more 46*> eigenvalues are close together or if there is a tiny element in the 47*> Z vector. For each such occurrence the order of the related secular 48*> equation problem is reduced by one. 49*> \endverbatim 50* 51* Arguments: 52* ========== 53* 54*> \param[in] ICOMPQ 55*> \verbatim 56*> ICOMPQ is INTEGER 57*> = 0: Compute eigenvalues only. 58*> = 1: Compute eigenvectors of original dense symmetric matrix 59*> also. On entry, Q contains the orthogonal matrix used 60*> to reduce the original matrix to tridiagonal form. 61*> \endverbatim 62*> 63*> \param[out] K 64*> \verbatim 65*> K is INTEGER 66*> The number of non-deflated eigenvalues, and the order of the 67*> related secular equation. 68*> \endverbatim 69*> 70*> \param[in] N 71*> \verbatim 72*> N is INTEGER 73*> The dimension of the symmetric tridiagonal matrix. N >= 0. 74*> \endverbatim 75*> 76*> \param[in] QSIZ 77*> \verbatim 78*> QSIZ is INTEGER 79*> The dimension of the orthogonal matrix used to reduce 80*> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. 81*> \endverbatim 82*> 83*> \param[in,out] D 84*> \verbatim 85*> D is DOUBLE PRECISION array, dimension (N) 86*> On entry, the eigenvalues of the two submatrices to be 87*> combined. On exit, the trailing (N-K) updated eigenvalues 88*> (those which were deflated) sorted into increasing order. 89*> \endverbatim 90*> 91*> \param[in,out] Q 92*> \verbatim 93*> Q is DOUBLE PRECISION array, dimension (LDQ,N) 94*> If ICOMPQ = 0, Q is not referenced. Otherwise, 95*> on entry, Q contains the eigenvectors of the partially solved 96*> system which has been previously updated in matrix 97*> multiplies with other partially solved eigensystems. 98*> On exit, Q contains the trailing (N-K) updated eigenvectors 99*> (those which were deflated) in its last N-K columns. 100*> \endverbatim 101*> 102*> \param[in] LDQ 103*> \verbatim 104*> LDQ is INTEGER 105*> The leading dimension of the array Q. LDQ >= max(1,N). 106*> \endverbatim 107*> 108*> \param[in] INDXQ 109*> \verbatim 110*> INDXQ is INTEGER array, dimension (N) 111*> The permutation which separately sorts the two sub-problems 112*> in D into ascending order. Note that elements in the second 113*> half of this permutation must first have CUTPNT added to 114*> their values in order to be accurate. 115*> \endverbatim 116*> 117*> \param[in,out] RHO 118*> \verbatim 119*> RHO is DOUBLE PRECISION 120*> On entry, the off-diagonal element associated with the rank-1 121*> cut which originally split the two submatrices which are now 122*> being recombined. 123*> On exit, RHO has been modified to the value required by 124*> DLAED3. 125*> \endverbatim 126*> 127*> \param[in] CUTPNT 128*> \verbatim 129*> CUTPNT is INTEGER 130*> The location of the last eigenvalue in the leading 131*> sub-matrix. min(1,N) <= CUTPNT <= N. 132*> \endverbatim 133*> 134*> \param[in] Z 135*> \verbatim 136*> Z is DOUBLE PRECISION array, dimension (N) 137*> On entry, Z contains the updating vector (the last row of 138*> the first sub-eigenvector matrix and the first row of the 139*> second sub-eigenvector matrix). 140*> On exit, the contents of Z are destroyed by the updating 141*> process. 142*> \endverbatim 143*> 144*> \param[out] DLAMDA 145*> \verbatim 146*> DLAMDA is DOUBLE PRECISION array, dimension (N) 147*> A copy of the first K eigenvalues which will be used by 148*> DLAED3 to form the secular equation. 149*> \endverbatim 150*> 151*> \param[out] Q2 152*> \verbatim 153*> Q2 is DOUBLE PRECISION array, dimension (LDQ2,N) 154*> If ICOMPQ = 0, Q2 is not referenced. Otherwise, 155*> a copy of the first K eigenvectors which will be used by 156*> DLAED7 in a matrix multiply (DGEMM) to update the new 157*> eigenvectors. 158*> \endverbatim 159*> 160*> \param[in] LDQ2 161*> \verbatim 162*> LDQ2 is INTEGER 163*> The leading dimension of the array Q2. LDQ2 >= max(1,N). 164*> \endverbatim 165*> 166*> \param[out] W 167*> \verbatim 168*> W is DOUBLE PRECISION array, dimension (N) 169*> The first k values of the final deflation-altered z-vector and 170*> will be passed to DLAED3. 171*> \endverbatim 172*> 173*> \param[out] PERM 174*> \verbatim 175*> PERM is INTEGER array, dimension (N) 176*> The permutations (from deflation and sorting) to be applied 177*> to each eigenblock. 178*> \endverbatim 179*> 180*> \param[out] GIVPTR 181*> \verbatim 182*> GIVPTR is INTEGER 183*> The number of Givens rotations which took place in this 184*> subproblem. 185*> \endverbatim 186*> 187*> \param[out] GIVCOL 188*> \verbatim 189*> GIVCOL is INTEGER array, dimension (2, N) 190*> Each pair of numbers indicates a pair of columns to take place 191*> in a Givens rotation. 192*> \endverbatim 193*> 194*> \param[out] GIVNUM 195*> \verbatim 196*> GIVNUM is DOUBLE PRECISION array, dimension (2, N) 197*> Each number indicates the S value to be used in the 198*> corresponding Givens rotation. 199*> \endverbatim 200*> 201*> \param[out] INDXP 202*> \verbatim 203*> INDXP is INTEGER array, dimension (N) 204*> The permutation used to place deflated values of D at the end 205*> of the array. INDXP(1:K) points to the nondeflated D-values 206*> and INDXP(K+1:N) points to the deflated eigenvalues. 207*> \endverbatim 208*> 209*> \param[out] INDX 210*> \verbatim 211*> INDX is INTEGER array, dimension (N) 212*> The permutation used to sort the contents of D into ascending 213*> order. 214*> \endverbatim 215*> 216*> \param[out] INFO 217*> \verbatim 218*> INFO is INTEGER 219*> = 0: successful exit. 220*> < 0: if INFO = -i, the i-th argument had an illegal value. 221*> \endverbatim 222* 223* Authors: 224* ======== 225* 226*> \author Univ. of Tennessee 227*> \author Univ. of California Berkeley 228*> \author Univ. of Colorado Denver 229*> \author NAG Ltd. 230* 231*> \ingroup auxOTHERcomputational 232* 233*> \par Contributors: 234* ================== 235*> 236*> Jeff Rutter, Computer Science Division, University of California 237*> at Berkeley, USA 238* 239* ===================================================================== 240 SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, 241 $ CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, 242 $ GIVCOL, GIVNUM, INDXP, INDX, INFO ) 243* 244* -- LAPACK computational routine -- 245* -- LAPACK is a software package provided by Univ. of Tennessee, -- 246* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 247* 248* .. Scalar Arguments .. 249 INTEGER CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N, 250 $ QSIZ 251 DOUBLE PRECISION RHO 252* .. 253* .. Array Arguments .. 254 INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), 255 $ INDXQ( * ), PERM( * ) 256 DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), 257 $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) 258* .. 259* 260* ===================================================================== 261* 262* .. Parameters .. 263 DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT 264 PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0, 265 $ TWO = 2.0D0, EIGHT = 8.0D0 ) 266* .. 267* .. Local Scalars .. 268* 269 INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2 270 DOUBLE PRECISION C, EPS, S, T, TAU, TOL 271* .. 272* .. External Functions .. 273 INTEGER IDAMAX 274 DOUBLE PRECISION DLAMCH, DLAPY2 275 EXTERNAL IDAMAX, DLAMCH, DLAPY2 276* .. 277* .. External Subroutines .. 278 EXTERNAL DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA 279* .. 280* .. Intrinsic Functions .. 281 INTRINSIC ABS, MAX, MIN, SQRT 282* .. 283* .. Executable Statements .. 284* 285* Test the input parameters. 286* 287 INFO = 0 288* 289 IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN 290 INFO = -1 291 ELSE IF( N.LT.0 ) THEN 292 INFO = -3 293 ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN 294 INFO = -4 295 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 296 INFO = -7 297 ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN 298 INFO = -10 299 ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN 300 INFO = -14 301 END IF 302 IF( INFO.NE.0 ) THEN 303 CALL XERBLA( 'DLAED8', -INFO ) 304 RETURN 305 END IF 306* 307* Need to initialize GIVPTR to O here in case of quick exit 308* to prevent an unspecified code behavior (usually sigfault) 309* when IWORK array on entry to *stedc is not zeroed 310* (or at least some IWORK entries which used in *laed7 for GIVPTR). 311* 312 GIVPTR = 0 313* 314* Quick return if possible 315* 316 IF( N.EQ.0 ) 317 $ RETURN 318* 319 N1 = CUTPNT 320 N2 = N - N1 321 N1P1 = N1 + 1 322* 323 IF( RHO.LT.ZERO ) THEN 324 CALL DSCAL( N2, MONE, Z( N1P1 ), 1 ) 325 END IF 326* 327* Normalize z so that norm(z) = 1 328* 329 T = ONE / SQRT( TWO ) 330 DO 10 J = 1, N 331 INDX( J ) = J 332 10 CONTINUE 333 CALL DSCAL( N, T, Z, 1 ) 334 RHO = ABS( TWO*RHO ) 335* 336* Sort the eigenvalues into increasing order 337* 338 DO 20 I = CUTPNT + 1, N 339 INDXQ( I ) = INDXQ( I ) + CUTPNT 340 20 CONTINUE 341 DO 30 I = 1, N 342 DLAMDA( I ) = D( INDXQ( I ) ) 343 W( I ) = Z( INDXQ( I ) ) 344 30 CONTINUE 345 I = 1 346 J = CUTPNT + 1 347 CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX ) 348 DO 40 I = 1, N 349 D( I ) = DLAMDA( INDX( I ) ) 350 Z( I ) = W( INDX( I ) ) 351 40 CONTINUE 352* 353* Calculate the allowable deflation tolerance 354* 355 IMAX = IDAMAX( N, Z, 1 ) 356 JMAX = IDAMAX( N, D, 1 ) 357 EPS = DLAMCH( 'Epsilon' ) 358 TOL = EIGHT*EPS*ABS( D( JMAX ) ) 359* 360* If the rank-1 modifier is small enough, no more needs to be done 361* except to reorganize Q so that its columns correspond with the 362* elements in D. 363* 364 IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN 365 K = 0 366 IF( ICOMPQ.EQ.0 ) THEN 367 DO 50 J = 1, N 368 PERM( J ) = INDXQ( INDX( J ) ) 369 50 CONTINUE 370 ELSE 371 DO 60 J = 1, N 372 PERM( J ) = INDXQ( INDX( J ) ) 373 CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) 374 60 CONTINUE 375 CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ), 376 $ LDQ ) 377 END IF 378 RETURN 379 END IF 380* 381* If there are multiple eigenvalues then the problem deflates. Here 382* the number of equal eigenvalues are found. As each equal 383* eigenvalue is found, an elementary reflector is computed to rotate 384* the corresponding eigensubspace so that the corresponding 385* components of Z are zero in this new basis. 386* 387 K = 0 388 K2 = N + 1 389 DO 70 J = 1, N 390 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN 391* 392* Deflate due to small z component. 393* 394 K2 = K2 - 1 395 INDXP( K2 ) = J 396 IF( J.EQ.N ) 397 $ GO TO 110 398 ELSE 399 JLAM = J 400 GO TO 80 401 END IF 402 70 CONTINUE 403 80 CONTINUE 404 J = J + 1 405 IF( J.GT.N ) 406 $ GO TO 100 407 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN 408* 409* Deflate due to small z component. 410* 411 K2 = K2 - 1 412 INDXP( K2 ) = J 413 ELSE 414* 415* Check if eigenvalues are close enough to allow deflation. 416* 417 S = Z( JLAM ) 418 C = Z( J ) 419* 420* Find sqrt(a**2+b**2) without overflow or 421* destructive underflow. 422* 423 TAU = DLAPY2( C, S ) 424 T = D( J ) - D( JLAM ) 425 C = C / TAU 426 S = -S / TAU 427 IF( ABS( T*C*S ).LE.TOL ) THEN 428* 429* Deflation is possible. 430* 431 Z( J ) = TAU 432 Z( JLAM ) = ZERO 433* 434* Record the appropriate Givens rotation 435* 436 GIVPTR = GIVPTR + 1 437 GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) ) 438 GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) ) 439 GIVNUM( 1, GIVPTR ) = C 440 GIVNUM( 2, GIVPTR ) = S 441 IF( ICOMPQ.EQ.1 ) THEN 442 CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1, 443 $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S ) 444 END IF 445 T = D( JLAM )*C*C + D( J )*S*S 446 D( J ) = D( JLAM )*S*S + D( J )*C*C 447 D( JLAM ) = T 448 K2 = K2 - 1 449 I = 1 450 90 CONTINUE 451 IF( K2+I.LE.N ) THEN 452 IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN 453 INDXP( K2+I-1 ) = INDXP( K2+I ) 454 INDXP( K2+I ) = JLAM 455 I = I + 1 456 GO TO 90 457 ELSE 458 INDXP( K2+I-1 ) = JLAM 459 END IF 460 ELSE 461 INDXP( K2+I-1 ) = JLAM 462 END IF 463 JLAM = J 464 ELSE 465 K = K + 1 466 W( K ) = Z( JLAM ) 467 DLAMDA( K ) = D( JLAM ) 468 INDXP( K ) = JLAM 469 JLAM = J 470 END IF 471 END IF 472 GO TO 80 473 100 CONTINUE 474* 475* Record the last eigenvalue. 476* 477 K = K + 1 478 W( K ) = Z( JLAM ) 479 DLAMDA( K ) = D( JLAM ) 480 INDXP( K ) = JLAM 481* 482 110 CONTINUE 483* 484* Sort the eigenvalues and corresponding eigenvectors into DLAMDA 485* and Q2 respectively. The eigenvalues/vectors which were not 486* deflated go into the first K slots of DLAMDA and Q2 respectively, 487* while those which were deflated go into the last N - K slots. 488* 489 IF( ICOMPQ.EQ.0 ) THEN 490 DO 120 J = 1, N 491 JP = INDXP( J ) 492 DLAMDA( J ) = D( JP ) 493 PERM( J ) = INDXQ( INDX( JP ) ) 494 120 CONTINUE 495 ELSE 496 DO 130 J = 1, N 497 JP = INDXP( J ) 498 DLAMDA( J ) = D( JP ) 499 PERM( J ) = INDXQ( INDX( JP ) ) 500 CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) 501 130 CONTINUE 502 END IF 503* 504* The deflated eigenvalues and their corresponding vectors go back 505* into the last N - K slots of D and Q respectively. 506* 507 IF( K.LT.N ) THEN 508 IF( ICOMPQ.EQ.0 ) THEN 509 CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) 510 ELSE 511 CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) 512 CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, 513 $ Q( 1, K+1 ), LDQ ) 514 END IF 515 END IF 516* 517 RETURN 518* 519* End of DLAED8 520* 521 END 522