LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
subroutine ssbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
SSBTRD
Definition ssbtrd.f:163
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
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
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
subroutine sstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
SSTEIN
Definition sstein.f:174
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 sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: