LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dspevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
double precision, dimension( * )  AP,
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 
)

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

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

Purpose:
 DSPEVX computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric matrix A in packed storage.  Eigenvalues/vectors
 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]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]AP
          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
          On entry, the upper or lower triangle of the symmetric matrix
          A, packed columnwise in a linear array.  The j-th column of A
          is stored in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.

          On exit, AP is overwritten by values generated during the
          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
          and first superdiagonal of the tridiagonal matrix T overwrite
          the corresponding elements of A, and if UPLO = 'L', the
          diagonal and first subdiagonal of T overwrite the
          corresponding elements of A.
[in]VL
          VL is DOUBLE PRECISION
          If RANGE='V', the lower bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]VU
          VU is DOUBLE PRECISION
          If RANGE='V', the upper bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
          If RANGE='I', the index of the
          smallest eigenvalue 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]IU
          IU is INTEGER
          If RANGE='I', the index of the
          largest eigenvalue 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 AP to tridiagonal form.

          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)
          If INFO = 0, 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, 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 (8*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
June 2016

Definition at line 236 of file dspevx.f.

236 *
237 * -- LAPACK driver routine (version 3.6.1) --
238 * -- LAPACK is a software package provided by Univ. of Tennessee, --
239 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240 * June 2016
241 *
242 * .. Scalar Arguments ..
243  CHARACTER jobz, range, uplo
244  INTEGER il, info, iu, ldz, m, n
245  DOUBLE PRECISION abstol, vl, vu
246 * ..
247 * .. Array Arguments ..
248  INTEGER ifail( * ), iwork( * )
249  DOUBLE PRECISION ap( * ), w( * ), work( * ), z( ldz, * )
250 * ..
251 *
252 * =====================================================================
253 *
254 * .. Parameters ..
255  DOUBLE PRECISION zero, one
256  parameter ( zero = 0.0d0, one = 1.0d0 )
257 * ..
258 * .. Local Scalars ..
259  LOGICAL alleig, indeig, test, valeig, wantz
260  CHARACTER order
261  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
262  $ indisp, indiwo, indtau, indwrk, iscale, itmp1,
263  $ j, jj, nsplit
264  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
265  $ sigma, smlnum, tmp1, vll, vuu
266 * ..
267 * .. External Functions ..
268  LOGICAL lsame
269  DOUBLE PRECISION dlamch, dlansp
270  EXTERNAL lsame, dlamch, dlansp
271 * ..
272 * .. External Subroutines ..
273  EXTERNAL dcopy, dopgtr, dopmtr, dscal, dsptrd, dstebz,
275 * ..
276 * .. Intrinsic Functions ..
277  INTRINSIC max, min, sqrt
278 * ..
279 * .. Executable Statements ..
280 *
281 * Test the input parameters.
282 *
283  wantz = lsame( jobz, 'V' )
284  alleig = lsame( range, 'A' )
285  valeig = lsame( range, 'V' )
286  indeig = lsame( range, 'I' )
287 *
288  info = 0
289  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
290  info = -1
291  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
292  info = -2
293  ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
294  $ THEN
295  info = -3
296  ELSE IF( n.LT.0 ) THEN
297  info = -4
298  ELSE
299  IF( valeig ) THEN
300  IF( n.GT.0 .AND. vu.LE.vl )
301  $ info = -7
302  ELSE IF( indeig ) THEN
303  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
304  info = -8
305  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
306  info = -9
307  END IF
308  END IF
309  END IF
310  IF( info.EQ.0 ) THEN
311  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
312  $ info = -14
313  END IF
314 *
315  IF( info.NE.0 ) THEN
316  CALL xerbla( 'DSPEVX', -info )
317  RETURN
318  END IF
319 *
320 * Quick return if possible
321 *
322  m = 0
323  IF( n.EQ.0 )
324  $ RETURN
325 *
326  IF( n.EQ.1 ) THEN
327  IF( alleig .OR. indeig ) THEN
328  m = 1
329  w( 1 ) = ap( 1 )
330  ELSE
331  IF( vl.LT.ap( 1 ) .AND. vu.GE.ap( 1 ) ) THEN
332  m = 1
333  w( 1 ) = ap( 1 )
334  END IF
335  END IF
336  IF( wantz )
337  $ z( 1, 1 ) = one
338  RETURN
339  END IF
340 *
341 * Get machine constants.
342 *
343  safmin = dlamch( 'Safe minimum' )
344  eps = dlamch( 'Precision' )
345  smlnum = safmin / eps
346  bignum = one / smlnum
347  rmin = sqrt( smlnum )
348  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
349 *
350 * Scale matrix to allowable range, if necessary.
351 *
352  iscale = 0
353  abstll = abstol
354  IF( valeig ) THEN
355  vll = vl
356  vuu = vu
357  ELSE
358  vll = zero
359  vuu = zero
360  END IF
361  anrm = dlansp( 'M', uplo, n, ap, work )
362  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
363  iscale = 1
364  sigma = rmin / anrm
365  ELSE IF( anrm.GT.rmax ) THEN
366  iscale = 1
367  sigma = rmax / anrm
368  END IF
369  IF( iscale.EQ.1 ) THEN
370  CALL dscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
371  IF( abstol.GT.0 )
372  $ abstll = abstol*sigma
373  IF( valeig ) THEN
374  vll = vl*sigma
375  vuu = vu*sigma
376  END IF
377  END IF
378 *
379 * Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
380 *
381  indtau = 1
382  inde = indtau + n
383  indd = inde + n
384  indwrk = indd + n
385  CALL dsptrd( uplo, n, ap, work( indd ), work( inde ),
386  $ work( indtau ), iinfo )
387 *
388 * If all eigenvalues are desired and ABSTOL is less than or equal
389 * to zero, then call DSTERF or DOPGTR and SSTEQR. If this fails
390 * for some eigenvalue, then try DSTEBZ.
391 *
392  test = .false.
393  IF (indeig) THEN
394  IF (il.EQ.1 .AND. iu.EQ.n) THEN
395  test = .true.
396  END IF
397  END IF
398  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
399  CALL dcopy( n, work( indd ), 1, w, 1 )
400  indee = indwrk + 2*n
401  IF( .NOT.wantz ) THEN
402  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
403  CALL dsterf( n, w, work( indee ), info )
404  ELSE
405  CALL dopgtr( uplo, n, ap, work( indtau ), z, ldz,
406  $ work( indwrk ), iinfo )
407  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
408  CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
409  $ work( indwrk ), info )
410  IF( info.EQ.0 ) THEN
411  DO 10 i = 1, n
412  ifail( i ) = 0
413  10 CONTINUE
414  END IF
415  END IF
416  IF( info.EQ.0 ) THEN
417  m = n
418  GO TO 20
419  END IF
420  info = 0
421  END IF
422 *
423 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
424 *
425  IF( wantz ) THEN
426  order = 'B'
427  ELSE
428  order = 'E'
429  END IF
430  indibl = 1
431  indisp = indibl + n
432  indiwo = indisp + n
433  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
434  $ work( indd ), work( inde ), m, nsplit, w,
435  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
436  $ iwork( indiwo ), info )
437 *
438  IF( wantz ) THEN
439  CALL dstein( n, work( indd ), work( inde ), m, w,
440  $ iwork( indibl ), iwork( indisp ), z, ldz,
441  $ work( indwrk ), iwork( indiwo ), ifail, info )
442 *
443 * Apply orthogonal matrix used in reduction to tridiagonal
444 * form to eigenvectors returned by DSTEIN.
445 *
446  CALL dopmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
447  $ work( indwrk ), iinfo )
448  END IF
449 *
450 * If matrix was scaled, then rescale eigenvalues appropriately.
451 *
452  20 CONTINUE
453  IF( iscale.EQ.1 ) THEN
454  IF( info.EQ.0 ) THEN
455  imax = m
456  ELSE
457  imax = info - 1
458  END IF
459  CALL dscal( imax, one / sigma, w, 1 )
460  END IF
461 *
462 * If eigenvalues are not in order, then sort them, along with
463 * eigenvectors.
464 *
465  IF( wantz ) THEN
466  DO 40 j = 1, m - 1
467  i = 0
468  tmp1 = w( j )
469  DO 30 jj = j + 1, m
470  IF( w( jj ).LT.tmp1 ) THEN
471  i = jj
472  tmp1 = w( jj )
473  END IF
474  30 CONTINUE
475 *
476  IF( i.NE.0 ) THEN
477  itmp1 = iwork( indibl+i-1 )
478  w( i ) = w( j )
479  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
480  w( j ) = tmp1
481  iwork( indibl+j-1 ) = itmp1
482  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
483  IF( info.NE.0 ) THEN
484  itmp1 = ifail( i )
485  ifail( i ) = ifail( j )
486  ifail( j ) = itmp1
487  END IF
488  END IF
489  40 CONTINUE
490  END IF
491 *
492  RETURN
493 *
494 * End of DSPEVX
495 *
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:275
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dopmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
DOPMTR
Definition: dopmtr.f:152
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:133
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dsptrd(UPLO, N, AP, D, E, TAU, INFO)
DSPTRD
Definition: dsptrd.f:152
subroutine dopgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
DOPGTR
Definition: dopgtr.f:116
double precision function dlansp(NORM, UPLO, N, AP, WORK)
DLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form.
Definition: dlansp.f:116
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:176
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: