LAPACK  3.10.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.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

Definition at line 186 of file sstedc.f.

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  REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
200 * ..
201 *
202 * =====================================================================
203 *
204 * .. Parameters ..
205  REAL ZERO, ONE, TWO
206  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
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  REAL EPS, ORGNRM, P, TINY
213 * ..
214 * .. External Functions ..
215  LOGICAL LSAME
216  INTEGER ILAENV
217  REAL SLAMCH, SLANST
218  EXTERNAL ilaenv, lsame, slamch, slanst
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL sgemm, slacpy, slaed0, slascl, slaset, slasrt,
223 * ..
224 * .. Intrinsic Functions ..
225  INTRINSIC abs, int, log, max, mod, real, 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, 'SSTEDC', ' ', 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( real( 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( 'SSTEDC', -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 SSTERF 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 SSTERF to compute the eigenvalues.
314 *
315  IF( icompz.EQ.0 ) THEN
316  CALL ssterf( 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 ssteqr( 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 slaset( 'Full', n, n, zero, one, z, ldz )
340  END IF
341 *
342 * Scale.
343 *
344  orgnrm = slanst( 'M', n, d, e )
345  IF( orgnrm.EQ.zero )
346  $ GO TO 50
347 *
348  eps = slamch( '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 = slanst( 'M', m, d( start ), e( start ) )
386  CALL slascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
387  $ info )
388  CALL slascl( '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 slaed0( 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 slascl( '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 ssteqr( 'I', m, d( start ), e( start ), work, m,
418  $ work( m*m+1 ), info )
419  CALL slacpy( 'A', n, m, z( 1, start ), ldz,
420  $ work( storez ), n )
421  CALL sgemm( '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 ssteqr( 'I', m, d( start ), e( start ),
426  $ z( start, start ), ldz, work, info )
427  ELSE
428  CALL ssterf( 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 slasrt( '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 sswap( 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 SSTEDC
478 *
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
real function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: slanst.f:100
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:131
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:88
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
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 sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: