LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ ssbevx()

subroutine ssbevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
integer  KD,
real, dimension( ldab, * )  AB,
integer  LDAB,
real, dimension( ldq, * )  Q,
integer  LDQ,
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 
)

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

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

Purpose:
 SSBEVX computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric band 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]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]KD
          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
[in,out]AB
          AB is REAL array, dimension (LDAB, N)
          On entry, the upper or lower triangle of the symmetric band
          matrix A, stored in the first KD+1 rows of the array.  The
          j-th column of A is stored in the j-th column of the array AB
          as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

          On exit, AB is overwritten by values generated during the
          reduction to tridiagonal form.  If UPLO = 'U', the first
          superdiagonal and the diagonal of the tridiagonal matrix T
          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
          the diagonal and first subdiagonal of T are returned in the
          first two rows of AB.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD + 1.
[out]Q
          Q is REAL array, dimension (LDQ, N)
          If JOBZ = 'V', the N-by-N orthogonal matrix used in the
                         reduction to tridiagonal form.
          If JOBZ = 'N', the array Q is not referenced.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.  If JOBZ = 'V', then
          LDQ >= max(1,N).
[in]VL
          VL is REAL
          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 REAL
          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 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 AB 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)
          The first M elements contain 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 (7*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.

Definition at line 262 of file ssbevx.f.

265 *
266 * -- LAPACK driver routine --
267 * -- LAPACK is a software package provided by Univ. of Tennessee, --
268 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269 *
270 * .. Scalar Arguments ..
271  CHARACTER JOBZ, RANGE, UPLO
272  INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
273  REAL ABSTOL, VL, VU
274 * ..
275 * .. Array Arguments ..
276  INTEGER IFAIL( * ), IWORK( * )
277  REAL AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
278  $ Z( LDZ, * )
279 * ..
280 *
281 * =====================================================================
282 *
283 * .. Parameters ..
284  REAL ZERO, ONE
285  parameter( zero = 0.0e0, one = 1.0e0 )
286 * ..
287 * .. Local Scalars ..
288  LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
289  CHARACTER ORDER
290  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
291  $ INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
292  $ NSPLIT
293  REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
294  $ SIGMA, SMLNUM, TMP1, VLL, VUU
295 * ..
296 * .. External Functions ..
297  LOGICAL LSAME
298  REAL SLAMCH, SLANSB
299  EXTERNAL lsame, slamch, slansb
300 * ..
301 * .. External Subroutines ..
302  EXTERNAL scopy, sgemv, slacpy, slascl, ssbtrd, sscal,
304 * ..
305 * .. Intrinsic Functions ..
306  INTRINSIC max, min, sqrt
307 * ..
308 * .. Executable Statements ..
309 *
310 * Test the input parameters.
311 *
312  wantz = lsame( jobz, 'V' )
313  alleig = lsame( range, 'A' )
314  valeig = lsame( range, 'V' )
315  indeig = lsame( range, 'I' )
316  lower = lsame( uplo, 'L' )
317 *
318  info = 0
319  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
320  info = -1
321  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
322  info = -2
323  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
324  info = -3
325  ELSE IF( n.LT.0 ) THEN
326  info = -4
327  ELSE IF( kd.LT.0 ) THEN
328  info = -5
329  ELSE IF( ldab.LT.kd+1 ) THEN
330  info = -7
331  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
332  info = -9
333  ELSE
334  IF( valeig ) THEN
335  IF( n.GT.0 .AND. vu.LE.vl )
336  $ info = -11
337  ELSE IF( indeig ) THEN
338  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
339  info = -12
340  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
341  info = -13
342  END IF
343  END IF
344  END IF
345  IF( info.EQ.0 ) THEN
346  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
347  $ info = -18
348  END IF
349 *
350  IF( info.NE.0 ) THEN
351  CALL xerbla( 'SSBEVX', -info )
352  RETURN
353  END IF
354 *
355 * Quick return if possible
356 *
357  m = 0
358  IF( n.EQ.0 )
359  $ RETURN
360 *
361  IF( n.EQ.1 ) THEN
362  m = 1
363  IF( lower ) THEN
364  tmp1 = ab( 1, 1 )
365  ELSE
366  tmp1 = ab( kd+1, 1 )
367  END IF
368  IF( valeig ) THEN
369  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
370  $ m = 0
371  END IF
372  IF( m.EQ.1 ) THEN
373  w( 1 ) = tmp1
374  IF( wantz )
375  $ z( 1, 1 ) = one
376  END IF
377  RETURN
378  END IF
379 *
380 * Get machine constants.
381 *
382  safmin = slamch( 'Safe minimum' )
383  eps = slamch( 'Precision' )
384  smlnum = safmin / eps
385  bignum = one / smlnum
386  rmin = sqrt( smlnum )
387  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
388 *
389 * Scale matrix to allowable range, if necessary.
390 *
391  iscale = 0
392  abstll = abstol
393  IF ( valeig ) THEN
394  vll = vl
395  vuu = vu
396  ELSE
397  vll = zero
398  vuu = zero
399  ENDIF
400  anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
401  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
402  iscale = 1
403  sigma = rmin / anrm
404  ELSE IF( anrm.GT.rmax ) THEN
405  iscale = 1
406  sigma = rmax / anrm
407  END IF
408  IF( iscale.EQ.1 ) THEN
409  IF( lower ) THEN
410  CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
411  ELSE
412  CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
413  END IF
414  IF( abstol.GT.0 )
415  $ abstll = abstol*sigma
416  IF( valeig ) THEN
417  vll = vl*sigma
418  vuu = vu*sigma
419  END IF
420  END IF
421 *
422 * Call SSBTRD to reduce symmetric band matrix to tridiagonal form.
423 *
424  indd = 1
425  inde = indd + n
426  indwrk = inde + n
427  CALL ssbtrd( jobz, uplo, n, kd, ab, ldab, work( indd ),
428  $ work( inde ), q, ldq, work( indwrk ), iinfo )
429 *
430 * If all eigenvalues are desired and ABSTOL is less than or equal
431 * to zero, then call SSTERF or SSTEQR. If this fails for some
432 * eigenvalue, then try SSTEBZ.
433 *
434  test = .false.
435  IF (indeig) THEN
436  IF (il.EQ.1 .AND. iu.EQ.n) THEN
437  test = .true.
438  END IF
439  END IF
440  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
441  CALL scopy( n, work( indd ), 1, w, 1 )
442  indee = indwrk + 2*n
443  IF( .NOT.wantz ) THEN
444  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
445  CALL ssterf( n, w, work( indee ), info )
446  ELSE
447  CALL slacpy( 'A', n, n, q, ldq, z, ldz )
448  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
449  CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
450  $ work( indwrk ), info )
451  IF( info.EQ.0 ) THEN
452  DO 10 i = 1, n
453  ifail( i ) = 0
454  10 CONTINUE
455  END IF
456  END IF
457  IF( info.EQ.0 ) THEN
458  m = n
459  GO TO 30
460  END IF
461  info = 0
462  END IF
463 *
464 * Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
465 *
466  IF( wantz ) THEN
467  order = 'B'
468  ELSE
469  order = 'E'
470  END IF
471  indibl = 1
472  indisp = indibl + n
473  indiwo = indisp + n
474  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
475  $ work( indd ), work( inde ), m, nsplit, w,
476  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
477  $ iwork( indiwo ), info )
478 *
479  IF( wantz ) THEN
480  CALL sstein( n, work( indd ), work( inde ), m, w,
481  $ iwork( indibl ), iwork( indisp ), z, ldz,
482  $ work( indwrk ), iwork( indiwo ), ifail, info )
483 *
484 * Apply orthogonal matrix used in reduction to tridiagonal
485 * form to eigenvectors returned by SSTEIN.
486 *
487  DO 20 j = 1, m
488  CALL scopy( n, z( 1, j ), 1, work( 1 ), 1 )
489  CALL sgemv( 'N', n, n, one, q, ldq, work, 1, zero,
490  $ z( 1, j ), 1 )
491  20 CONTINUE
492  END IF
493 *
494 * If matrix was scaled, then rescale eigenvalues appropriately.
495 *
496  30 CONTINUE
497  IF( iscale.EQ.1 ) THEN
498  IF( info.EQ.0 ) THEN
499  imax = m
500  ELSE
501  imax = info - 1
502  END IF
503  CALL sscal( imax, one / sigma, w, 1 )
504  END IF
505 *
506 * If eigenvalues are not in order, then sort them, along with
507 * eigenvectors.
508 *
509  IF( wantz ) THEN
510  DO 50 j = 1, m - 1
511  i = 0
512  tmp1 = w( j )
513  DO 40 jj = j + 1, m
514  IF( w( jj ).LT.tmp1 ) THEN
515  i = jj
516  tmp1 = w( jj )
517  END IF
518  40 CONTINUE
519 *
520  IF( i.NE.0 ) THEN
521  itmp1 = iwork( indibl+i-1 )
522  w( i ) = w( j )
523  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
524  w( j ) = tmp1
525  iwork( indibl+j-1 ) = itmp1
526  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
527  IF( info.NE.0 ) THEN
528  itmp1 = ifail( i )
529  ifail( i ) = ifail( j )
530  ifail( j ) = itmp1
531  END IF
532  END IF
533  50 CONTINUE
534  END IF
535 *
536  RETURN
537 *
538 * End of SSBEVX
539 *
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:131
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:273
real function slansb(NORM, UPLO, N, K, AB, LDAB, WORK)
SLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: slansb.f:129
subroutine ssbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
SSBTRD
Definition: ssbtrd.f:163
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
Definition: sstein.f:174
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: