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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

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

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

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

Purpose:
 DSTEVR 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, DSTEVR calls DSTEMR to compute the
 eigenspectrum using Relatively Robust Representations.  DSTEMR
 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 : DSTEVR calls DSTEMR when the full spectrum is requested
 on machines which conform to the ieee-754 floating point standard.
 DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
 when partial spectrum requests are made.

 Normal execution of DSTEMR 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, DSTEBZ and
          DSTEIN are called
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in,out]D
          D is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
[in]VU
          VU is DOUBLE PRECISION
          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 DOUBLE PRECISION
          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
          DLAMCH( '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 DOUBLE PRECISION array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.
[out]Z
          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 >= max(1,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 >= max(1,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

Definition at line 299 of file dstevr.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: