001: SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, 002: $ CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, 003: $ GIVCOL, GIVNUM, INDXP, INDX, INFO ) 004: * 005: * -- LAPACK routine (version 3.2) -- 006: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 007: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 008: * November 2006 009: * 010: * .. Scalar Arguments .. 011: INTEGER CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N, 012: $ QSIZ 013: DOUBLE PRECISION RHO 014: * .. 015: * .. Array Arguments .. 016: INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), 017: $ INDXQ( * ), PERM( * ) 018: DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), 019: $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) 020: * .. 021: * 022: * Purpose 023: * ======= 024: * 025: * DLAED8 merges the two sets of eigenvalues together into a single 026: * sorted set. Then it tries to deflate the size of the problem. 027: * There are two ways in which deflation can occur: when two or more 028: * eigenvalues are close together or if there is a tiny element in the 029: * Z vector. For each such occurrence the order of the related secular 030: * equation problem is reduced by one. 031: * 032: * Arguments 033: * ========= 034: * 035: * ICOMPQ (input) INTEGER 036: * = 0: Compute eigenvalues only. 037: * = 1: Compute eigenvectors of original dense symmetric matrix 038: * also. On entry, Q contains the orthogonal matrix used 039: * to reduce the original matrix to tridiagonal form. 040: * 041: * K (output) INTEGER 042: * The number of non-deflated eigenvalues, and the order of the 043: * related secular equation. 044: * 045: * N (input) INTEGER 046: * The dimension of the symmetric tridiagonal matrix. N >= 0. 047: * 048: * QSIZ (input) INTEGER 049: * The dimension of the orthogonal matrix used to reduce 050: * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. 051: * 052: * D (input/output) DOUBLE PRECISION array, dimension (N) 053: * On entry, the eigenvalues of the two submatrices to be 054: * combined. On exit, the trailing (N-K) updated eigenvalues 055: * (those which were deflated) sorted into increasing order. 056: * 057: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) 058: * If ICOMPQ = 0, Q is not referenced. Otherwise, 059: * on entry, Q contains the eigenvectors of the partially solved 060: * system which has been previously updated in matrix 061: * multiplies with other partially solved eigensystems. 062: * On exit, Q contains the trailing (N-K) updated eigenvectors 063: * (those which were deflated) in its last N-K columns. 064: * 065: * LDQ (input) INTEGER 066: * The leading dimension of the array Q. LDQ >= max(1,N). 067: * 068: * INDXQ (input) INTEGER array, dimension (N) 069: * The permutation which separately sorts the two sub-problems 070: * in D into ascending order. Note that elements in the second 071: * half of this permutation must first have CUTPNT added to 072: * their values in order to be accurate. 073: * 074: * RHO (input/output) DOUBLE PRECISION 075: * On entry, the off-diagonal element associated with the rank-1 076: * cut which originally split the two submatrices which are now 077: * being recombined. 078: * On exit, RHO has been modified to the value required by 079: * DLAED3. 080: * 081: * CUTPNT (input) INTEGER 082: * The location of the last eigenvalue in the leading 083: * sub-matrix. min(1,N) <= CUTPNT <= N. 084: * 085: * Z (input) DOUBLE PRECISION array, dimension (N) 086: * On entry, Z contains the updating vector (the last row of 087: * the first sub-eigenvector matrix and the first row of the 088: * second sub-eigenvector matrix). 089: * On exit, the contents of Z are destroyed by the updating 090: * process. 091: * 092: * DLAMDA (output) DOUBLE PRECISION array, dimension (N) 093: * A copy of the first K eigenvalues which will be used by 094: * DLAED3 to form the secular equation. 095: * 096: * Q2 (output) DOUBLE PRECISION array, dimension (LDQ2,N) 097: * If ICOMPQ = 0, Q2 is not referenced. Otherwise, 098: * a copy of the first K eigenvectors which will be used by 099: * DLAED7 in a matrix multiply (DGEMM) to update the new 100: * eigenvectors. 101: * 102: * LDQ2 (input) INTEGER 103: * The leading dimension of the array Q2. LDQ2 >= max(1,N). 104: * 105: * W (output) DOUBLE PRECISION array, dimension (N) 106: * The first k values of the final deflation-altered z-vector and 107: * will be passed to DLAED3. 108: * 109: * PERM (output) INTEGER array, dimension (N) 110: * The permutations (from deflation and sorting) to be applied 111: * to each eigenblock. 112: * 113: * GIVPTR (output) INTEGER 114: * The number of Givens rotations which took place in this 115: * subproblem. 116: * 117: * GIVCOL (output) INTEGER array, dimension (2, N) 118: * Each pair of numbers indicates a pair of columns to take place 119: * in a Givens rotation. 120: * 121: * GIVNUM (output) DOUBLE PRECISION array, dimension (2, N) 122: * Each number indicates the S value to be used in the 123: * corresponding Givens rotation. 124: * 125: * INDXP (workspace) INTEGER array, dimension (N) 126: * The permutation used to place deflated values of D at the end 127: * of the array. INDXP(1:K) points to the nondeflated D-values 128: * and INDXP(K+1:N) points to the deflated eigenvalues. 129: * 130: * INDX (workspace) INTEGER array, dimension (N) 131: * The permutation used to sort the contents of D into ascending 132: * order. 133: * 134: * INFO (output) INTEGER 135: * = 0: successful exit. 136: * < 0: if INFO = -i, the i-th argument had an illegal value. 137: * 138: * Further Details 139: * =============== 140: * 141: * Based on contributions by 142: * Jeff Rutter, Computer Science Division, University of California 143: * at Berkeley, USA 144: * 145: * ===================================================================== 146: * 147: * .. Parameters .. 148: DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT 149: PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0, 150: $ TWO = 2.0D0, EIGHT = 8.0D0 ) 151: * .. 152: * .. Local Scalars .. 153: * 154: INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2 155: DOUBLE PRECISION C, EPS, S, T, TAU, TOL 156: * .. 157: * .. External Functions .. 158: INTEGER IDAMAX 159: DOUBLE PRECISION DLAMCH, DLAPY2 160: EXTERNAL IDAMAX, DLAMCH, DLAPY2 161: * .. 162: * .. External Subroutines .. 163: EXTERNAL DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA 164: * .. 165: * .. Intrinsic Functions .. 166: INTRINSIC ABS, MAX, MIN, SQRT 167: * .. 168: * .. Executable Statements .. 169: * 170: * Test the input parameters. 171: * 172: INFO = 0 173: * 174: IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN 175: INFO = -1 176: ELSE IF( N.LT.0 ) THEN 177: INFO = -3 178: ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN 179: INFO = -4 180: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 181: INFO = -7 182: ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN 183: INFO = -10 184: ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN 185: INFO = -14 186: END IF 187: IF( INFO.NE.0 ) THEN 188: CALL XERBLA( 'DLAED8', -INFO ) 189: RETURN 190: END IF 191: * 192: * Quick return if possible 193: * 194: IF( N.EQ.0 ) 195: $ RETURN 196: * 197: N1 = CUTPNT 198: N2 = N - N1 199: N1P1 = N1 + 1 200: * 201: IF( RHO.LT.ZERO ) THEN 202: CALL DSCAL( N2, MONE, Z( N1P1 ), 1 ) 203: END IF 204: * 205: * Normalize z so that norm(z) = 1 206: * 207: T = ONE / SQRT( TWO ) 208: DO 10 J = 1, N 209: INDX( J ) = J 210: 10 CONTINUE 211: CALL DSCAL( N, T, Z, 1 ) 212: RHO = ABS( TWO*RHO ) 213: * 214: * Sort the eigenvalues into increasing order 215: * 216: DO 20 I = CUTPNT + 1, N 217: INDXQ( I ) = INDXQ( I ) + CUTPNT 218: 20 CONTINUE 219: DO 30 I = 1, N 220: DLAMDA( I ) = D( INDXQ( I ) ) 221: W( I ) = Z( INDXQ( I ) ) 222: 30 CONTINUE 223: I = 1 224: J = CUTPNT + 1 225: CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX ) 226: DO 40 I = 1, N 227: D( I ) = DLAMDA( INDX( I ) ) 228: Z( I ) = W( INDX( I ) ) 229: 40 CONTINUE 230: * 231: * Calculate the allowable deflation tolerence 232: * 233: IMAX = IDAMAX( N, Z, 1 ) 234: JMAX = IDAMAX( N, D, 1 ) 235: EPS = DLAMCH( 'Epsilon' ) 236: TOL = EIGHT*EPS*ABS( D( JMAX ) ) 237: * 238: * If the rank-1 modifier is small enough, no more needs to be done 239: * except to reorganize Q so that its columns correspond with the 240: * elements in D. 241: * 242: IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN 243: K = 0 244: IF( ICOMPQ.EQ.0 ) THEN 245: DO 50 J = 1, N 246: PERM( J ) = INDXQ( INDX( J ) ) 247: 50 CONTINUE 248: ELSE 249: DO 60 J = 1, N 250: PERM( J ) = INDXQ( INDX( J ) ) 251: CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) 252: 60 CONTINUE 253: CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ), 254: $ LDQ ) 255: END IF 256: RETURN 257: END IF 258: * 259: * If there are multiple eigenvalues then the problem deflates. Here 260: * the number of equal eigenvalues are found. As each equal 261: * eigenvalue is found, an elementary reflector is computed to rotate 262: * the corresponding eigensubspace so that the corresponding 263: * components of Z are zero in this new basis. 264: * 265: K = 0 266: GIVPTR = 0 267: K2 = N + 1 268: DO 70 J = 1, N 269: IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN 270: * 271: * Deflate due to small z component. 272: * 273: K2 = K2 - 1 274: INDXP( K2 ) = J 275: IF( J.EQ.N ) 276: $ GO TO 110 277: ELSE 278: JLAM = J 279: GO TO 80 280: END IF 281: 70 CONTINUE 282: 80 CONTINUE 283: J = J + 1 284: IF( J.GT.N ) 285: $ GO TO 100 286: IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN 287: * 288: * Deflate due to small z component. 289: * 290: K2 = K2 - 1 291: INDXP( K2 ) = J 292: ELSE 293: * 294: * Check if eigenvalues are close enough to allow deflation. 295: * 296: S = Z( JLAM ) 297: C = Z( J ) 298: * 299: * Find sqrt(a**2+b**2) without overflow or 300: * destructive underflow. 301: * 302: TAU = DLAPY2( C, S ) 303: T = D( J ) - D( JLAM ) 304: C = C / TAU 305: S = -S / TAU 306: IF( ABS( T*C*S ).LE.TOL ) THEN 307: * 308: * Deflation is possible. 309: * 310: Z( J ) = TAU 311: Z( JLAM ) = ZERO 312: * 313: * Record the appropriate Givens rotation 314: * 315: GIVPTR = GIVPTR + 1 316: GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) ) 317: GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) ) 318: GIVNUM( 1, GIVPTR ) = C 319: GIVNUM( 2, GIVPTR ) = S 320: IF( ICOMPQ.EQ.1 ) THEN 321: CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1, 322: $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S ) 323: END IF 324: T = D( JLAM )*C*C + D( J )*S*S 325: D( J ) = D( JLAM )*S*S + D( J )*C*C 326: D( JLAM ) = T 327: K2 = K2 - 1 328: I = 1 329: 90 CONTINUE 330: IF( K2+I.LE.N ) THEN 331: IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN 332: INDXP( K2+I-1 ) = INDXP( K2+I ) 333: INDXP( K2+I ) = JLAM 334: I = I + 1 335: GO TO 90 336: ELSE 337: INDXP( K2+I-1 ) = JLAM 338: END IF 339: ELSE 340: INDXP( K2+I-1 ) = JLAM 341: END IF 342: JLAM = J 343: ELSE 344: K = K + 1 345: W( K ) = Z( JLAM ) 346: DLAMDA( K ) = D( JLAM ) 347: INDXP( K ) = JLAM 348: JLAM = J 349: END IF 350: END IF 351: GO TO 80 352: 100 CONTINUE 353: * 354: * Record the last eigenvalue. 355: * 356: K = K + 1 357: W( K ) = Z( JLAM ) 358: DLAMDA( K ) = D( JLAM ) 359: INDXP( K ) = JLAM 360: * 361: 110 CONTINUE 362: * 363: * Sort the eigenvalues and corresponding eigenvectors into DLAMDA 364: * and Q2 respectively. The eigenvalues/vectors which were not 365: * deflated go into the first K slots of DLAMDA and Q2 respectively, 366: * while those which were deflated go into the last N - K slots. 367: * 368: IF( ICOMPQ.EQ.0 ) THEN 369: DO 120 J = 1, N 370: JP = INDXP( J ) 371: DLAMDA( J ) = D( JP ) 372: PERM( J ) = INDXQ( INDX( JP ) ) 373: 120 CONTINUE 374: ELSE 375: DO 130 J = 1, N 376: JP = INDXP( J ) 377: DLAMDA( J ) = D( JP ) 378: PERM( J ) = INDXQ( INDX( JP ) ) 379: CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) 380: 130 CONTINUE 381: END IF 382: * 383: * The deflated eigenvalues and their corresponding vectors go back 384: * into the last N - K slots of D and Q respectively. 385: * 386: IF( K.LT.N ) THEN 387: IF( ICOMPQ.EQ.0 ) THEN 388: CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) 389: ELSE 390: CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) 391: CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, 392: $ Q( 1, K+1 ), LDQ ) 393: END IF 394: END IF 395: * 396: RETURN 397: * 398: * End of DLAED8 399: * 400: END 401: