Annotation of rpl/lapack/lapack/zlaed7.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZLAED7
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZLAED7 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaed7.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaed7.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaed7.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
! 22: * LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
! 23: * GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
! 24: * INFO )
! 25: *
! 26: * .. Scalar Arguments ..
! 27: * INTEGER CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
! 28: * $ TLVLS
! 29: * DOUBLE PRECISION RHO
! 30: * ..
! 31: * .. Array Arguments ..
! 32: * INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
! 33: * $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
! 34: * DOUBLE PRECISION D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
! 35: * COMPLEX*16 Q( LDQ, * ), WORK( * )
! 36: * ..
! 37: *
! 38: *
! 39: *> \par Purpose:
! 40: * =============
! 41: *>
! 42: *> \verbatim
! 43: *>
! 44: *> ZLAED7 computes the updated eigensystem of a diagonal
! 45: *> matrix after modification by a rank-one symmetric matrix. This
! 46: *> routine is used only for the eigenproblem which requires all
! 47: *> eigenvalues and optionally eigenvectors of a dense or banded
! 48: *> Hermitian matrix that has been reduced to tridiagonal form.
! 49: *>
! 50: *> T = Q(in) ( D(in) + RHO * Z*Z**H ) Q**H(in) = Q(out) * D(out) * Q**H(out)
! 51: *>
! 52: *> where Z = Q**Hu, u is a vector of length N with ones in the
! 53: *> CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
! 54: *>
! 55: *> The eigenvectors of the original matrix are stored in Q, and the
! 56: *> eigenvalues are in D. The algorithm consists of three stages:
! 57: *>
! 58: *> The first stage consists of deflating the size of the problem
! 59: *> when there are multiple eigenvalues or if there is a zero in
! 60: *> the Z vector. For each such occurence the dimension of the
! 61: *> secular equation problem is reduced by one. This stage is
! 62: *> performed by the routine DLAED2.
! 63: *>
! 64: *> The second stage consists of calculating the updated
! 65: *> eigenvalues. This is done by finding the roots of the secular
! 66: *> equation via the routine DLAED4 (as called by SLAED3).
! 67: *> This routine also calculates the eigenvectors of the current
! 68: *> problem.
! 69: *>
! 70: *> The final stage consists of computing the updated eigenvectors
! 71: *> directly using the updated eigenvalues. The eigenvectors for
! 72: *> the current problem are multiplied with the eigenvectors from
! 73: *> the overall problem.
! 74: *> \endverbatim
! 75: *
! 76: * Arguments:
! 77: * ==========
! 78: *
! 79: *> \param[in] N
! 80: *> \verbatim
! 81: *> N is INTEGER
! 82: *> The dimension of the symmetric tridiagonal matrix. N >= 0.
! 83: *> \endverbatim
! 84: *>
! 85: *> \param[in] CUTPNT
! 86: *> \verbatim
! 87: *> CUTPNT is INTEGER
! 88: *> Contains the location of the last eigenvalue in the leading
! 89: *> sub-matrix. min(1,N) <= CUTPNT <= N.
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] QSIZ
! 93: *> \verbatim
! 94: *> QSIZ is INTEGER
! 95: *> The dimension of the unitary matrix used to reduce
! 96: *> the full matrix to tridiagonal form. QSIZ >= N.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] TLVLS
! 100: *> \verbatim
! 101: *> TLVLS is INTEGER
! 102: *> The total number of merging levels in the overall divide and
! 103: *> conquer tree.
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[in] CURLVL
! 107: *> \verbatim
! 108: *> CURLVL is INTEGER
! 109: *> The current level in the overall merge routine,
! 110: *> 0 <= curlvl <= tlvls.
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[in] CURPBM
! 114: *> \verbatim
! 115: *> CURPBM is INTEGER
! 116: *> The current problem in the current level in the overall
! 117: *> merge routine (counting from upper left to lower right).
! 118: *> \endverbatim
! 119: *>
! 120: *> \param[in,out] D
! 121: *> \verbatim
! 122: *> D is DOUBLE PRECISION array, dimension (N)
! 123: *> On entry, the eigenvalues of the rank-1-perturbed matrix.
! 124: *> On exit, the eigenvalues of the repaired matrix.
! 125: *> \endverbatim
! 126: *>
! 127: *> \param[in,out] Q
! 128: *> \verbatim
! 129: *> Q is COMPLEX*16 array, dimension (LDQ,N)
! 130: *> On entry, the eigenvectors of the rank-1-perturbed matrix.
! 131: *> On exit, the eigenvectors of the repaired tridiagonal matrix.
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[in] LDQ
! 135: *> \verbatim
! 136: *> LDQ is INTEGER
! 137: *> The leading dimension of the array Q. LDQ >= max(1,N).
! 138: *> \endverbatim
! 139: *>
! 140: *> \param[in] RHO
! 141: *> \verbatim
! 142: *> RHO is DOUBLE PRECISION
! 143: *> Contains the subdiagonal element used to create the rank-1
! 144: *> modification.
! 145: *> \endverbatim
! 146: *>
! 147: *> \param[out] INDXQ
! 148: *> \verbatim
! 149: *> INDXQ is INTEGER array, dimension (N)
! 150: *> This contains the permutation which will reintegrate the
! 151: *> subproblem just solved back into sorted order,
! 152: *> ie. D( INDXQ( I = 1, N ) ) will be in ascending order.
! 153: *> \endverbatim
! 154: *>
! 155: *> \param[out] IWORK
! 156: *> \verbatim
! 157: *> IWORK is INTEGER array, dimension (4*N)
! 158: *> \endverbatim
! 159: *>
! 160: *> \param[out] RWORK
! 161: *> \verbatim
! 162: *> RWORK is DOUBLE PRECISION array,
! 163: *> dimension (3*N+2*QSIZ*N)
! 164: *> \endverbatim
! 165: *>
! 166: *> \param[out] WORK
! 167: *> \verbatim
! 168: *> WORK is COMPLEX*16 array, dimension (QSIZ*N)
! 169: *> \endverbatim
! 170: *>
! 171: *> \param[in,out] QSTORE
! 172: *> \verbatim
! 173: *> QSTORE is DOUBLE PRECISION array, dimension (N**2+1)
! 174: *> Stores eigenvectors of submatrices encountered during
! 175: *> divide and conquer, packed together. QPTR points to
! 176: *> beginning of the submatrices.
! 177: *> \endverbatim
! 178: *>
! 179: *> \param[in,out] QPTR
! 180: *> \verbatim
! 181: *> QPTR is INTEGER array, dimension (N+2)
! 182: *> List of indices pointing to beginning of submatrices stored
! 183: *> in QSTORE. The submatrices are numbered starting at the
! 184: *> bottom left of the divide and conquer tree, from left to
! 185: *> right and bottom to top.
! 186: *> \endverbatim
! 187: *>
! 188: *> \param[in] PRMPTR
! 189: *> \verbatim
! 190: *> PRMPTR is INTEGER array, dimension (N lg N)
! 191: *> Contains a list of pointers which indicate where in PERM a
! 192: *> level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
! 193: *> indicates the size of the permutation and also the size of
! 194: *> the full, non-deflated problem.
! 195: *> \endverbatim
! 196: *>
! 197: *> \param[in] PERM
! 198: *> \verbatim
! 199: *> PERM is INTEGER array, dimension (N lg N)
! 200: *> Contains the permutations (from deflation and sorting) to be
! 201: *> applied to each eigenblock.
! 202: *> \endverbatim
! 203: *>
! 204: *> \param[in] GIVPTR
! 205: *> \verbatim
! 206: *> GIVPTR is INTEGER array, dimension (N lg N)
! 207: *> Contains a list of pointers which indicate where in GIVCOL a
! 208: *> level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
! 209: *> indicates the number of Givens rotations.
! 210: *> \endverbatim
! 211: *>
! 212: *> \param[in] GIVCOL
! 213: *> \verbatim
! 214: *> GIVCOL is INTEGER array, dimension (2, N lg N)
! 215: *> Each pair of numbers indicates a pair of columns to take place
! 216: *> in a Givens rotation.
! 217: *> \endverbatim
! 218: *>
! 219: *> \param[in] GIVNUM
! 220: *> \verbatim
! 221: *> GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)
! 222: *> Each number indicates the S value to be used in the
! 223: *> corresponding Givens rotation.
! 224: *> \endverbatim
! 225: *>
! 226: *> \param[out] INFO
! 227: *> \verbatim
! 228: *> INFO is INTEGER
! 229: *> = 0: successful exit.
! 230: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 231: *> > 0: if INFO = 1, an eigenvalue did not converge
! 232: *> \endverbatim
! 233: *
! 234: * Authors:
! 235: * ========
! 236: *
! 237: *> \author Univ. of Tennessee
! 238: *> \author Univ. of California Berkeley
! 239: *> \author Univ. of Colorado Denver
! 240: *> \author NAG Ltd.
! 241: *
! 242: *> \date November 2011
! 243: *
! 244: *> \ingroup complex16OTHERcomputational
! 245: *
! 246: * =====================================================================
1.1 bertrand 247: SUBROUTINE ZLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
248: $ LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
249: $ GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
250: $ INFO )
251: *
1.9 ! bertrand 252: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 253: * -- LAPACK is a software package provided by Univ. of Tennessee, --
254: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 255: * November 2011
1.1 bertrand 256: *
257: * .. Scalar Arguments ..
258: INTEGER CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
259: $ TLVLS
260: DOUBLE PRECISION RHO
261: * ..
262: * .. Array Arguments ..
263: INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
264: $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
265: DOUBLE PRECISION D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
266: COMPLEX*16 Q( LDQ, * ), WORK( * )
267: * ..
268: *
269: * =====================================================================
270: *
271: * .. Local Scalars ..
272: INTEGER COLTYP, CURR, I, IDLMDA, INDX,
273: $ INDXC, INDXP, IQ, IW, IZ, K, N1, N2, PTR
274: * ..
275: * .. External Subroutines ..
276: EXTERNAL DLAED9, DLAEDA, DLAMRG, XERBLA, ZLACRM, ZLAED8
277: * ..
278: * .. Intrinsic Functions ..
279: INTRINSIC MAX, MIN
280: * ..
281: * .. Executable Statements ..
282: *
283: * Test the input parameters.
284: *
285: INFO = 0
286: *
287: * IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
288: * INFO = -1
289: * ELSE IF( N.LT.0 ) THEN
290: IF( N.LT.0 ) THEN
291: INFO = -1
292: ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
293: INFO = -2
294: ELSE IF( QSIZ.LT.N ) THEN
295: INFO = -3
296: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
297: INFO = -9
298: END IF
299: IF( INFO.NE.0 ) THEN
300: CALL XERBLA( 'ZLAED7', -INFO )
301: RETURN
302: END IF
303: *
304: * Quick return if possible
305: *
306: IF( N.EQ.0 )
307: $ RETURN
308: *
309: * The following values are for bookkeeping purposes only. They are
310: * integer pointers which indicate the portion of the workspace
311: * used by a particular array in DLAED2 and SLAED3.
312: *
313: IZ = 1
314: IDLMDA = IZ + N
315: IW = IDLMDA + N
316: IQ = IW + N
317: *
318: INDX = 1
319: INDXC = INDX + N
320: COLTYP = INDXC + N
321: INDXP = COLTYP + N
322: *
323: * Form the z-vector which consists of the last row of Q_1 and the
324: * first row of Q_2.
325: *
326: PTR = 1 + 2**TLVLS
327: DO 10 I = 1, CURLVL - 1
328: PTR = PTR + 2**( TLVLS-I )
329: 10 CONTINUE
330: CURR = PTR + CURPBM
331: CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
332: $ GIVCOL, GIVNUM, QSTORE, QPTR, RWORK( IZ ),
333: $ RWORK( IZ+N ), INFO )
334: *
335: * When solving the final problem, we no longer need the stored data,
336: * so we will overwrite the data from this level onto the previously
337: * used storage space.
338: *
339: IF( CURLVL.EQ.TLVLS ) THEN
340: QPTR( CURR ) = 1
341: PRMPTR( CURR ) = 1
342: GIVPTR( CURR ) = 1
343: END IF
344: *
345: * Sort and Deflate eigenvalues.
346: *
347: CALL ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, RWORK( IZ ),
348: $ RWORK( IDLMDA ), WORK, QSIZ, RWORK( IW ),
349: $ IWORK( INDXP ), IWORK( INDX ), INDXQ,
350: $ PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
351: $ GIVCOL( 1, GIVPTR( CURR ) ),
352: $ GIVNUM( 1, GIVPTR( CURR ) ), INFO )
353: PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
354: GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
355: *
356: * Solve Secular Equation.
357: *
358: IF( K.NE.0 ) THEN
359: CALL DLAED9( K, 1, K, N, D, RWORK( IQ ), K, RHO,
360: $ RWORK( IDLMDA ), RWORK( IW ),
361: $ QSTORE( QPTR( CURR ) ), K, INFO )
362: CALL ZLACRM( QSIZ, K, WORK, QSIZ, QSTORE( QPTR( CURR ) ), K, Q,
363: $ LDQ, RWORK( IQ ) )
364: QPTR( CURR+1 ) = QPTR( CURR ) + K**2
365: IF( INFO.NE.0 ) THEN
366: RETURN
367: END IF
368: *
369: * Prepare the INDXQ sorting premutation.
370: *
371: N1 = K
372: N2 = N - K
373: CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
374: ELSE
375: QPTR( CURR+1 ) = QPTR( CURR )
376: DO 20 I = 1, N
377: INDXQ( I ) = I
378: 20 CONTINUE
379: END IF
380: *
381: RETURN
382: *
383: * End of ZLAED7
384: *
385: END
CVSweb interface <joel.bertrand@systella.fr>