LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dstedc.f
Go to the documentation of this file.
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
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*> \ingroup auxOTHERcomputational
177*
178*> \par Contributors:
179* ==================
180*>
181*> Jeff Rutter, Computer Science Division, University of California
182*> at Berkeley, USA \n
183*> Modified by Francoise Tisseur, University of Tennessee
184*>
185* =====================================================================
186 SUBROUTINE dstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
187 \$ LIWORK, INFO )
188*
189* -- LAPACK computational routine --
190* -- LAPACK is a software package provided by Univ. of Tennessee, --
191* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
192*
193* .. Scalar Arguments ..
194 CHARACTER COMPZ
195 INTEGER INFO, LDZ, LIWORK, LWORK, N
196* ..
197* .. Array Arguments ..
198 INTEGER IWORK( * )
199 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
200* ..
201*
202* =====================================================================
203*
204* .. Parameters ..
205 DOUBLE PRECISION ZERO, ONE, TWO
206 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
207* ..
208* .. Local Scalars ..
209 LOGICAL LQUERY
210 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
211 \$ lwmin, m, smlsiz, start, storez, strtrw
212 DOUBLE PRECISION EPS, ORGNRM, P, TINY
213* ..
214* .. External Functions ..
215 LOGICAL LSAME
216 INTEGER ILAENV
217 DOUBLE PRECISION DLAMCH, DLANST
218 EXTERNAL lsame, ilaenv, dlamch, dlanst
219* ..
220* .. External Subroutines ..
221 EXTERNAL dgemm, dlacpy, dlaed0, dlascl, dlaset, dlasrt,
223* ..
224* .. Intrinsic Functions ..
225 INTRINSIC abs, dble, int, log, max, mod, sqrt
226* ..
227* .. Executable Statements ..
228*
229* Test the input parameters.
230*
231 info = 0
232 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
233*
234 IF( lsame( compz, 'N' ) ) THEN
235 icompz = 0
236 ELSE IF( lsame( compz, 'V' ) ) THEN
237 icompz = 1
238 ELSE IF( lsame( compz, 'I' ) ) THEN
239 icompz = 2
240 ELSE
241 icompz = -1
242 END IF
243 IF( icompz.LT.0 ) THEN
244 info = -1
245 ELSE IF( n.LT.0 ) THEN
246 info = -2
247 ELSE IF( ( ldz.LT.1 ) .OR.
248 \$ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
249 info = -6
250 END IF
251*
252 IF( info.EQ.0 ) THEN
253*
254* Compute the workspace requirements
255*
256 smlsiz = ilaenv( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
257 IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
258 liwmin = 1
259 lwmin = 1
260 ELSE IF( n.LE.smlsiz ) THEN
261 liwmin = 1
262 lwmin = 2*( n - 1 )
263 ELSE
264 lgn = int( log( dble( n ) )/log( two ) )
265 IF( 2**lgn.LT.n )
266 \$ lgn = lgn + 1
267 IF( 2**lgn.LT.n )
268 \$ lgn = lgn + 1
269 IF( icompz.EQ.1 ) THEN
270 lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
271 liwmin = 6 + 6*n + 5*n*lgn
272 ELSE IF( icompz.EQ.2 ) THEN
273 lwmin = 1 + 4*n + n**2
274 liwmin = 3 + 5*n
275 END IF
276 END IF
277 work( 1 ) = lwmin
278 iwork( 1 ) = liwmin
279*
280 IF( lwork.LT.lwmin .AND. .NOT. lquery ) THEN
281 info = -8
282 ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery ) THEN
283 info = -10
284 END IF
285 END IF
286*
287 IF( info.NE.0 ) THEN
288 CALL xerbla( 'DSTEDC', -info )
289 RETURN
290 ELSE IF (lquery) THEN
291 RETURN
292 END IF
293*
294* Quick return if possible
295*
296 IF( n.EQ.0 )
297 \$ RETURN
298 IF( n.EQ.1 ) THEN
299 IF( icompz.NE.0 )
300 \$ z( 1, 1 ) = one
301 RETURN
302 END IF
303*
304* If the following conditional clause is removed, then the routine
305* will use the Divide and Conquer routine to compute only the
306* eigenvalues, which requires (3N + 3N**2) real workspace and
307* (2 + 5N + 2N lg(N)) integer workspace.
308* Since on many architectures DSTERF is much faster than any other
309* algorithm for finding eigenvalues only, it is used here
310* as the default. If the conditional clause is removed, then
311* information on the size of workspace needs to be changed.
312*
313* If COMPZ = 'N', use DSTERF to compute the eigenvalues.
314*
315 IF( icompz.EQ.0 ) THEN
316 CALL dsterf( n, d, e, info )
317 GO TO 50
318 END IF
319*
320* If N is smaller than the minimum divide size (SMLSIZ+1), then
321* solve the problem with another solver.
322*
323 IF( n.LE.smlsiz ) THEN
324*
325 CALL dsteqr( compz, n, d, e, z, ldz, work, info )
326*
327 ELSE
328*
329* If COMPZ = 'V', the Z matrix must be stored elsewhere for later
330* use.
331*
332 IF( icompz.EQ.1 ) THEN
333 storez = 1 + n*n
334 ELSE
335 storez = 1
336 END IF
337*
338 IF( icompz.EQ.2 ) THEN
339 CALL dlaset( 'Full', n, n, zero, one, z, ldz )
340 END IF
341*
342* Scale.
343*
344 orgnrm = dlanst( 'M', n, d, e )
345 IF( orgnrm.EQ.zero )
346 \$ GO TO 50
347*
348 eps = dlamch( 'Epsilon' )
349*
350 start = 1
351*
352* while ( START <= N )
353*
354 10 CONTINUE
355 IF( start.LE.n ) THEN
356*
357* Let FINISH be the position of the next subdiagonal entry
358* such that E( FINISH ) <= TINY or FINISH = N if no such
359* subdiagonal exists. The matrix identified by the elements
360* between START and FINISH constitutes an independent
361* sub-problem.
362*
363 finish = start
364 20 CONTINUE
365 IF( finish.LT.n ) THEN
366 tiny = eps*sqrt( abs( d( finish ) ) )*
367 \$ sqrt( abs( d( finish+1 ) ) )
368 IF( abs( e( finish ) ).GT.tiny ) THEN
369 finish = finish + 1
370 GO TO 20
371 END IF
372 END IF
373*
374* (Sub) Problem determined. Compute its size and solve it.
375*
376 m = finish - start + 1
377 IF( m.EQ.1 ) THEN
378 start = finish + 1
379 GO TO 10
380 END IF
381 IF( m.GT.smlsiz ) THEN
382*
383* Scale.
384*
385 orgnrm = dlanst( 'M', m, d( start ), e( start ) )
386 CALL dlascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
387 \$ info )
388 CALL dlascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
389 \$ m-1, info )
390*
391 IF( icompz.EQ.1 ) THEN
392 strtrw = 1
393 ELSE
394 strtrw = start
395 END IF
396 CALL dlaed0( icompz, n, m, d( start ), e( start ),
397 \$ z( strtrw, start ), ldz, work( 1 ), n,
398 \$ work( storez ), iwork, info )
399 IF( info.NE.0 ) THEN
400 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
401 \$ mod( info, ( m+1 ) ) + start - 1
402 GO TO 50
403 END IF
404*
405* Scale back.
406*
407 CALL dlascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
408 \$ info )
409*
410 ELSE
411 IF( icompz.EQ.1 ) THEN
412*
413* Since QR won't update a Z matrix which is larger than
414* the length of D, we must solve the sub-problem in a
415* workspace and then multiply back into Z.
416*
417 CALL dsteqr( 'I', m, d( start ), e( start ), work, m,
418 \$ work( m*m+1 ), info )
419 CALL dlacpy( 'A', n, m, z( 1, start ), ldz,
420 \$ work( storez ), n )
421 CALL dgemm( 'N', 'N', n, m, m, one,
422 \$ work( storez ), n, work, m, zero,
423 \$ z( 1, start ), ldz )
424 ELSE IF( icompz.EQ.2 ) THEN
425 CALL dsteqr( 'I', m, d( start ), e( start ),
426 \$ z( start, start ), ldz, work, info )
427 ELSE
428 CALL dsterf( m, d( start ), e( start ), info )
429 END IF
430 IF( info.NE.0 ) THEN
431 info = start*( n+1 ) + finish
432 GO TO 50
433 END IF
434 END IF
435*
436 start = finish + 1
437 GO TO 10
438 END IF
439*
440* endwhile
441*
442 IF( icompz.EQ.0 ) THEN
443*
444* Use Quick Sort
445*
446 CALL dlasrt( 'I', n, d, info )
447*
448 ELSE
449*
450* Use Selection Sort to minimize swaps of eigenvectors
451*
452 DO 40 ii = 2, n
453 i = ii - 1
454 k = i
455 p = d( i )
456 DO 30 j = ii, n
457 IF( d( j ).LT.p ) THEN
458 k = j
459 p = d( j )
460 END IF
461 30 CONTINUE
462 IF( k.NE.i ) THEN
463 d( k ) = d( i )
464 d( i ) = p
465 CALL dswap( n, z( 1, i ), 1, z( 1, k ), 1 )
466 END IF
467 40 CONTINUE
468 END IF
469 END IF
470*
471 50 CONTINUE
472 work( 1 ) = lwmin
473 iwork( 1 ) = liwmin
474*
475 RETURN
476*
477* End of DSTEDC
478*
479 END
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:188
subroutine dlaed0(ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)
DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: dlaed0.f:172
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187