Annotation of rpl/lapack/lapack/dlaed8.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLAED8
! 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: *> \date November 2011
! 232: *
! 233: *> \ingroup auxOTHERcomputational
! 234: *
! 235: *> \par Contributors:
! 236: * ==================
! 237: *>
! 238: *> Jeff Rutter, Computer Science Division, University of California
! 239: *> at Berkeley, USA
! 240: *
! 241: * =====================================================================
1.1 bertrand 242: SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO,
243: $ CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR,
244: $ GIVCOL, GIVNUM, INDXP, INDX, INFO )
245: *
1.9 ! bertrand 246: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 247: * -- LAPACK is a software package provided by Univ. of Tennessee, --
248: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 249: * November 2011
1.1 bertrand 250: *
251: * .. Scalar Arguments ..
252: INTEGER CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N,
253: $ QSIZ
254: DOUBLE PRECISION RHO
255: * ..
256: * .. Array Arguments ..
257: INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
258: $ INDXQ( * ), PERM( * )
259: DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ),
260: $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * )
261: * ..
262: *
263: * =====================================================================
264: *
265: * .. Parameters ..
266: DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
267: PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
268: $ TWO = 2.0D0, EIGHT = 8.0D0 )
269: * ..
270: * .. Local Scalars ..
271: *
272: INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
273: DOUBLE PRECISION C, EPS, S, T, TAU, TOL
274: * ..
275: * .. External Functions ..
276: INTEGER IDAMAX
277: DOUBLE PRECISION DLAMCH, DLAPY2
278: EXTERNAL IDAMAX, DLAMCH, DLAPY2
279: * ..
280: * .. External Subroutines ..
281: EXTERNAL DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA
282: * ..
283: * .. Intrinsic Functions ..
284: INTRINSIC ABS, MAX, MIN, SQRT
285: * ..
286: * .. Executable Statements ..
287: *
288: * Test the input parameters.
289: *
290: INFO = 0
291: *
292: IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
293: INFO = -1
294: ELSE IF( N.LT.0 ) THEN
295: INFO = -3
296: ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
297: INFO = -4
298: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
299: INFO = -7
300: ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
301: INFO = -10
302: ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
303: INFO = -14
304: END IF
305: IF( INFO.NE.0 ) THEN
306: CALL XERBLA( 'DLAED8', -INFO )
307: RETURN
308: END IF
309: *
1.5 bertrand 310: * Need to initialize GIVPTR to O here in case of quick exit
311: * to prevent an unspecified code behavior (usually sigfault)
312: * when IWORK array on entry to *stedc is not zeroed
313: * (or at least some IWORK entries which used in *laed7 for GIVPTR).
314: *
315: GIVPTR = 0
316: *
1.1 bertrand 317: * Quick return if possible
318: *
319: IF( N.EQ.0 )
320: $ RETURN
321: *
322: N1 = CUTPNT
323: N2 = N - N1
324: N1P1 = N1 + 1
325: *
326: IF( RHO.LT.ZERO ) THEN
327: CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
328: END IF
329: *
330: * Normalize z so that norm(z) = 1
331: *
332: T = ONE / SQRT( TWO )
333: DO 10 J = 1, N
334: INDX( J ) = J
335: 10 CONTINUE
336: CALL DSCAL( N, T, Z, 1 )
337: RHO = ABS( TWO*RHO )
338: *
339: * Sort the eigenvalues into increasing order
340: *
341: DO 20 I = CUTPNT + 1, N
342: INDXQ( I ) = INDXQ( I ) + CUTPNT
343: 20 CONTINUE
344: DO 30 I = 1, N
345: DLAMDA( I ) = D( INDXQ( I ) )
346: W( I ) = Z( INDXQ( I ) )
347: 30 CONTINUE
348: I = 1
349: J = CUTPNT + 1
350: CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
351: DO 40 I = 1, N
352: D( I ) = DLAMDA( INDX( I ) )
353: Z( I ) = W( INDX( I ) )
354: 40 CONTINUE
355: *
356: * Calculate the allowable deflation tolerence
357: *
358: IMAX = IDAMAX( N, Z, 1 )
359: JMAX = IDAMAX( N, D, 1 )
360: EPS = DLAMCH( 'Epsilon' )
361: TOL = EIGHT*EPS*ABS( D( JMAX ) )
362: *
363: * If the rank-1 modifier is small enough, no more needs to be done
364: * except to reorganize Q so that its columns correspond with the
365: * elements in D.
366: *
367: IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
368: K = 0
369: IF( ICOMPQ.EQ.0 ) THEN
370: DO 50 J = 1, N
371: PERM( J ) = INDXQ( INDX( J ) )
372: 50 CONTINUE
373: ELSE
374: DO 60 J = 1, N
375: PERM( J ) = INDXQ( INDX( J ) )
376: CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
377: 60 CONTINUE
378: CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ),
379: $ LDQ )
380: END IF
381: RETURN
382: END IF
383: *
384: * If there are multiple eigenvalues then the problem deflates. Here
385: * the number of equal eigenvalues are found. As each equal
386: * eigenvalue is found, an elementary reflector is computed to rotate
387: * the corresponding eigensubspace so that the corresponding
388: * components of Z are zero in this new basis.
389: *
390: K = 0
391: K2 = N + 1
392: DO 70 J = 1, N
393: IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
394: *
395: * Deflate due to small z component.
396: *
397: K2 = K2 - 1
398: INDXP( K2 ) = J
399: IF( J.EQ.N )
400: $ GO TO 110
401: ELSE
402: JLAM = J
403: GO TO 80
404: END IF
405: 70 CONTINUE
406: 80 CONTINUE
407: J = J + 1
408: IF( J.GT.N )
409: $ GO TO 100
410: IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
411: *
412: * Deflate due to small z component.
413: *
414: K2 = K2 - 1
415: INDXP( K2 ) = J
416: ELSE
417: *
418: * Check if eigenvalues are close enough to allow deflation.
419: *
420: S = Z( JLAM )
421: C = Z( J )
422: *
423: * Find sqrt(a**2+b**2) without overflow or
424: * destructive underflow.
425: *
426: TAU = DLAPY2( C, S )
427: T = D( J ) - D( JLAM )
428: C = C / TAU
429: S = -S / TAU
430: IF( ABS( T*C*S ).LE.TOL ) THEN
431: *
432: * Deflation is possible.
433: *
434: Z( J ) = TAU
435: Z( JLAM ) = ZERO
436: *
437: * Record the appropriate Givens rotation
438: *
439: GIVPTR = GIVPTR + 1
440: GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
441: GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
442: GIVNUM( 1, GIVPTR ) = C
443: GIVNUM( 2, GIVPTR ) = S
444: IF( ICOMPQ.EQ.1 ) THEN
445: CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
446: $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
447: END IF
448: T = D( JLAM )*C*C + D( J )*S*S
449: D( J ) = D( JLAM )*S*S + D( J )*C*C
450: D( JLAM ) = T
451: K2 = K2 - 1
452: I = 1
453: 90 CONTINUE
454: IF( K2+I.LE.N ) THEN
455: IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
456: INDXP( K2+I-1 ) = INDXP( K2+I )
457: INDXP( K2+I ) = JLAM
458: I = I + 1
459: GO TO 90
460: ELSE
461: INDXP( K2+I-1 ) = JLAM
462: END IF
463: ELSE
464: INDXP( K2+I-1 ) = JLAM
465: END IF
466: JLAM = J
467: ELSE
468: K = K + 1
469: W( K ) = Z( JLAM )
470: DLAMDA( K ) = D( JLAM )
471: INDXP( K ) = JLAM
472: JLAM = J
473: END IF
474: END IF
475: GO TO 80
476: 100 CONTINUE
477: *
478: * Record the last eigenvalue.
479: *
480: K = K + 1
481: W( K ) = Z( JLAM )
482: DLAMDA( K ) = D( JLAM )
483: INDXP( K ) = JLAM
484: *
485: 110 CONTINUE
486: *
487: * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
488: * and Q2 respectively. The eigenvalues/vectors which were not
489: * deflated go into the first K slots of DLAMDA and Q2 respectively,
490: * while those which were deflated go into the last N - K slots.
491: *
492: IF( ICOMPQ.EQ.0 ) THEN
493: DO 120 J = 1, N
494: JP = INDXP( J )
495: DLAMDA( J ) = D( JP )
496: PERM( J ) = INDXQ( INDX( JP ) )
497: 120 CONTINUE
498: ELSE
499: DO 130 J = 1, N
500: JP = INDXP( J )
501: DLAMDA( J ) = D( JP )
502: PERM( J ) = INDXQ( INDX( JP ) )
503: CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
504: 130 CONTINUE
505: END IF
506: *
507: * The deflated eigenvalues and their corresponding vectors go back
508: * into the last N - K slots of D and Q respectively.
509: *
510: IF( K.LT.N ) THEN
511: IF( ICOMPQ.EQ.0 ) THEN
512: CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
513: ELSE
514: CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
515: CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2,
516: $ Q( 1, K+1 ), LDQ )
517: END IF
518: END IF
519: *
520: RETURN
521: *
522: * End of DLAED8
523: *
524: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed8.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack