LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sstedc.f
Go to the documentation of this file.
1*> \brief \b SSTEDC
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SSTEDC + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstedc.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstedc.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstedc.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
20* LIWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER COMPZ
24* INTEGER INFO, LDZ, LIWORK, LWORK, N
25* ..
26* .. Array Arguments ..
27* INTEGER IWORK( * )
28* REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SSTEDC computes all eigenvalues and, optionally, eigenvectors of a
38*> symmetric tridiagonal matrix using the divide and conquer method.
39*> The eigenvectors of a full or band real symmetric matrix can also be
40*> found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this
41*> matrix to tridiagonal form.
42*>
43*> \endverbatim
44*
45* Arguments:
46* ==========
47*
48*> \param[in] COMPZ
49*> \verbatim
50*> COMPZ is CHARACTER*1
51*> = 'N': Compute eigenvalues only.
52*> = 'I': Compute eigenvectors of tridiagonal matrix also.
53*> = 'V': Compute eigenvectors of original dense symmetric
54*> matrix also. On entry, Z contains the orthogonal
55*> matrix used to reduce the original matrix to
56*> tridiagonal form.
57*> \endverbatim
58*>
59*> \param[in] N
60*> \verbatim
61*> N is INTEGER
62*> The dimension of the symmetric tridiagonal matrix. N >= 0.
63*> \endverbatim
64*>
65*> \param[in,out] D
66*> \verbatim
67*> D is REAL array, dimension (N)
68*> On entry, the diagonal elements of the tridiagonal matrix.
69*> On exit, if INFO = 0, the eigenvalues in ascending order.
70*> \endverbatim
71*>
72*> \param[in,out] E
73*> \verbatim
74*> E is REAL array, dimension (N-1)
75*> On entry, the subdiagonal elements of the tridiagonal matrix.
76*> On exit, E has been destroyed.
77*> \endverbatim
78*>
79*> \param[in,out] Z
80*> \verbatim
81*> Z is REAL array, dimension (LDZ,N)
82*> On entry, if COMPZ = 'V', then Z contains the orthogonal
83*> matrix used in the reduction to tridiagonal form.
84*> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
85*> orthonormal eigenvectors of the original symmetric matrix,
86*> and if COMPZ = 'I', Z contains the orthonormal eigenvectors
87*> of the symmetric tridiagonal matrix.
88*> If COMPZ = 'N', then Z is not referenced.
89*> \endverbatim
90*>
91*> \param[in] LDZ
92*> \verbatim
93*> LDZ is INTEGER
94*> The leading dimension of the array Z. LDZ >= 1.
95*> If eigenvectors are desired, then LDZ >= max(1,N).
96*> \endverbatim
97*>
98*> \param[out] WORK
99*> \verbatim
100*> WORK is REAL array, dimension (MAX(1,LWORK))
101*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
102*> \endverbatim
103*>
104*> \param[in] LWORK
105*> \verbatim
106*> LWORK is INTEGER
107*> The dimension of the array WORK.
108*> If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
109*> If COMPZ = 'V' and N > 1 then LWORK must be at least
110*> ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
111*> where lg( N ) = smallest integer k such
112*> that 2**k >= N.
113*> If COMPZ = 'I' and N > 1 then LWORK must be at least
114*> ( 1 + 4*N + N**2 ).
115*> Note that for COMPZ = 'I' or 'V', then if N is less than or
116*> equal to the minimum divide size, usually 25, then LWORK need
117*> only be max(1,2*(N-1)).
118*>
119*> If LWORK = -1, then a workspace query is assumed; the routine
120*> only calculates the optimal size of the WORK array, returns
121*> this value as the first entry of the WORK array, and no error
122*> message related to LWORK is issued by XERBLA.
123*> \endverbatim
124*>
125*> \param[out] IWORK
126*> \verbatim
127*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
128*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
129*> \endverbatim
130*>
131*> \param[in] LIWORK
132*> \verbatim
133*> LIWORK is INTEGER
134*> The dimension of the array IWORK.
135*> If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
136*> If COMPZ = 'V' and N > 1 then LIWORK must be at least
137*> ( 6 + 6*N + 5*N*lg N ).
138*> If COMPZ = 'I' and N > 1 then LIWORK must be at least
139*> ( 3 + 5*N ).
140*> Note that for COMPZ = 'I' or 'V', then if N is less than or
141*> equal to the minimum divide size, usually 25, then LIWORK
142*> need only be 1.
143*>
144*> If LIWORK = -1, then a workspace query is assumed; the
145*> routine only calculates the optimal size of the IWORK array,
146*> returns this value as the first entry of the IWORK array, and
147*> no error message related to LIWORK is issued by XERBLA.
148*> \endverbatim
149*>
150*> \param[out] INFO
151*> \verbatim
152*> INFO is INTEGER
153*> = 0: successful exit.
154*> < 0: if INFO = -i, the i-th argument had an illegal value.
155*> > 0: The algorithm failed to compute an eigenvalue while
156*> working on the submatrix lying in rows and columns
157*> INFO/(N+1) through mod(INFO,N+1).
158*> \endverbatim
159*
160* Authors:
161* ========
162*
163*> \author Univ. of Tennessee
164*> \author Univ. of California Berkeley
165*> \author Univ. of Colorado Denver
166*> \author NAG Ltd.
167*
168*> \ingroup stedc
169*
170*> \par Contributors:
171* ==================
172*>
173*> Jeff Rutter, Computer Science Division, University of California
174*> at Berkeley, USA \n
175*> Modified by Francoise Tisseur, University of Tennessee
176*>
177* =====================================================================
178 SUBROUTINE sstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
179 $ LIWORK, INFO )
180*
181* -- LAPACK computational routine --
182* -- LAPACK is a software package provided by Univ. of Tennessee, --
183* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184*
185* .. Scalar Arguments ..
186 CHARACTER COMPZ
187 INTEGER INFO, LDZ, LIWORK, LWORK, N
188* ..
189* .. Array Arguments ..
190 INTEGER IWORK( * )
191 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
192* ..
193*
194* =====================================================================
195*
196* .. Parameters ..
197 REAL ZERO, ONE, TWO
198 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
199* ..
200* .. Local Scalars ..
201 LOGICAL LQUERY
202 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
203 $ lwmin, m, smlsiz, start, storez, strtrw
204 REAL EPS, ORGNRM, P, TINY
205* ..
206* .. External Functions ..
207 LOGICAL LSAME
208 INTEGER ILAENV
209 REAL SLAMCH, SLANST, SROUNDUP_LWORK
210 EXTERNAL ilaenv, lsame, slamch, slanst,
211 $ sroundup_lwork
212* ..
213* .. External Subroutines ..
214 EXTERNAL sgemm, slacpy, slaed0, slascl, slaset,
215 $ slasrt,
217* ..
218* .. Intrinsic Functions ..
219 INTRINSIC abs, int, log, max, mod, real, sqrt
220* ..
221* .. Executable Statements ..
222*
223* Test the input parameters.
224*
225 info = 0
226 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
227*
228 IF( lsame( compz, 'N' ) ) THEN
229 icompz = 0
230 ELSE IF( lsame( compz, 'V' ) ) THEN
231 icompz = 1
232 ELSE IF( lsame( compz, 'I' ) ) THEN
233 icompz = 2
234 ELSE
235 icompz = -1
236 END IF
237 IF( icompz.LT.0 ) THEN
238 info = -1
239 ELSE IF( n.LT.0 ) THEN
240 info = -2
241 ELSE IF( ( ldz.LT.1 ) .OR.
242 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
243 info = -6
244 END IF
245*
246 IF( info.EQ.0 ) THEN
247*
248* Compute the workspace requirements
249*
250 smlsiz = ilaenv( 9, 'SSTEDC', ' ', 0, 0, 0, 0 )
251 IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
252 liwmin = 1
253 lwmin = 1
254 ELSE IF( n.LE.smlsiz ) THEN
255 liwmin = 1
256 lwmin = 2*( n - 1 )
257 ELSE
258 lgn = int( log( real( n ) )/log( two ) )
259 IF( 2**lgn.LT.n )
260 $ lgn = lgn + 1
261 IF( 2**lgn.LT.n )
262 $ lgn = lgn + 1
263 IF( icompz.EQ.1 ) THEN
264 lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
265 liwmin = 6 + 6*n + 5*n*lgn
266 ELSE IF( icompz.EQ.2 ) THEN
267 lwmin = 1 + 4*n + n**2
268 liwmin = 3 + 5*n
269 END IF
270 END IF
271 work( 1 ) = sroundup_lwork(lwmin)
272 iwork( 1 ) = liwmin
273*
274 IF( lwork.LT.lwmin .AND. .NOT. lquery ) THEN
275 info = -8
276 ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery ) THEN
277 info = -10
278 END IF
279 END IF
280*
281 IF( info.NE.0 ) THEN
282 CALL xerbla( 'SSTEDC', -info )
283 RETURN
284 ELSE IF (lquery) THEN
285 RETURN
286 END IF
287*
288* Quick return if possible
289*
290 IF( n.EQ.0 )
291 $ RETURN
292 IF( n.EQ.1 ) THEN
293 IF( icompz.NE.0 )
294 $ z( 1, 1 ) = one
295 RETURN
296 END IF
297*
298* If the following conditional clause is removed, then the routine
299* will use the Divide and Conquer routine to compute only the
300* eigenvalues, which requires (3N + 3N**2) real workspace and
301* (2 + 5N + 2N lg(N)) integer workspace.
302* Since on many architectures SSTERF is much faster than any other
303* algorithm for finding eigenvalues only, it is used here
304* as the default. If the conditional clause is removed, then
305* information on the size of workspace needs to be changed.
306*
307* If COMPZ = 'N', use SSTERF to compute the eigenvalues.
308*
309 IF( icompz.EQ.0 ) THEN
310 CALL ssterf( n, d, e, info )
311 GO TO 50
312 END IF
313*
314* If N is smaller than the minimum divide size (SMLSIZ+1), then
315* solve the problem with another solver.
316*
317 IF( n.LE.smlsiz ) THEN
318*
319 CALL ssteqr( compz, n, d, e, z, ldz, work, info )
320*
321 ELSE
322*
323* If COMPZ = 'V', the Z matrix must be stored elsewhere for later
324* use.
325*
326 IF( icompz.EQ.1 ) THEN
327 storez = 1 + n*n
328 ELSE
329 storez = 1
330 END IF
331*
332 IF( icompz.EQ.2 ) THEN
333 CALL slaset( 'Full', n, n, zero, one, z, ldz )
334 END IF
335*
336* Scale.
337*
338 orgnrm = slanst( 'M', n, d, e )
339 IF( orgnrm.EQ.zero )
340 $ GO TO 50
341*
342 eps = slamch( 'Epsilon' )
343*
344 start = 1
345*
346* while ( START <= N )
347*
348 10 CONTINUE
349 IF( start.LE.n ) THEN
350*
351* Let FINISH be the position of the next subdiagonal entry
352* such that E( FINISH ) <= TINY or FINISH = N if no such
353* subdiagonal exists. The matrix identified by the elements
354* between START and FINISH constitutes an independent
355* sub-problem.
356*
357 finish = start
358 20 CONTINUE
359 IF( finish.LT.n ) THEN
360 tiny = eps*sqrt( abs( d( finish ) ) )*
361 $ sqrt( abs( d( finish+1 ) ) )
362 IF( abs( e( finish ) ).GT.tiny ) THEN
363 finish = finish + 1
364 GO TO 20
365 END IF
366 END IF
367*
368* (Sub) Problem determined. Compute its size and solve it.
369*
370 m = finish - start + 1
371 IF( m.EQ.1 ) THEN
372 start = finish + 1
373 GO TO 10
374 END IF
375 IF( m.GT.smlsiz ) THEN
376*
377* Scale.
378*
379 orgnrm = slanst( 'M', m, d( start ), e( start ) )
380 CALL slascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ),
381 $ m,
382 $ info )
383 CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1,
384 $ e( start ),
385 $ m-1, info )
386*
387 IF( icompz.EQ.1 ) THEN
388 strtrw = 1
389 ELSE
390 strtrw = start
391 END IF
392 CALL slaed0( icompz, n, m, d( start ), e( start ),
393 $ z( strtrw, start ), ldz, work( 1 ), n,
394 $ work( storez ), iwork, info )
395 IF( info.NE.0 ) THEN
396 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
397 $ mod( info, ( m+1 ) ) + start - 1
398 GO TO 50
399 END IF
400*
401* Scale back.
402*
403 CALL slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ),
404 $ m,
405 $ info )
406*
407 ELSE
408 IF( icompz.EQ.1 ) THEN
409*
410* Since QR won't update a Z matrix which is larger than
411* the length of D, we must solve the sub-problem in a
412* workspace and then multiply back into Z.
413*
414 CALL ssteqr( 'I', m, d( start ), e( start ), work,
415 $ m,
416 $ work( m*m+1 ), info )
417 CALL slacpy( 'A', n, m, z( 1, start ), ldz,
418 $ work( storez ), n )
419 CALL sgemm( 'N', 'N', n, m, m, one,
420 $ work( storez ), n, work, m, zero,
421 $ z( 1, start ), ldz )
422 ELSE IF( icompz.EQ.2 ) THEN
423 CALL ssteqr( 'I', m, d( start ), e( start ),
424 $ z( start, start ), ldz, work, info )
425 ELSE
426 CALL ssterf( m, d( start ), e( start ), info )
427 END IF
428 IF( info.NE.0 ) THEN
429 info = start*( n+1 ) + finish
430 GO TO 50
431 END IF
432 END IF
433*
434 start = finish + 1
435 GO TO 10
436 END IF
437*
438* endwhile
439*
440 IF( icompz.EQ.0 ) THEN
441*
442* Use Quick Sort
443*
444 CALL slasrt( 'I', n, d, info )
445*
446 ELSE
447*
448* Use Selection Sort to minimize swaps of eigenvectors
449*
450 DO 40 ii = 2, n
451 i = ii - 1
452 k = i
453 p = d( i )
454 DO 30 j = ii, n
455 IF( d( j ).LT.p ) THEN
456 k = j
457 p = d( j )
458 END IF
459 30 CONTINUE
460 IF( k.NE.i ) THEN
461 d( k ) = d( i )
462 d( i ) = p
463 CALL sswap( n, z( 1, i ), 1, z( 1, k ), 1 )
464 END IF
465 40 CONTINUE
466 END IF
467 END IF
468*
469 50 CONTINUE
470 work( 1 ) = sroundup_lwork(lwmin)
471 iwork( 1 ) = liwmin
472*
473 RETURN
474*
475* End of SSTEDC
476*
477 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
subroutine slaed0(icompq, qsiz, n, d, e, q, ldq, qstore, ldqs, work, iwork, info)
SLAED0 used by SSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition slaed0.f:170
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:86
subroutine sstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
SSTEDC
Definition sstedc.f:180
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:129
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:84
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82