LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
sstevr.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine sstevr (JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK, LIWORK, INFO)
  SSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices More...
 

Function/Subroutine Documentation

subroutine sstevr ( character  JOBZ,
character  RANGE,
integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
real  VL,
real  VU,
integer  IL,
integer  IU,
real  ABSTOL,
integer  M,
real, dimension( * )  W,
real, dimension( ldz, * )  Z,
integer  LDZ,
integer, dimension( * )  ISUPPZ,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
integer  INFO 
)

SSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

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

Purpose:
 SSTEVR computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric tridiagonal matrix T.  Eigenvalues and
 eigenvectors can be selected by specifying either a range of values
 or a range of indices for the desired eigenvalues.

 Whenever possible, SSTEVR calls SSTEMR to compute the
 eigenspectrum using Relatively Robust Representations.  SSTEMR
 computes eigenvalues by the dqds algorithm, while orthogonal
 eigenvectors are computed from various "good" L D L^T representations
 (also known as Relatively Robust Representations). Gram-Schmidt
 orthogonalization is avoided as far as possible. More specifically,
 the various steps of the algorithm are as follows. For the i-th
 unreduced block of T,
    (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
         is a relatively robust representation,
    (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
        relative accuracy by the dqds algorithm,
    (c) If there is a cluster of close eigenvalues, "choose" sigma_i
        close to the cluster, and go to step (a),
    (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
        compute the corresponding eigenvector by forming a
        rank-revealing twisted factorization.
 The desired accuracy of the output can be specified by the input
 parameter ABSTOL.

 For more details, see "A new O(n^2) algorithm for the symmetric
 tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
 Computer Science Division Technical Report No. UCB//CSD-97-971,
 UC Berkeley, May 1997.


 Note 1 : SSTEVR calls SSTEMR when the full spectrum is requested
 on machines which conform to the ieee-754 floating point standard.
 SSTEVR calls SSTEBZ and SSTEIN on non-ieee machines and
 when partial spectrum requests are made.

 Normal execution of SSTEMR may create NaNs and infinities and
 hence may abort due to a floating point exception in environments
 which do not handle NaNs and infinities in the ieee standard default
 manner.
Parameters
[in]JOBZ
          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': all eigenvalues will be found.
          = 'V': all eigenvalues in the half-open interval (VL,VU]
                 will be found.
          = 'I': the IL-th through IU-th eigenvalues will be found.
          For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
          SSTEIN are called
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in,out]D
          D is REAL array, dimension (N)
          On entry, the n diagonal elements of the tridiagonal matrix
          A.
          On exit, D may be multiplied by a constant factor chosen
          to avoid over/underflow in computing the eigenvalues.
[in,out]E
          E is REAL array, dimension (max(1,N-1))
          On entry, the (n-1) subdiagonal elements of the tridiagonal
          matrix A in elements 1 to N-1 of E.
          On exit, E may be multiplied by a constant factor chosen
          to avoid over/underflow in computing the eigenvalues.
[in]VL
          VL is REAL
[in]VU
          VU is REAL
          If RANGE='V', the lower and upper bounds of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
[in]IU
          IU is INTEGER
          If RANGE='I', the indices (in ascending order) of the
          smallest and largest eigenvalues to be returned.
          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]ABSTOL
          ABSTOL is REAL
          The absolute error tolerance for the eigenvalues.
          An approximate eigenvalue is accepted as converged
          when it is determined to lie in an interval [a,b]
          of width less than or equal to

                  ABSTOL + EPS *   max( |a|,|b| ) ,

          where EPS is the machine precision.  If ABSTOL is less than
          or equal to zero, then  EPS*|T|  will be used in its place,
          where |T| is the 1-norm of the tridiagonal matrix obtained
          by reducing A to tridiagonal form.

          See "Computing Small Singular Values of Bidiagonal Matrices
          with Guaranteed High Relative Accuracy," by Demmel and
          Kahan, LAPACK Working Note #3.

          If high relative accuracy is important, set ABSTOL to
          SLAMCH( 'Safe minimum' ).  Doing so will guarantee that
          eigenvalues are computed to high relative accuracy when
          possible in future releases.  The current code does not
          make any guarantees about high relative accuracy, but
          future releases will. See J. Barlow and J. Demmel,
          "Computing Accurate Eigensystems of Scaled Diagonally
          Dominant Matrices", LAPACK Working Note #7, for a discussion
          of which matrices define their eigenvalues to high relative
          accuracy.
[out]M
          M is INTEGER
          The total number of eigenvalues found.  0 <= M <= N.
          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
[out]W
          W is REAL array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.
[out]Z
          Z is REAL array, dimension (LDZ, max(1,M) )
          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
          contain the orthonormal eigenvectors of the matrix A
          corresponding to the selected eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z; if RANGE = 'V', the exact value of M
          is not known in advance and an upper bound must be used.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(1,N).
[out]ISUPPZ
          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
          The support of the eigenvectors in Z, i.e., the indices
          indicating the nonzero elements in Z. The i-th eigenvector
          is nonzero only in elements ISUPPZ( 2*i-1 ) through
          ISUPPZ( 2*i ).
          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal (and
          minimal) LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= 20*N.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal sizes of the WORK and IWORK
          arrays, returns these values as the first entries of the WORK
          and IWORK arrays, and no error message related to LWORK or
          LIWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
          On exit, if INFO = 0, IWORK(1) returns the optimal (and
          minimal) LIWORK.
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK.  LIWORK >= 10*N.

          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates the optimal sizes of the WORK and
          IWORK arrays, returns these values as the first entries of
          the WORK and IWORK arrays, and no error message related to
          LWORK or 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:  Internal error
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Contributors:
Inderjit Dhillon, IBM Almaden, USA
Osni Marques, LBNL/NERSC, USA
Ken Stanley, Computer Science Division, University of California at Berkeley, USA
Jason Riedy, Computer Science Division, University of California at Berkeley, USA

Definition at line 301 of file sstevr.f.

301 *
302 * -- LAPACK driver routine (version 3.4.0) --
303 * -- LAPACK is a software package provided by Univ. of Tennessee, --
304 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
305 * November 2011
306 *
307 * .. Scalar Arguments ..
308  CHARACTER jobz, range
309  INTEGER il, info, iu, ldz, liwork, lwork, m, n
310  REAL abstol, vl, vu
311 * ..
312 * .. Array Arguments ..
313  INTEGER isuppz( * ), iwork( * )
314  REAL d( * ), e( * ), w( * ), work( * ), z( ldz, * )
315 * ..
316 *
317 * =====================================================================
318 *
319 * .. Parameters ..
320  REAL zero, one, two
321  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
322 * ..
323 * .. Local Scalars ..
324  LOGICAL alleig, indeig, test, lquery, valeig, wantz,
325  $ tryrac
326  CHARACTER order
327  INTEGER i, ieeeok, imax, indibl, indifl, indisp,
328  $ indiwo, iscale, j, jj, liwmin, lwmin, nsplit
329  REAL bignum, eps, rmax, rmin, safmin, sigma, smlnum,
330  $ tmp1, tnrm, vll, vuu
331 * ..
332 * .. External Functions ..
333  LOGICAL lsame
334  INTEGER ilaenv
335  REAL slamch, slanst
336  EXTERNAL lsame, ilaenv, slamch, slanst
337 * ..
338 * .. External Subroutines ..
339  EXTERNAL scopy, sscal, sstebz, sstemr, sstein, ssterf,
340  $ sswap, xerbla
341 * ..
342 * .. Intrinsic Functions ..
343  INTRINSIC max, min, sqrt
344 * ..
345 * .. Executable Statements ..
346 *
347 *
348 * Test the input parameters.
349 *
350  ieeeok = ilaenv( 10, 'SSTEVR', 'N', 1, 2, 3, 4 )
351 *
352  wantz = lsame( jobz, 'V' )
353  alleig = lsame( range, 'A' )
354  valeig = lsame( range, 'V' )
355  indeig = lsame( range, 'I' )
356 *
357  lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
358  lwmin = max( 1, 20*n )
359  liwmin = max(1, 10*n )
360 *
361 *
362  info = 0
363  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
364  info = -1
365  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
366  info = -2
367  ELSE IF( n.LT.0 ) THEN
368  info = -3
369  ELSE
370  IF( valeig ) THEN
371  IF( n.GT.0 .AND. vu.LE.vl )
372  $ info = -7
373  ELSE IF( indeig ) THEN
374  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
375  info = -8
376  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
377  info = -9
378  END IF
379  END IF
380  END IF
381  IF( info.EQ.0 ) THEN
382  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
383  info = -14
384  END IF
385  END IF
386 *
387  IF( info.EQ.0 ) THEN
388  work( 1 ) = lwmin
389  iwork( 1 ) = liwmin
390 *
391  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
392  info = -17
393  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
394  info = -19
395  END IF
396  END IF
397 *
398  IF( info.NE.0 ) THEN
399  CALL xerbla( 'SSTEVR', -info )
400  RETURN
401  ELSE IF( lquery ) THEN
402  RETURN
403  END IF
404 *
405 * Quick return if possible
406 *
407  m = 0
408  IF( n.EQ.0 )
409  $ RETURN
410 *
411  IF( n.EQ.1 ) THEN
412  IF( alleig .OR. indeig ) THEN
413  m = 1
414  w( 1 ) = d( 1 )
415  ELSE
416  IF( vl.LT.d( 1 ) .AND. vu.GE.d( 1 ) ) THEN
417  m = 1
418  w( 1 ) = d( 1 )
419  END IF
420  END IF
421  IF( wantz )
422  $ z( 1, 1 ) = one
423  RETURN
424  END IF
425 *
426 * Get machine constants.
427 *
428  safmin = slamch( 'Safe minimum' )
429  eps = slamch( 'Precision' )
430  smlnum = safmin / eps
431  bignum = one / smlnum
432  rmin = sqrt( smlnum )
433  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
434 *
435 *
436 * Scale matrix to allowable range, if necessary.
437 *
438  iscale = 0
439  vll = vl
440  vuu = vu
441 *
442  tnrm = slanst( 'M', n, d, e )
443  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
444  iscale = 1
445  sigma = rmin / tnrm
446  ELSE IF( tnrm.GT.rmax ) THEN
447  iscale = 1
448  sigma = rmax / tnrm
449  END IF
450  IF( iscale.EQ.1 ) THEN
451  CALL sscal( n, sigma, d, 1 )
452  CALL sscal( n-1, sigma, e( 1 ), 1 )
453  IF( valeig ) THEN
454  vll = vl*sigma
455  vuu = vu*sigma
456  END IF
457  END IF
458 
459 * Initialize indices into workspaces. Note: These indices are used only
460 * if SSTERF or SSTEMR fail.
461 
462 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and
463 * stores the block indices of each of the M<=N eigenvalues.
464  indibl = 1
465 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and
466 * stores the starting and finishing indices of each block.
467  indisp = indibl + n
468 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
469 * that corresponding to eigenvectors that fail to converge in
470 * SSTEIN. This information is discarded; if any fail, the driver
471 * returns INFO > 0.
472  indifl = indisp + n
473 * INDIWO is the offset of the remaining integer workspace.
474  indiwo = indisp + n
475 *
476 * If all eigenvalues are desired, then
477 * call SSTERF or SSTEMR. If this fails for some eigenvalue, then
478 * try SSTEBZ.
479 *
480 *
481  test = .false.
482  IF( indeig ) THEN
483  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
484  test = .true.
485  END IF
486  END IF
487  IF( ( alleig .OR. test ) .AND. ieeeok.EQ.1 ) THEN
488  CALL scopy( n-1, e( 1 ), 1, work( 1 ), 1 )
489  IF( .NOT.wantz ) THEN
490  CALL scopy( n, d, 1, w, 1 )
491  CALL ssterf( n, w, work, info )
492  ELSE
493  CALL scopy( n, d, 1, work( n+1 ), 1 )
494  IF (abstol .LE. two*n*eps) THEN
495  tryrac = .true.
496  ELSE
497  tryrac = .false.
498  END IF
499  CALL sstemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,
500  $ iu, m, w, z, ldz, n, isuppz, tryrac,
501  $ work( 2*n+1 ), lwork-2*n, iwork, liwork, info )
502 *
503  END IF
504  IF( info.EQ.0 ) THEN
505  m = n
506  GO TO 10
507  END IF
508  info = 0
509  END IF
510 *
511 * Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
512 *
513  IF( wantz ) THEN
514  order = 'B'
515  ELSE
516  order = 'E'
517  END IF
518 
519  CALL sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,
520  $ nsplit, w, iwork( indibl ), iwork( indisp ), work,
521  $ iwork( indiwo ), info )
522 *
523  IF( wantz ) THEN
524  CALL sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),
525  $ z, ldz, work, iwork( indiwo ), iwork( indifl ),
526  $ info )
527  END IF
528 *
529 * If matrix was scaled, then rescale eigenvalues appropriately.
530 *
531  10 CONTINUE
532  IF( iscale.EQ.1 ) THEN
533  IF( info.EQ.0 ) THEN
534  imax = m
535  ELSE
536  imax = info - 1
537  END IF
538  CALL sscal( imax, one / sigma, w, 1 )
539  END IF
540 *
541 * If eigenvalues are not in order, then sort them, along with
542 * eigenvectors.
543 *
544  IF( wantz ) THEN
545  DO 30 j = 1, m - 1
546  i = 0
547  tmp1 = w( j )
548  DO 20 jj = j + 1, m
549  IF( w( jj ).LT.tmp1 ) THEN
550  i = jj
551  tmp1 = w( jj )
552  END IF
553  20 CONTINUE
554 *
555  IF( i.NE.0 ) THEN
556  w( i ) = w( j )
557  w( j ) = tmp1
558  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
559  END IF
560  30 CONTINUE
561  END IF
562 *
563 * Causes problems with tests 19 & 20:
564 * IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002
565 *
566 *
567  work( 1 ) = lwmin
568  iwork( 1 ) = liwmin
569  RETURN
570 *
571 * End of SSTEVR
572 *
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
Definition: sstein.f:176
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
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:265
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine sstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEMR
Definition: sstemr.f:314
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: