LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 *> 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 SLAED3 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 REAL 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 REAL 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 REAL 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 REAL 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 December 2016
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 sstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
189  $ LIWORK, INFO )
190 *
191 * -- LAPACK computational routine (version 3.7.0) --
192 * -- LAPACK is a software package provided by Univ. of Tennessee, --
193 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194 * December 2016
195 *
196 * .. Scalar Arguments ..
197  CHARACTER COMPZ
198  INTEGER INFO, LDZ, LIWORK, LWORK, N
199 * ..
200 * .. Array Arguments ..
201  INTEGER IWORK( * )
202  REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  REAL ZERO, ONE, TWO
209  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
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  REAL EPS, ORGNRM, P, TINY
216 * ..
217 * .. External Functions ..
218  LOGICAL LSAME
219  INTEGER ILAENV
220  REAL SLAMCH, SLANST
221  EXTERNAL ilaenv, lsame, slamch, slanst
222 * ..
223 * .. External Subroutines ..
224  EXTERNAL sgemm, slacpy, slaed0, slascl, slaset, slasrt,
225  $ ssteqr, ssterf, sswap, xerbla
226 * ..
227 * .. Intrinsic Functions ..
228  INTRINSIC abs, int, log, max, mod, real, 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, 'SSTEDC', ' ', 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( real( 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( 'SSTEDC', -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 SSTERF 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 SSTERF to compute the eigenvalues.
317 *
318  IF( icompz.EQ.0 ) THEN
319  CALL ssterf( 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 ssteqr( 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 slaset( 'Full', n, n, zero, one, z, ldz )
343  END IF
344 *
345 * Scale.
346 *
347  orgnrm = slanst( 'M', n, d, e )
348  IF( orgnrm.EQ.zero )
349  $ GO TO 50
350 *
351  eps = slamch( '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 = slanst( 'M', m, d( start ), e( start ) )
389  CALL slascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
390  $ info )
391  CALL slascl( '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 slaed0( 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 slascl( '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 ssteqr( 'I', m, d( start ), e( start ), work, m,
421  $ work( m*m+1 ), info )
422  CALL slacpy( 'A', n, m, z( 1, start ), ldz,
423  $ work( storez ), n )
424  CALL sgemm( '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 ssteqr( 'I', m, d( start ), e( start ),
429  $ z( start, start ), ldz, work, info )
430  ELSE
431  CALL ssterf( 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 slasrt( '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 sswap( 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 SSTEDC
481 *
482  END
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
sstedc
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:190
slascl
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:145
slacpy
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
ssterf
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
slaset
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:112
ssteqr
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
slaed0
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:174
slasrt
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:90