Annotation of rpl/lapack/lapack/dstedc.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b DSTEBZ
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DSTEDC + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstedc.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstedc.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstedc.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
! 22: * LIWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER COMPZ
! 26: * INTEGER INFO, LDZ, LIWORK, LWORK, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IWORK( * )
! 30: * DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
! 40: *> symmetric tridiagonal matrix using the divide and conquer method.
! 41: *> The eigenvectors of a full or band real symmetric matrix can also be
! 42: *> found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
! 43: *> matrix to tridiagonal form.
! 44: *>
! 45: *> This code makes very mild assumptions about floating point
! 46: *> arithmetic. It will work on machines with a guard digit in
! 47: *> add/subtract, or on those binary machines without guard digits
! 48: *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
! 49: *> It could conceivably fail on hexadecimal or decimal machines
! 50: *> without guard digits, but we know of none. See DLAED3 for details.
! 51: *> \endverbatim
! 52: *
! 53: * Arguments:
! 54: * ==========
! 55: *
! 56: *> \param[in] COMPZ
! 57: *> \verbatim
! 58: *> COMPZ is CHARACTER*1
! 59: *> = 'N': Compute eigenvalues only.
! 60: *> = 'I': Compute eigenvectors of tridiagonal matrix also.
! 61: *> = 'V': Compute eigenvectors of original dense symmetric
! 62: *> matrix also. On entry, Z contains the orthogonal
! 63: *> matrix used to reduce the original matrix to
! 64: *> tridiagonal form.
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in] N
! 68: *> \verbatim
! 69: *> N is INTEGER
! 70: *> The dimension of the symmetric tridiagonal matrix. N >= 0.
! 71: *> \endverbatim
! 72: *>
! 73: *> \param[in,out] D
! 74: *> \verbatim
! 75: *> D is DOUBLE PRECISION array, dimension (N)
! 76: *> On entry, the diagonal elements of the tridiagonal matrix.
! 77: *> On exit, if INFO = 0, the eigenvalues in ascending order.
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in,out] E
! 81: *> \verbatim
! 82: *> E is DOUBLE PRECISION array, dimension (N-1)
! 83: *> On entry, the subdiagonal elements of the tridiagonal matrix.
! 84: *> On exit, E has been destroyed.
! 85: *> \endverbatim
! 86: *>
! 87: *> \param[in,out] Z
! 88: *> \verbatim
! 89: *> Z is DOUBLE PRECISION array, dimension (LDZ,N)
! 90: *> On entry, if COMPZ = 'V', then Z contains the orthogonal
! 91: *> matrix used in the reduction to tridiagonal form.
! 92: *> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
! 93: *> orthonormal eigenvectors of the original symmetric matrix,
! 94: *> and if COMPZ = 'I', Z contains the orthonormal eigenvectors
! 95: *> of the symmetric tridiagonal matrix.
! 96: *> If COMPZ = 'N', then Z is not referenced.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] LDZ
! 100: *> \verbatim
! 101: *> LDZ is INTEGER
! 102: *> The leading dimension of the array Z. LDZ >= 1.
! 103: *> If eigenvectors are desired, then LDZ >= max(1,N).
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[out] WORK
! 107: *> \verbatim
! 108: *> WORK is DOUBLE PRECISION array,
! 109: *> dimension (LWORK)
! 110: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[in] LWORK
! 114: *> \verbatim
! 115: *> LWORK is INTEGER
! 116: *> The dimension of the array WORK.
! 117: *> If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
! 118: *> If COMPZ = 'V' and N > 1 then LWORK must be at least
! 119: *> ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
! 120: *> where lg( N ) = smallest integer k such
! 121: *> that 2**k >= N.
! 122: *> If COMPZ = 'I' and N > 1 then LWORK must be at least
! 123: *> ( 1 + 4*N + N**2 ).
! 124: *> Note that for COMPZ = 'I' or 'V', then if N is less than or
! 125: *> equal to the minimum divide size, usually 25, then LWORK need
! 126: *> only be max(1,2*(N-1)).
! 127: *>
! 128: *> If LWORK = -1, then a workspace query is assumed; the routine
! 129: *> only calculates the optimal size of the WORK array, returns
! 130: *> this value as the first entry of the WORK array, and no error
! 131: *> message related to LWORK is issued by XERBLA.
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[out] IWORK
! 135: *> \verbatim
! 136: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
! 137: *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
! 138: *> \endverbatim
! 139: *>
! 140: *> \param[in] LIWORK
! 141: *> \verbatim
! 142: *> LIWORK is INTEGER
! 143: *> The dimension of the array IWORK.
! 144: *> If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
! 145: *> If COMPZ = 'V' and N > 1 then LIWORK must be at least
! 146: *> ( 6 + 6*N + 5*N*lg N ).
! 147: *> If COMPZ = 'I' and N > 1 then LIWORK must be at least
! 148: *> ( 3 + 5*N ).
! 149: *> Note that for COMPZ = 'I' or 'V', then if N is less than or
! 150: *> equal to the minimum divide size, usually 25, then LIWORK
! 151: *> need only be 1.
! 152: *>
! 153: *> If LIWORK = -1, then a workspace query is assumed; the
! 154: *> routine only calculates the optimal size of the IWORK array,
! 155: *> returns this value as the first entry of the IWORK array, and
! 156: *> no error message related to LIWORK is issued by XERBLA.
! 157: *> \endverbatim
! 158: *>
! 159: *> \param[out] INFO
! 160: *> \verbatim
! 161: *> INFO is INTEGER
! 162: *> = 0: successful exit.
! 163: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 164: *> > 0: The algorithm failed to compute an eigenvalue while
! 165: *> working on the submatrix lying in rows and columns
! 166: *> INFO/(N+1) through mod(INFO,N+1).
! 167: *> \endverbatim
! 168: *
! 169: * Authors:
! 170: * ========
! 171: *
! 172: *> \author Univ. of Tennessee
! 173: *> \author Univ. of California Berkeley
! 174: *> \author Univ. of Colorado Denver
! 175: *> \author NAG Ltd.
! 176: *
! 177: *> \date November 2011
! 178: *
! 179: *> \ingroup auxOTHERcomputational
! 180: *
! 181: *> \par Contributors:
! 182: * ==================
! 183: *>
! 184: *> Jeff Rutter, Computer Science Division, University of California
! 185: *> at Berkeley, USA \n
! 186: *> Modified by Francoise Tisseur, University of Tennessee
! 187: *>
! 188: * =====================================================================
1.1 bertrand 189: SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
190: $ LIWORK, INFO )
191: *
1.8 ! bertrand 192: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 193: * -- LAPACK is a software package provided by Univ. of Tennessee, --
194: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 195: * November 2011
1.1 bertrand 196: *
197: * .. Scalar Arguments ..
198: CHARACTER COMPZ
199: INTEGER INFO, LDZ, LIWORK, LWORK, N
200: * ..
201: * .. Array Arguments ..
202: INTEGER IWORK( * )
203: DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
204: * ..
205: *
206: * =====================================================================
207: *
208: * .. Parameters ..
209: DOUBLE PRECISION ZERO, ONE, TWO
210: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
211: * ..
212: * .. Local Scalars ..
213: LOGICAL LQUERY
214: INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
215: $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
216: DOUBLE PRECISION EPS, ORGNRM, P, TINY
217: * ..
218: * .. External Functions ..
219: LOGICAL LSAME
220: INTEGER ILAENV
221: DOUBLE PRECISION DLAMCH, DLANST
222: EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
223: * ..
224: * .. External Subroutines ..
225: EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
226: $ DSTEQR, DSTERF, DSWAP, XERBLA
227: * ..
228: * .. Intrinsic Functions ..
229: INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT
230: * ..
231: * .. Executable Statements ..
232: *
233: * Test the input parameters.
234: *
235: INFO = 0
236: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
237: *
238: IF( LSAME( COMPZ, 'N' ) ) THEN
239: ICOMPZ = 0
240: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
241: ICOMPZ = 1
242: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
243: ICOMPZ = 2
244: ELSE
245: ICOMPZ = -1
246: END IF
247: IF( ICOMPZ.LT.0 ) THEN
248: INFO = -1
249: ELSE IF( N.LT.0 ) THEN
250: INFO = -2
251: ELSE IF( ( LDZ.LT.1 ) .OR.
252: $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
253: INFO = -6
254: END IF
255: *
256: IF( INFO.EQ.0 ) THEN
257: *
258: * Compute the workspace requirements
259: *
260: SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
261: IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
262: LIWMIN = 1
263: LWMIN = 1
264: ELSE IF( N.LE.SMLSIZ ) THEN
265: LIWMIN = 1
266: LWMIN = 2*( N - 1 )
267: ELSE
268: LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
269: IF( 2**LGN.LT.N )
270: $ LGN = LGN + 1
271: IF( 2**LGN.LT.N )
272: $ LGN = LGN + 1
273: IF( ICOMPZ.EQ.1 ) THEN
1.8 ! bertrand 274: LWMIN = 1 + 3*N + 2*N*LGN + 4*N**2
1.1 bertrand 275: LIWMIN = 6 + 6*N + 5*N*LGN
276: ELSE IF( ICOMPZ.EQ.2 ) THEN
277: LWMIN = 1 + 4*N + N**2
278: LIWMIN = 3 + 5*N
279: END IF
280: END IF
281: WORK( 1 ) = LWMIN
282: IWORK( 1 ) = LIWMIN
283: *
284: IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
285: INFO = -8
286: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
287: INFO = -10
288: END IF
289: END IF
290: *
291: IF( INFO.NE.0 ) THEN
292: CALL XERBLA( 'DSTEDC', -INFO )
293: RETURN
294: ELSE IF (LQUERY) THEN
295: RETURN
296: END IF
297: *
298: * Quick return if possible
299: *
300: IF( N.EQ.0 )
301: $ RETURN
302: IF( N.EQ.1 ) THEN
303: IF( ICOMPZ.NE.0 )
304: $ Z( 1, 1 ) = ONE
305: RETURN
306: END IF
307: *
308: * If the following conditional clause is removed, then the routine
309: * will use the Divide and Conquer routine to compute only the
310: * eigenvalues, which requires (3N + 3N**2) real workspace and
311: * (2 + 5N + 2N lg(N)) integer workspace.
312: * Since on many architectures DSTERF is much faster than any other
313: * algorithm for finding eigenvalues only, it is used here
314: * as the default. If the conditional clause is removed, then
315: * information on the size of workspace needs to be changed.
316: *
317: * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
318: *
319: IF( ICOMPZ.EQ.0 ) THEN
320: CALL DSTERF( N, D, E, INFO )
321: GO TO 50
322: END IF
323: *
324: * If N is smaller than the minimum divide size (SMLSIZ+1), then
325: * solve the problem with another solver.
326: *
327: IF( N.LE.SMLSIZ ) THEN
328: *
329: CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
330: *
331: ELSE
332: *
333: * If COMPZ = 'V', the Z matrix must be stored elsewhere for later
334: * use.
335: *
336: IF( ICOMPZ.EQ.1 ) THEN
337: STOREZ = 1 + N*N
338: ELSE
339: STOREZ = 1
340: END IF
341: *
342: IF( ICOMPZ.EQ.2 ) THEN
343: CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
344: END IF
345: *
346: * Scale.
347: *
348: ORGNRM = DLANST( 'M', N, D, E )
349: IF( ORGNRM.EQ.ZERO )
350: $ GO TO 50
351: *
352: EPS = DLAMCH( 'Epsilon' )
353: *
354: START = 1
355: *
356: * while ( START <= N )
357: *
358: 10 CONTINUE
359: IF( START.LE.N ) THEN
360: *
361: * Let FINISH be the position of the next subdiagonal entry
362: * such that E( FINISH ) <= TINY or FINISH = N if no such
363: * subdiagonal exists. The matrix identified by the elements
364: * between START and FINISH constitutes an independent
365: * sub-problem.
366: *
367: FINISH = START
368: 20 CONTINUE
369: IF( FINISH.LT.N ) THEN
370: TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
371: $ SQRT( ABS( D( FINISH+1 ) ) )
372: IF( ABS( E( FINISH ) ).GT.TINY ) THEN
373: FINISH = FINISH + 1
374: GO TO 20
375: END IF
376: END IF
377: *
378: * (Sub) Problem determined. Compute its size and solve it.
379: *
380: M = FINISH - START + 1
381: IF( M.EQ.1 ) THEN
382: START = FINISH + 1
383: GO TO 10
384: END IF
385: IF( M.GT.SMLSIZ ) THEN
386: *
387: * Scale.
388: *
389: ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
390: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
391: $ INFO )
392: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
393: $ M-1, INFO )
394: *
395: IF( ICOMPZ.EQ.1 ) THEN
396: STRTRW = 1
397: ELSE
398: STRTRW = START
399: END IF
400: CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
401: $ Z( STRTRW, START ), LDZ, WORK( 1 ), N,
402: $ WORK( STOREZ ), IWORK, INFO )
403: IF( INFO.NE.0 ) THEN
404: INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
405: $ MOD( INFO, ( M+1 ) ) + START - 1
406: GO TO 50
407: END IF
408: *
409: * Scale back.
410: *
411: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
412: $ INFO )
413: *
414: ELSE
415: IF( ICOMPZ.EQ.1 ) THEN
416: *
417: * Since QR won't update a Z matrix which is larger than
418: * the length of D, we must solve the sub-problem in a
419: * workspace and then multiply back into Z.
420: *
421: CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
422: $ WORK( M*M+1 ), INFO )
423: CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
424: $ WORK( STOREZ ), N )
425: CALL DGEMM( 'N', 'N', N, M, M, ONE,
426: $ WORK( STOREZ ), N, WORK, M, ZERO,
427: $ Z( 1, START ), LDZ )
428: ELSE IF( ICOMPZ.EQ.2 ) THEN
429: CALL DSTEQR( 'I', M, D( START ), E( START ),
430: $ Z( START, START ), LDZ, WORK, INFO )
431: ELSE
432: CALL DSTERF( M, D( START ), E( START ), INFO )
433: END IF
434: IF( INFO.NE.0 ) THEN
435: INFO = START*( N+1 ) + FINISH
436: GO TO 50
437: END IF
438: END IF
439: *
440: START = FINISH + 1
441: GO TO 10
442: END IF
443: *
444: * endwhile
445: *
446: * If the problem split any number of times, then the eigenvalues
447: * will not be properly ordered. Here we permute the eigenvalues
448: * (and the associated eigenvectors) into ascending order.
449: *
450: IF( M.NE.N ) THEN
451: IF( ICOMPZ.EQ.0 ) THEN
452: *
453: * Use Quick Sort
454: *
455: CALL DLASRT( 'I', N, D, INFO )
456: *
457: ELSE
458: *
459: * Use Selection Sort to minimize swaps of eigenvectors
460: *
461: DO 40 II = 2, N
462: I = II - 1
463: K = I
464: P = D( I )
465: DO 30 J = II, N
466: IF( D( J ).LT.P ) THEN
467: K = J
468: P = D( J )
469: END IF
470: 30 CONTINUE
471: IF( K.NE.I ) THEN
472: D( K ) = D( I )
473: D( I ) = P
474: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
475: END IF
476: 40 CONTINUE
477: END IF
478: END IF
479: END IF
480: *
481: 50 CONTINUE
482: WORK( 1 ) = LWMIN
483: IWORK( 1 ) = LIWMIN
484: *
485: RETURN
486: *
487: * End of DSTEDC
488: *
489: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dstedc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack