LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ sstedc()

subroutine sstedc ( character  COMPZ,
integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
integer  INFO 
)

SSTEDC

Download SSTEDC + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SSTEDC computes all eigenvalues and, optionally, eigenvectors of a
 symmetric tridiagonal matrix using the divide and conquer method.
 The eigenvectors of a full or band real symmetric matrix can also be
 found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this
 matrix to tridiagonal form.

 This code makes very mild assumptions about floating point
 arithmetic. It will work on machines with a guard digit in
 add/subtract, or on those binary machines without guard digits
 which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
 It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.  See SLAED3 for details.
Parameters
[in]COMPZ
          COMPZ is CHARACTER*1
          = 'N':  Compute eigenvalues only.
          = 'I':  Compute eigenvectors of tridiagonal matrix also.
          = 'V':  Compute eigenvectors of original dense symmetric
                  matrix also.  On entry, Z contains the orthogonal
                  matrix used to reduce the original matrix to
                  tridiagonal form.
[in]N
          N is INTEGER
          The dimension of the symmetric tridiagonal matrix.  N >= 0.
[in,out]D
          D is REAL array, dimension (N)
          On entry, the diagonal elements of the tridiagonal matrix.
          On exit, if INFO = 0, the eigenvalues in ascending order.
[in,out]E
          E is REAL array, dimension (N-1)
          On entry, the subdiagonal elements of the tridiagonal matrix.
          On exit, E has been destroyed.
[in,out]Z
          Z is REAL array, dimension (LDZ,N)
          On entry, if COMPZ = 'V', then Z contains the orthogonal
          matrix used in the reduction to tridiagonal form.
          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
          orthonormal eigenvectors of the original symmetric matrix,
          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
          of the symmetric tridiagonal matrix.
          If  COMPZ = 'N', then Z is not referenced.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1.
          If eigenvectors are desired, then LDZ >= max(1,N).
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
          If COMPZ = 'V' and N > 1 then LWORK must be at least
                         ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
                         where lg( N ) = smallest integer k such
                         that 2**k >= N.
          If COMPZ = 'I' and N > 1 then LWORK must be at least
                         ( 1 + 4*N + N**2 ).
          Note that for COMPZ = 'I' or 'V', then if N is less than or
          equal to the minimum divide size, usually 25, then LWORK need
          only be max(1,2*(N-1)).

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK.
          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
          If COMPZ = 'V' and N > 1 then LIWORK must be at least
                         ( 6 + 6*N + 5*N*lg N ).
          If COMPZ = 'I' and N > 1 then LIWORK must be at least
                         ( 3 + 5*N ).
          Note that for COMPZ = 'I' or 'V', then if N is less than or
          equal to the minimum divide size, usually 25, then LIWORK
          need only be 1.

          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates the optimal size of the IWORK array,
          returns this value as the first entry of the IWORK array, and
          no error message related to LIWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  The algorithm failed to compute an eigenvalue while
                working on the submatrix lying in rows and columns
                INFO/(N+1) through mod(INFO,N+1).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

Definition at line 190 of file sstedc.f.

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,
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 *
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:90
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
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
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
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
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
real function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
Definition: slanst.f:102
Here is the call graph for this function:
Here is the caller graph for this function: