Annotation of rpl/lapack/lapack/zlaed0.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b ZLAED0
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZLAED0 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaed0.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaed0.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaed0.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
! 22: * IWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER INFO, LDQ, LDQS, N, QSIZ
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * INTEGER IWORK( * )
! 29: * DOUBLE PRECISION D( * ), E( * ), RWORK( * )
! 30: * COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> Using the divide and conquer method, ZLAED0 computes all eigenvalues
! 40: *> of a symmetric tridiagonal matrix which is one diagonal block of
! 41: *> those from reducing a dense or band Hermitian matrix and
! 42: *> corresponding eigenvectors of the dense or band matrix.
! 43: *> \endverbatim
! 44: *
! 45: * Arguments:
! 46: * ==========
! 47: *
! 48: *> \param[in] QSIZ
! 49: *> \verbatim
! 50: *> QSIZ is INTEGER
! 51: *> The dimension of the unitary matrix used to reduce
! 52: *> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
! 53: *> \endverbatim
! 54: *>
! 55: *> \param[in] N
! 56: *> \verbatim
! 57: *> N is INTEGER
! 58: *> The dimension of the symmetric tridiagonal matrix. N >= 0.
! 59: *> \endverbatim
! 60: *>
! 61: *> \param[in,out] D
! 62: *> \verbatim
! 63: *> D is DOUBLE PRECISION array, dimension (N)
! 64: *> On entry, the diagonal elements of the tridiagonal matrix.
! 65: *> On exit, the eigenvalues in ascending order.
! 66: *> \endverbatim
! 67: *>
! 68: *> \param[in,out] E
! 69: *> \verbatim
! 70: *> E is DOUBLE PRECISION array, dimension (N-1)
! 71: *> On entry, the off-diagonal elements of the tridiagonal matrix.
! 72: *> On exit, E has been destroyed.
! 73: *> \endverbatim
! 74: *>
! 75: *> \param[in,out] Q
! 76: *> \verbatim
! 77: *> Q is COMPLEX*16 array, dimension (LDQ,N)
! 78: *> On entry, Q must contain an QSIZ x N matrix whose columns
! 79: *> unitarily orthonormal. It is a part of the unitary matrix
! 80: *> that reduces the full dense Hermitian matrix to a
! 81: *> (reducible) symmetric tridiagonal matrix.
! 82: *> \endverbatim
! 83: *>
! 84: *> \param[in] LDQ
! 85: *> \verbatim
! 86: *> LDQ is INTEGER
! 87: *> The leading dimension of the array Q. LDQ >= max(1,N).
! 88: *> \endverbatim
! 89: *>
! 90: *> \param[out] IWORK
! 91: *> \verbatim
! 92: *> IWORK is INTEGER array,
! 93: *> the dimension of IWORK must be at least
! 94: *> 6 + 6*N + 5*N*lg N
! 95: *> ( lg( N ) = smallest integer k
! 96: *> such that 2^k >= N )
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[out] RWORK
! 100: *> \verbatim
! 101: *> RWORK is DOUBLE PRECISION array,
! 102: *> dimension (1 + 3*N + 2*N*lg N + 3*N**2)
! 103: *> ( lg( N ) = smallest integer k
! 104: *> such that 2^k >= N )
! 105: *> \endverbatim
! 106: *>
! 107: *> \param[out] QSTORE
! 108: *> \verbatim
! 109: *> QSTORE is COMPLEX*16 array, dimension (LDQS, N)
! 110: *> Used to store parts of
! 111: *> the eigenvector matrix when the updating matrix multiplies
! 112: *> take place.
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in] LDQS
! 116: *> \verbatim
! 117: *> LDQS is INTEGER
! 118: *> The leading dimension of the array QSTORE.
! 119: *> LDQS >= max(1,N).
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[out] INFO
! 123: *> \verbatim
! 124: *> INFO is INTEGER
! 125: *> = 0: successful exit.
! 126: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 127: *> > 0: The algorithm failed to compute an eigenvalue while
! 128: *> working on the submatrix lying in rows and columns
! 129: *> INFO/(N+1) through mod(INFO,N+1).
! 130: *> \endverbatim
! 131: *
! 132: * Authors:
! 133: * ========
! 134: *
! 135: *> \author Univ. of Tennessee
! 136: *> \author Univ. of California Berkeley
! 137: *> \author Univ. of Colorado Denver
! 138: *> \author NAG Ltd.
! 139: *
! 140: *> \date November 2011
! 141: *
! 142: *> \ingroup complex16OTHERcomputational
! 143: *
! 144: * =====================================================================
1.1 bertrand 145: SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
146: $ IWORK, INFO )
147: *
1.8 ! bertrand 148: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 149: * -- LAPACK is a software package provided by Univ. of Tennessee, --
150: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 151: * November 2011
1.1 bertrand 152: *
153: * .. Scalar Arguments ..
154: INTEGER INFO, LDQ, LDQS, N, QSIZ
155: * ..
156: * .. Array Arguments ..
157: INTEGER IWORK( * )
158: DOUBLE PRECISION D( * ), E( * ), RWORK( * )
159: COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * )
160: * ..
161: *
162: * =====================================================================
163: *
164: * Warning: N could be as big as QSIZ!
165: *
166: * .. Parameters ..
167: DOUBLE PRECISION TWO
168: PARAMETER ( TWO = 2.D+0 )
169: * ..
170: * .. Local Scalars ..
171: INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
172: $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
173: $ J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
174: $ SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
175: DOUBLE PRECISION TEMP
176: * ..
177: * .. External Subroutines ..
178: EXTERNAL DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
179: * ..
180: * .. External Functions ..
181: INTEGER ILAENV
182: EXTERNAL ILAENV
183: * ..
184: * .. Intrinsic Functions ..
185: INTRINSIC ABS, DBLE, INT, LOG, MAX
186: * ..
187: * .. Executable Statements ..
188: *
189: * Test the input parameters.
190: *
191: INFO = 0
192: *
193: * IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
194: * INFO = -1
195: * ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
196: * $ THEN
197: IF( QSIZ.LT.MAX( 0, N ) ) THEN
198: INFO = -1
199: ELSE IF( N.LT.0 ) THEN
200: INFO = -2
201: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
202: INFO = -6
203: ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
204: INFO = -8
205: END IF
206: IF( INFO.NE.0 ) THEN
207: CALL XERBLA( 'ZLAED0', -INFO )
208: RETURN
209: END IF
210: *
211: * Quick return if possible
212: *
213: IF( N.EQ.0 )
214: $ RETURN
215: *
216: SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
217: *
218: * Determine the size and placement of the submatrices, and save in
219: * the leading elements of IWORK.
220: *
221: IWORK( 1 ) = N
222: SUBPBS = 1
223: TLVLS = 0
224: 10 CONTINUE
225: IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
226: DO 20 J = SUBPBS, 1, -1
227: IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
228: IWORK( 2*J-1 ) = IWORK( J ) / 2
229: 20 CONTINUE
230: TLVLS = TLVLS + 1
231: SUBPBS = 2*SUBPBS
232: GO TO 10
233: END IF
234: DO 30 J = 2, SUBPBS
235: IWORK( J ) = IWORK( J ) + IWORK( J-1 )
236: 30 CONTINUE
237: *
238: * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
239: * using rank-1 modifications (cuts).
240: *
241: SPM1 = SUBPBS - 1
242: DO 40 I = 1, SPM1
243: SUBMAT = IWORK( I ) + 1
244: SMM1 = SUBMAT - 1
245: D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
246: D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
247: 40 CONTINUE
248: *
249: INDXQ = 4*N + 3
250: *
251: * Set up workspaces for eigenvalues only/accumulate new vectors
252: * routine
253: *
254: TEMP = LOG( DBLE( N ) ) / LOG( TWO )
255: LGN = INT( TEMP )
256: IF( 2**LGN.LT.N )
257: $ LGN = LGN + 1
258: IF( 2**LGN.LT.N )
259: $ LGN = LGN + 1
260: IPRMPT = INDXQ + N + 1
261: IPERM = IPRMPT + N*LGN
262: IQPTR = IPERM + N*LGN
263: IGIVPT = IQPTR + N + 2
264: IGIVCL = IGIVPT + N*LGN
265: *
266: IGIVNM = 1
267: IQ = IGIVNM + 2*N*LGN
268: IWREM = IQ + N**2 + 1
269: * Initialize pointers
270: DO 50 I = 0, SUBPBS
271: IWORK( IPRMPT+I ) = 1
272: IWORK( IGIVPT+I ) = 1
273: 50 CONTINUE
274: IWORK( IQPTR ) = 1
275: *
276: * Solve each submatrix eigenproblem at the bottom of the divide and
277: * conquer tree.
278: *
279: CURR = 0
280: DO 70 I = 0, SPM1
281: IF( I.EQ.0 ) THEN
282: SUBMAT = 1
283: MATSIZ = IWORK( 1 )
284: ELSE
285: SUBMAT = IWORK( I ) + 1
286: MATSIZ = IWORK( I+1 ) - IWORK( I )
287: END IF
288: LL = IQ - 1 + IWORK( IQPTR+CURR )
289: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
290: $ RWORK( LL ), MATSIZ, RWORK, INFO )
291: CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
292: $ MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
293: $ RWORK( IWREM ) )
294: IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
295: CURR = CURR + 1
296: IF( INFO.GT.0 ) THEN
297: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
298: RETURN
299: END IF
300: K = 1
301: DO 60 J = SUBMAT, IWORK( I+1 )
302: IWORK( INDXQ+J ) = K
303: K = K + 1
304: 60 CONTINUE
305: 70 CONTINUE
306: *
307: * Successively merge eigensystems of adjacent submatrices
308: * into eigensystem for the corresponding larger matrix.
309: *
310: * while ( SUBPBS > 1 )
311: *
312: CURLVL = 1
313: 80 CONTINUE
314: IF( SUBPBS.GT.1 ) THEN
315: SPM2 = SUBPBS - 2
316: DO 90 I = 0, SPM2, 2
317: IF( I.EQ.0 ) THEN
318: SUBMAT = 1
319: MATSIZ = IWORK( 2 )
320: MSD2 = IWORK( 1 )
321: CURPRB = 0
322: ELSE
323: SUBMAT = IWORK( I ) + 1
324: MATSIZ = IWORK( I+2 ) - IWORK( I )
325: MSD2 = MATSIZ / 2
326: CURPRB = CURPRB + 1
327: END IF
328: *
329: * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
330: * into an eigensystem of size MATSIZ. ZLAED7 handles the case
331: * when the eigenvectors of a full or band Hermitian matrix (which
332: * was reduced to tridiagonal form) are desired.
333: *
334: * I am free to use Q as a valuable working space until Loop 150.
335: *
336: CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
337: $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
338: $ E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
339: $ RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
340: $ IWORK( IPERM ), IWORK( IGIVPT ),
341: $ IWORK( IGIVCL ), RWORK( IGIVNM ),
342: $ Q( 1, SUBMAT ), RWORK( IWREM ),
343: $ IWORK( SUBPBS+1 ), INFO )
344: IF( INFO.GT.0 ) THEN
345: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
346: RETURN
347: END IF
348: IWORK( I / 2+1 ) = IWORK( I+2 )
349: 90 CONTINUE
350: SUBPBS = SUBPBS / 2
351: CURLVL = CURLVL + 1
352: GO TO 80
353: END IF
354: *
355: * end while
356: *
357: * Re-merge the eigenvalues/vectors which were deflated at the final
358: * merge step.
359: *
360: DO 100 I = 1, N
361: J = IWORK( INDXQ+I )
362: RWORK( I ) = D( J )
363: CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
364: 100 CONTINUE
365: CALL DCOPY( N, RWORK, 1, D, 1 )
366: *
367: RETURN
368: *
369: * End of ZLAED0
370: *
371: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaed0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack