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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

subroutine dstevx ( 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,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

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

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

Purpose:
 DSTEVX computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric tridiagonal matrix A.  Eigenvalues and
 eigenvectors can be selected by specifying either a range of values
 or a range of indices for the desired eigenvalues.
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.
[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.

          Eigenvalues will be computed most accurately when ABSTOL is
          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
          If this routine returns with INFO>0, indicating that some
          eigenvectors did not converge, try setting ABSTOL to
          2*DLAMCH('S').

          See "Computing Small Singular Values of Bidiagonal Matrices
          with Guaranteed High Relative Accuracy," by Demmel and
          Kahan, LAPACK Working Note #3.
[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).
          If an eigenvector fails to converge (INFO > 0), then that
          column of Z contains the latest approximation to the
          eigenvector, and the index of the eigenvector is returned
          in IFAIL.  If JOBZ = 'N', then Z is not referenced.
          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]WORK
          WORK is DOUBLE PRECISION array, dimension (5*N)
[out]IWORK
          IWORK is INTEGER array, dimension (5*N)
[out]IFAIL
          IFAIL is INTEGER array, dimension (N)
          If JOBZ = 'V', then if INFO = 0, the first M elements of
          IFAIL are zero.  If INFO > 0, then IFAIL contains the
          indices of the eigenvectors that failed to converge.
          If JOBZ = 'N', then IFAIL is not referenced.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, then i eigenvectors failed to converge.
                Their indices are stored in array IFAIL.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 222 of file dstevx.f.

222 *
223 * -- LAPACK driver routine (version 3.4.0) --
224 * -- LAPACK is a software package provided by Univ. of Tennessee, --
225 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
226 * November 2011
227 *
228 * .. Scalar Arguments ..
229  CHARACTER jobz, range
230  INTEGER il, info, iu, ldz, m, n
231  DOUBLE PRECISION abstol, vl, vu
232 * ..
233 * .. Array Arguments ..
234  INTEGER ifail( * ), iwork( * )
235  DOUBLE PRECISION d( * ), e( * ), w( * ), work( * ), z( ldz, * )
236 * ..
237 *
238 * =====================================================================
239 *
240 * .. Parameters ..
241  DOUBLE PRECISION zero, one
242  parameter( zero = 0.0d0, one = 1.0d0 )
243 * ..
244 * .. Local Scalars ..
245  LOGICAL alleig, indeig, test, valeig, wantz
246  CHARACTER order
247  INTEGER i, imax, indibl, indisp, indiwo, indwrk,
248  $ iscale, itmp1, j, jj, nsplit
249  DOUBLE PRECISION bignum, eps, rmax, rmin, safmin, sigma, smlnum,
250  $ tmp1, tnrm, vll, vuu
251 * ..
252 * .. External Functions ..
253  LOGICAL lsame
254  DOUBLE PRECISION dlamch, dlanst
255  EXTERNAL lsame, dlamch, dlanst
256 * ..
257 * .. External Subroutines ..
258  EXTERNAL dcopy, dscal, dstebz, dstein, dsteqr, dsterf,
259  $ dswap, xerbla
260 * ..
261 * .. Intrinsic Functions ..
262  INTRINSIC max, min, sqrt
263 * ..
264 * .. Executable Statements ..
265 *
266 * Test the input parameters.
267 *
268  wantz = lsame( jobz, 'V' )
269  alleig = lsame( range, 'A' )
270  valeig = lsame( range, 'V' )
271  indeig = lsame( range, 'I' )
272 *
273  info = 0
274  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
275  info = -1
276  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
277  info = -2
278  ELSE IF( n.LT.0 ) THEN
279  info = -3
280  ELSE
281  IF( valeig ) THEN
282  IF( n.GT.0 .AND. vu.LE.vl )
283  $ info = -7
284  ELSE IF( indeig ) THEN
285  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
286  info = -8
287  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
288  info = -9
289  END IF
290  END IF
291  END IF
292  IF( info.EQ.0 ) THEN
293  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
294  $ info = -14
295  END IF
296 *
297  IF( info.NE.0 ) THEN
298  CALL xerbla( 'DSTEVX', -info )
299  RETURN
300  END IF
301 *
302 * Quick return if possible
303 *
304  m = 0
305  IF( n.EQ.0 )
306  $ RETURN
307 *
308  IF( n.EQ.1 ) THEN
309  IF( alleig .OR. indeig ) THEN
310  m = 1
311  w( 1 ) = d( 1 )
312  ELSE
313  IF( vl.LT.d( 1 ) .AND. vu.GE.d( 1 ) ) THEN
314  m = 1
315  w( 1 ) = d( 1 )
316  END IF
317  END IF
318  IF( wantz )
319  $ z( 1, 1 ) = one
320  RETURN
321  END IF
322 *
323 * Get machine constants.
324 *
325  safmin = dlamch( 'Safe minimum' )
326  eps = dlamch( 'Precision' )
327  smlnum = safmin / eps
328  bignum = one / smlnum
329  rmin = sqrt( smlnum )
330  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
331 *
332 * Scale matrix to allowable range, if necessary.
333 *
334  iscale = 0
335  IF( valeig ) THEN
336  vll = vl
337  vuu = vu
338  ELSE
339  vll = zero
340  vuu = zero
341  END IF
342  tnrm = dlanst( 'M', n, d, e )
343  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
344  iscale = 1
345  sigma = rmin / tnrm
346  ELSE IF( tnrm.GT.rmax ) THEN
347  iscale = 1
348  sigma = rmax / tnrm
349  END IF
350  IF( iscale.EQ.1 ) THEN
351  CALL dscal( n, sigma, d, 1 )
352  CALL dscal( n-1, sigma, e( 1 ), 1 )
353  IF( valeig ) THEN
354  vll = vl*sigma
355  vuu = vu*sigma
356  END IF
357  END IF
358 *
359 * If all eigenvalues are desired and ABSTOL is less than zero, then
360 * call DSTERF or SSTEQR. If this fails for some eigenvalue, then
361 * try DSTEBZ.
362 *
363  test = .false.
364  IF( indeig ) THEN
365  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
366  test = .true.
367  END IF
368  END IF
369  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
370  CALL dcopy( n, d, 1, w, 1 )
371  CALL dcopy( n-1, e( 1 ), 1, work( 1 ), 1 )
372  indwrk = n + 1
373  IF( .NOT.wantz ) THEN
374  CALL dsterf( n, w, work, info )
375  ELSE
376  CALL dsteqr( 'I', n, w, work, z, ldz, work( indwrk ), info )
377  IF( info.EQ.0 ) THEN
378  DO 10 i = 1, n
379  ifail( i ) = 0
380  10 CONTINUE
381  END IF
382  END IF
383  IF( info.EQ.0 ) THEN
384  m = n
385  GO TO 20
386  END IF
387  info = 0
388  END IF
389 *
390 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
391 *
392  IF( wantz ) THEN
393  order = 'B'
394  ELSE
395  order = 'E'
396  END IF
397  indwrk = 1
398  indibl = 1
399  indisp = indibl + n
400  indiwo = indisp + n
401  CALL dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,
402  $ nsplit, w, iwork( indibl ), iwork( indisp ),
403  $ work( indwrk ), iwork( indiwo ), info )
404 *
405  IF( wantz ) THEN
406  CALL dstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),
407  $ z, ldz, work( indwrk ), iwork( indiwo ), ifail,
408  $ info )
409  END IF
410 *
411 * If matrix was scaled, then rescale eigenvalues appropriately.
412 *
413  20 CONTINUE
414  IF( iscale.EQ.1 ) THEN
415  IF( info.EQ.0 ) THEN
416  imax = m
417  ELSE
418  imax = info - 1
419  END IF
420  CALL dscal( imax, one / sigma, w, 1 )
421  END IF
422 *
423 * If eigenvalues are not in order, then sort them, along with
424 * eigenvectors.
425 *
426  IF( wantz ) THEN
427  DO 40 j = 1, m - 1
428  i = 0
429  tmp1 = w( j )
430  DO 30 jj = j + 1, m
431  IF( w( jj ).LT.tmp1 ) THEN
432  i = jj
433  tmp1 = w( jj )
434  END IF
435  30 CONTINUE
436 *
437  IF( i.NE.0 ) THEN
438  itmp1 = iwork( indibl+i-1 )
439  w( i ) = w( j )
440  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
441  w( j ) = tmp1
442  iwork( indibl+j-1 ) = itmp1
443  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
444  IF( info.NE.0 ) THEN
445  itmp1 = ifail( i )
446  ifail( i ) = ifail( j )
447  ifail( j ) = itmp1
448  END IF
449  END IF
450  40 CONTINUE
451  END IF
452 *
453  RETURN
454 *
455 * End of DSTEVX
456 *
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:133
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
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
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: