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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

subroutine sspevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
real, dimension( * )  AP,
real  VL,
real  VU,
integer  IL,
integer  IU,
real  ABSTOL,
integer  M,
real, dimension( * )  W,
real, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

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

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

Purpose:
 SSPEVX 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 REAL 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 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 AP to tridiagonal form.

          Eigenvalues will be computed most accurately when ABSTOL is
          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
          If this routine returns with INFO>0, indicating that some
          eigenvectors did not converge, try setting ABSTOL to
          2*SLAMCH('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 REAL array, dimension (N)
          If INFO = 0, 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).
          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 REAL 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 sspevx.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  REAL abstol, vl, vu
239 * ..
240 * .. Array Arguments ..
241  INTEGER ifail( * ), iwork( * )
242  REAL ap( * ), w( * ), work( * ), z( ldz, * )
243 * ..
244 *
245 * =====================================================================
246 *
247 * .. Parameters ..
248  REAL zero, one
249  parameter( zero = 0.0e0, one = 1.0e0 )
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  REAL abstll, anrm, bignum, eps, rmax, rmin, safmin,
258  $ sigma, smlnum, tmp1, vll, vuu
259 * ..
260 * .. External Functions ..
261  LOGICAL lsame
262  REAL slamch, slansp
263  EXTERNAL lsame, slamch, slansp
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL scopy, sopgtr, sopmtr, sscal, ssptrd, sstebz,
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( 'SSPEVX', -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 = slamch( 'Safe minimum' )
337  eps = slamch( '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  ENDIF
354  anrm = slansp( '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 sscal( ( 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 SSPTRD 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 ssptrd( 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 SSTERF or SOPGTR and SSTEQR. If this fails
383 * for some eigenvalue, then try SSTEBZ.
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 scopy( n, work( indd ), 1, w, 1 )
393  indee = indwrk + 2*n
394  IF( .NOT.wantz ) THEN
395  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
396  CALL ssterf( n, w, work( indee ), info )
397  ELSE
398  CALL sopgtr( uplo, n, ap, work( indtau ), z, ldz,
399  $ work( indwrk ), iinfo )
400  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
401  CALL ssteqr( 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 SSTEBZ 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 sstebz( 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 sstein( 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 SSTEIN.
438 *
439  CALL sopmtr( '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 sscal( 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 sswap( 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 SSPEVX
488 *
subroutine sopmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
SOPMTR
Definition: sopmtr.f:152
real function slansp(NORM, UPLO, N, AP, WORK)
SLANSP 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: slansp.f:116
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
subroutine ssptrd(UPLO, N, AP, D, E, TAU, INFO)
SSPTRD
Definition: ssptrd.f:152
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sopgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
SOPGTR
Definition: sopgtr.f:116
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
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
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
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: