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