LAPACK 3.12.0
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*> \htmlonly
9*> Download SSTEDC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstedc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstedc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstedc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSTEDC( 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* REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SSTEDC 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 SSYTRD or SSPTRD or SSBTRD has been used to reduce this
43*> matrix to tridiagonal form.
44*>
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] COMPZ
51*> \verbatim
52*> COMPZ is CHARACTER*1
53*> = 'N': Compute eigenvalues only.
54*> = 'I': Compute eigenvectors of tridiagonal matrix also.
55*> = 'V': Compute eigenvectors of original dense symmetric
56*> matrix also. On entry, Z contains the orthogonal
57*> matrix used to reduce the original matrix to
58*> tridiagonal form.
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*> N is INTEGER
64*> The dimension of the symmetric tridiagonal matrix. N >= 0.
65*> \endverbatim
66*>
67*> \param[in,out] D
68*> \verbatim
69*> D is REAL array, dimension (N)
70*> On entry, the diagonal elements of the tridiagonal matrix.
71*> On exit, if INFO = 0, the eigenvalues in ascending order.
72*> \endverbatim
73*>
74*> \param[in,out] E
75*> \verbatim
76*> E is REAL array, dimension (N-1)
77*> On entry, the subdiagonal elements of the tridiagonal matrix.
78*> On exit, E has been destroyed.
79*> \endverbatim
80*>
81*> \param[in,out] Z
82*> \verbatim
83*> Z is REAL array, dimension (LDZ,N)
84*> On entry, if COMPZ = 'V', then Z contains the orthogonal
85*> matrix used in the reduction to tridiagonal form.
86*> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
87*> orthonormal eigenvectors of the original symmetric matrix,
88*> and if COMPZ = 'I', Z contains the orthonormal eigenvectors
89*> of the symmetric tridiagonal matrix.
90*> If COMPZ = 'N', then Z is not referenced.
91*> \endverbatim
92*>
93*> \param[in] LDZ
94*> \verbatim
95*> LDZ is INTEGER
96*> The leading dimension of the array Z. LDZ >= 1.
97*> If eigenvectors are desired, then LDZ >= max(1,N).
98*> \endverbatim
99*>
100*> \param[out] WORK
101*> \verbatim
102*> WORK is REAL array, dimension (MAX(1,LWORK))
103*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
104*> \endverbatim
105*>
106*> \param[in] LWORK
107*> \verbatim
108*> LWORK is INTEGER
109*> The dimension of the array WORK.
110*> If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
111*> If COMPZ = 'V' and N > 1 then LWORK must be at least
112*> ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
113*> where lg( N ) = smallest integer k such
114*> that 2**k >= N.
115*> If COMPZ = 'I' and N > 1 then LWORK must be at least
116*> ( 1 + 4*N + N**2 ).
117*> Note that for COMPZ = 'I' or 'V', then if N is less than or
118*> equal to the minimum divide size, usually 25, then LWORK need
119*> only be max(1,2*(N-1)).
120*>
121*> If LWORK = -1, then a workspace query is assumed; the routine
122*> only calculates the optimal size of the WORK array, returns
123*> this value as the first entry of the WORK array, and no error
124*> message related to LWORK is issued by XERBLA.
125*> \endverbatim
126*>
127*> \param[out] IWORK
128*> \verbatim
129*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
130*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
131*> \endverbatim
132*>
133*> \param[in] LIWORK
134*> \verbatim
135*> LIWORK is INTEGER
136*> The dimension of the array IWORK.
137*> If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
138*> If COMPZ = 'V' and N > 1 then LIWORK must be at least
139*> ( 6 + 6*N + 5*N*lg N ).
140*> If COMPZ = 'I' and N > 1 then LIWORK must be at least
141*> ( 3 + 5*N ).
142*> Note that for COMPZ = 'I' or 'V', then if N is less than or
143*> equal to the minimum divide size, usually 25, then LIWORK
144*> need only be 1.
145*>
146*> If LIWORK = -1, then a workspace query is assumed; the
147*> routine only calculates the optimal size of the IWORK array,
148*> returns this value as the first entry of the IWORK array, and
149*> no error message related to LIWORK is issued by XERBLA.
150*> \endverbatim
151*>
152*> \param[out] INFO
153*> \verbatim
154*> INFO is INTEGER
155*> = 0: successful exit.
156*> < 0: if INFO = -i, the i-th argument had an illegal value.
157*> > 0: The algorithm failed to compute an eigenvalue while
158*> working on the submatrix lying in rows and columns
159*> INFO/(N+1) through mod(INFO,N+1).
160*> \endverbatim
161*
162* Authors:
163* ========
164*
165*> \author Univ. of Tennessee
166*> \author Univ. of California Berkeley
167*> \author Univ. of Colorado Denver
168*> \author NAG Ltd.
169*
170*> \ingroup stedc
171*
172*> \par Contributors:
173* ==================
174*>
175*> Jeff Rutter, Computer Science Division, University of California
176*> at Berkeley, USA \n
177*> Modified by Francoise Tisseur, University of Tennessee
178*>
179* =====================================================================
180 SUBROUTINE sstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
181 $ LIWORK, INFO )
182*
183* -- LAPACK computational routine --
184* -- LAPACK is a software package provided by Univ. of Tennessee, --
185* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187* .. Scalar Arguments ..
188 CHARACTER COMPZ
189 INTEGER INFO, LDZ, LIWORK, LWORK, N
190* ..
191* .. Array Arguments ..
192 INTEGER IWORK( * )
193 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 REAL ZERO, ONE, TWO
200 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
201* ..
202* .. Local Scalars ..
203 LOGICAL LQUERY
204 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
205 $ lwmin, m, smlsiz, start, storez, strtrw
206 REAL EPS, ORGNRM, P, TINY
207* ..
208* .. External Functions ..
209 LOGICAL LSAME
210 INTEGER ILAENV
211 REAL SLAMCH, SLANST, SROUNDUP_LWORK
212 EXTERNAL ilaenv, lsame, slamch, slanst, sroundup_lwork
213* ..
214* .. External Subroutines ..
215 EXTERNAL sgemm, slacpy, slaed0, slascl, slaset, 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 ), m,
381 $ info )
382 CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
383 $ m-1, info )
384*
385 IF( icompz.EQ.1 ) THEN
386 strtrw = 1
387 ELSE
388 strtrw = start
389 END IF
390 CALL slaed0( icompz, n, m, d( start ), e( start ),
391 $ z( strtrw, start ), ldz, work( 1 ), n,
392 $ work( storez ), iwork, info )
393 IF( info.NE.0 ) THEN
394 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
395 $ mod( info, ( m+1 ) ) + start - 1
396 GO TO 50
397 END IF
398*
399* Scale back.
400*
401 CALL slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
402 $ info )
403*
404 ELSE
405 IF( icompz.EQ.1 ) THEN
406*
407* Since QR won't update a Z matrix which is larger than
408* the length of D, we must solve the sub-problem in a
409* workspace and then multiply back into Z.
410*
411 CALL ssteqr( 'I', m, d( start ), e( start ), work, m,
412 $ work( m*m+1 ), info )
413 CALL slacpy( 'A', n, m, z( 1, start ), ldz,
414 $ work( storez ), n )
415 CALL sgemm( 'N', 'N', n, m, m, one,
416 $ work( storez ), n, work, m, zero,
417 $ z( 1, start ), ldz )
418 ELSE IF( icompz.EQ.2 ) THEN
419 CALL ssteqr( 'I', m, d( start ), e( start ),
420 $ z( start, start ), ldz, work, info )
421 ELSE
422 CALL ssterf( m, d( start ), e( start ), info )
423 END IF
424 IF( info.NE.0 ) THEN
425 info = start*( n+1 ) + finish
426 GO TO 50
427 END IF
428 END IF
429*
430 start = finish + 1
431 GO TO 10
432 END IF
433*
434* endwhile
435*
436 IF( icompz.EQ.0 ) THEN
437*
438* Use Quick Sort
439*
440 CALL slasrt( 'I', n, d, info )
441*
442 ELSE
443*
444* Use Selection Sort to minimize swaps of eigenvectors
445*
446 DO 40 ii = 2, n
447 i = ii - 1
448 k = i
449 p = d( i )
450 DO 30 j = ii, n
451 IF( d( j ).LT.p ) THEN
452 k = j
453 p = d( j )
454 END IF
455 30 CONTINUE
456 IF( k.NE.i ) THEN
457 d( k ) = d( i )
458 d( i ) = p
459 CALL sswap( n, z( 1, i ), 1, z( 1, k ), 1 )
460 END IF
461 40 CONTINUE
462 END IF
463 END IF
464*
465 50 CONTINUE
466 work( 1 ) = sroundup_lwork(lwmin)
467 iwork( 1 ) = liwmin
468*
469 RETURN
470*
471* End of SSTEDC
472*
473 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:103
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:172
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:143
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:110
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:88
subroutine sstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
SSTEDC
Definition sstedc.f:182
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:131
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82