1: *> \brief \b DSTEDC
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, dimension (MAX(1,LWORK))
109: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
110: *> \endverbatim
111: *>
112: *> \param[in] LWORK
113: *> \verbatim
114: *> LWORK is INTEGER
115: *> The dimension of the array WORK.
116: *> If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
117: *> If COMPZ = 'V' and N > 1 then LWORK must be at least
118: *> ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
119: *> where lg( N ) = smallest integer k such
120: *> that 2**k >= N.
121: *> If COMPZ = 'I' and N > 1 then LWORK must be at least
122: *> ( 1 + 4*N + N**2 ).
123: *> Note that for COMPZ = 'I' or 'V', then if N is less than or
124: *> equal to the minimum divide size, usually 25, then LWORK need
125: *> only be max(1,2*(N-1)).
126: *>
127: *> If LWORK = -1, then a workspace query is assumed; the routine
128: *> only calculates the optimal size of the WORK array, returns
129: *> this value as the first entry of the WORK array, and no error
130: *> message related to LWORK is issued by XERBLA.
131: *> \endverbatim
132: *>
133: *> \param[out] IWORK
134: *> \verbatim
135: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
136: *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
137: *> \endverbatim
138: *>
139: *> \param[in] LIWORK
140: *> \verbatim
141: *> LIWORK is INTEGER
142: *> The dimension of the array IWORK.
143: *> If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
144: *> If COMPZ = 'V' and N > 1 then LIWORK must be at least
145: *> ( 6 + 6*N + 5*N*lg N ).
146: *> If COMPZ = 'I' and N > 1 then LIWORK must be at least
147: *> ( 3 + 5*N ).
148: *> Note that for COMPZ = 'I' or 'V', then if N is less than or
149: *> equal to the minimum divide size, usually 25, then LIWORK
150: *> need only be 1.
151: *>
152: *> If LIWORK = -1, then a workspace query is assumed; the
153: *> routine only calculates the optimal size of the IWORK array,
154: *> returns this value as the first entry of the IWORK array, and
155: *> no error message related to LIWORK is issued by XERBLA.
156: *> \endverbatim
157: *>
158: *> \param[out] INFO
159: *> \verbatim
160: *> INFO is INTEGER
161: *> = 0: successful exit.
162: *> < 0: if INFO = -i, the i-th argument had an illegal value.
163: *> > 0: The algorithm failed to compute an eigenvalue while
164: *> working on the submatrix lying in rows and columns
165: *> INFO/(N+1) through mod(INFO,N+1).
166: *> \endverbatim
167: *
168: * Authors:
169: * ========
170: *
171: *> \author Univ. of Tennessee
172: *> \author Univ. of California Berkeley
173: *> \author Univ. of Colorado Denver
174: *> \author NAG Ltd.
175: *
176: *> \date June 2017
177: *
178: *> \ingroup auxOTHERcomputational
179: *
180: *> \par Contributors:
181: * ==================
182: *>
183: *> Jeff Rutter, Computer Science Division, University of California
184: *> at Berkeley, USA \n
185: *> Modified by Francoise Tisseur, University of Tennessee
186: *>
187: * =====================================================================
188: SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
189: $ LIWORK, INFO )
190: *
191: * -- LAPACK computational routine (version 3.7.1) --
192: * -- LAPACK is a software package provided by Univ. of Tennessee, --
193: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194: * June 2017
195: *
196: * .. Scalar Arguments ..
197: CHARACTER COMPZ
198: INTEGER INFO, LDZ, LIWORK, LWORK, N
199: * ..
200: * .. Array Arguments ..
201: INTEGER IWORK( * )
202: DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
203: * ..
204: *
205: * =====================================================================
206: *
207: * .. Parameters ..
208: DOUBLE PRECISION ZERO, ONE, TWO
209: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
210: * ..
211: * .. Local Scalars ..
212: LOGICAL LQUERY
213: INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
214: $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
215: DOUBLE PRECISION EPS, ORGNRM, P, TINY
216: * ..
217: * .. External Functions ..
218: LOGICAL LSAME
219: INTEGER ILAENV
220: DOUBLE PRECISION DLAMCH, DLANST
221: EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
222: * ..
223: * .. External Subroutines ..
224: EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
225: $ DSTEQR, DSTERF, DSWAP, XERBLA
226: * ..
227: * .. Intrinsic Functions ..
228: INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT
229: * ..
230: * .. Executable Statements ..
231: *
232: * Test the input parameters.
233: *
234: INFO = 0
235: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
236: *
237: IF( LSAME( COMPZ, 'N' ) ) THEN
238: ICOMPZ = 0
239: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
240: ICOMPZ = 1
241: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
242: ICOMPZ = 2
243: ELSE
244: ICOMPZ = -1
245: END IF
246: IF( ICOMPZ.LT.0 ) THEN
247: INFO = -1
248: ELSE IF( N.LT.0 ) THEN
249: INFO = -2
250: ELSE IF( ( LDZ.LT.1 ) .OR.
251: $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
252: INFO = -6
253: END IF
254: *
255: IF( INFO.EQ.0 ) THEN
256: *
257: * Compute the workspace requirements
258: *
259: SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
260: IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
261: LIWMIN = 1
262: LWMIN = 1
263: ELSE IF( N.LE.SMLSIZ ) THEN
264: LIWMIN = 1
265: LWMIN = 2*( N - 1 )
266: ELSE
267: LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
268: IF( 2**LGN.LT.N )
269: $ LGN = LGN + 1
270: IF( 2**LGN.LT.N )
271: $ LGN = LGN + 1
272: IF( ICOMPZ.EQ.1 ) THEN
273: LWMIN = 1 + 3*N + 2*N*LGN + 4*N**2
274: LIWMIN = 6 + 6*N + 5*N*LGN
275: ELSE IF( ICOMPZ.EQ.2 ) THEN
276: LWMIN = 1 + 4*N + N**2
277: LIWMIN = 3 + 5*N
278: END IF
279: END IF
280: WORK( 1 ) = LWMIN
281: IWORK( 1 ) = LIWMIN
282: *
283: IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
284: INFO = -8
285: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
286: INFO = -10
287: END IF
288: END IF
289: *
290: IF( INFO.NE.0 ) THEN
291: CALL XERBLA( 'DSTEDC', -INFO )
292: RETURN
293: ELSE IF (LQUERY) THEN
294: RETURN
295: END IF
296: *
297: * Quick return if possible
298: *
299: IF( N.EQ.0 )
300: $ RETURN
301: IF( N.EQ.1 ) THEN
302: IF( ICOMPZ.NE.0 )
303: $ Z( 1, 1 ) = ONE
304: RETURN
305: END IF
306: *
307: * If the following conditional clause is removed, then the routine
308: * will use the Divide and Conquer routine to compute only the
309: * eigenvalues, which requires (3N + 3N**2) real workspace and
310: * (2 + 5N + 2N lg(N)) integer workspace.
311: * Since on many architectures DSTERF is much faster than any other
312: * algorithm for finding eigenvalues only, it is used here
313: * as the default. If the conditional clause is removed, then
314: * information on the size of workspace needs to be changed.
315: *
316: * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
317: *
318: IF( ICOMPZ.EQ.0 ) THEN
319: CALL DSTERF( N, D, E, INFO )
320: GO TO 50
321: END IF
322: *
323: * If N is smaller than the minimum divide size (SMLSIZ+1), then
324: * solve the problem with another solver.
325: *
326: IF( N.LE.SMLSIZ ) THEN
327: *
328: CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
329: *
330: ELSE
331: *
332: * If COMPZ = 'V', the Z matrix must be stored elsewhere for later
333: * use.
334: *
335: IF( ICOMPZ.EQ.1 ) THEN
336: STOREZ = 1 + N*N
337: ELSE
338: STOREZ = 1
339: END IF
340: *
341: IF( ICOMPZ.EQ.2 ) THEN
342: CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
343: END IF
344: *
345: * Scale.
346: *
347: ORGNRM = DLANST( 'M', N, D, E )
348: IF( ORGNRM.EQ.ZERO )
349: $ GO TO 50
350: *
351: EPS = DLAMCH( 'Epsilon' )
352: *
353: START = 1
354: *
355: * while ( START <= N )
356: *
357: 10 CONTINUE
358: IF( START.LE.N ) THEN
359: *
360: * Let FINISH be the position of the next subdiagonal entry
361: * such that E( FINISH ) <= TINY or FINISH = N if no such
362: * subdiagonal exists. The matrix identified by the elements
363: * between START and FINISH constitutes an independent
364: * sub-problem.
365: *
366: FINISH = START
367: 20 CONTINUE
368: IF( FINISH.LT.N ) THEN
369: TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
370: $ SQRT( ABS( D( FINISH+1 ) ) )
371: IF( ABS( E( FINISH ) ).GT.TINY ) THEN
372: FINISH = FINISH + 1
373: GO TO 20
374: END IF
375: END IF
376: *
377: * (Sub) Problem determined. Compute its size and solve it.
378: *
379: M = FINISH - START + 1
380: IF( M.EQ.1 ) THEN
381: START = FINISH + 1
382: GO TO 10
383: END IF
384: IF( M.GT.SMLSIZ ) THEN
385: *
386: * Scale.
387: *
388: ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
389: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
390: $ INFO )
391: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
392: $ M-1, INFO )
393: *
394: IF( ICOMPZ.EQ.1 ) THEN
395: STRTRW = 1
396: ELSE
397: STRTRW = START
398: END IF
399: CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
400: $ Z( STRTRW, START ), LDZ, WORK( 1 ), N,
401: $ WORK( STOREZ ), IWORK, INFO )
402: IF( INFO.NE.0 ) THEN
403: INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
404: $ MOD( INFO, ( M+1 ) ) + START - 1
405: GO TO 50
406: END IF
407: *
408: * Scale back.
409: *
410: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
411: $ INFO )
412: *
413: ELSE
414: IF( ICOMPZ.EQ.1 ) THEN
415: *
416: * Since QR won't update a Z matrix which is larger than
417: * the length of D, we must solve the sub-problem in a
418: * workspace and then multiply back into Z.
419: *
420: CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
421: $ WORK( M*M+1 ), INFO )
422: CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
423: $ WORK( STOREZ ), N )
424: CALL DGEMM( 'N', 'N', N, M, M, ONE,
425: $ WORK( STOREZ ), N, WORK, M, ZERO,
426: $ Z( 1, START ), LDZ )
427: ELSE IF( ICOMPZ.EQ.2 ) THEN
428: CALL DSTEQR( 'I', M, D( START ), E( START ),
429: $ Z( START, START ), LDZ, WORK, INFO )
430: ELSE
431: CALL DSTERF( M, D( START ), E( START ), INFO )
432: END IF
433: IF( INFO.NE.0 ) THEN
434: INFO = START*( N+1 ) + FINISH
435: GO TO 50
436: END IF
437: END IF
438: *
439: START = FINISH + 1
440: GO TO 10
441: END IF
442: *
443: * endwhile
444: *
445: IF( ICOMPZ.EQ.0 ) THEN
446: *
447: * Use Quick Sort
448: *
449: CALL DLASRT( 'I', N, D, INFO )
450: *
451: ELSE
452: *
453: * Use Selection Sort to minimize swaps of eigenvectors
454: *
455: DO 40 II = 2, N
456: I = II - 1
457: K = I
458: P = D( I )
459: DO 30 J = II, N
460: IF( D( J ).LT.P ) THEN
461: K = J
462: P = D( J )
463: END IF
464: 30 CONTINUE
465: IF( K.NE.I ) THEN
466: D( K ) = D( I )
467: D( I ) = P
468: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
469: END IF
470: 40 CONTINUE
471: END IF
472: END IF
473: *
474: 50 CONTINUE
475: WORK( 1 ) = LWMIN
476: IWORK( 1 ) = LIWMIN
477: *
478: RETURN
479: *
480: * End of DSTEDC
481: *
482: 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