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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

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
[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 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
November 2011

Definition at line 229 of file dspevx.f.

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