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

Go to the source code of this file.

Functions/Subroutines

subroutine dsbevx (JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO)
  DSBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices More...
 

Function/Subroutine Documentation

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

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

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

Purpose:
 DSBEVX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 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 AB 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)
          The first M elements contain 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 (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.
Date
November 2011

Definition at line 260 of file dsbevx.f.

260 *
261 * -- LAPACK driver routine (version 3.4.0) --
262 * -- LAPACK is a software package provided by Univ. of Tennessee, --
263 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
264 * November 2011
265 *
266 * .. Scalar Arguments ..
267  CHARACTER jobz, range, uplo
268  INTEGER il, info, iu, kd, ldab, ldq, ldz, m, n
269  DOUBLE PRECISION abstol, vl, vu
270 * ..
271 * .. Array Arguments ..
272  INTEGER ifail( * ), iwork( * )
273  DOUBLE PRECISION ab( ldab, * ), q( ldq, * ), w( * ), work( * ),
274  $ z( ldz, * )
275 * ..
276 *
277 * =====================================================================
278 *
279 * .. Parameters ..
280  DOUBLE PRECISION zero, one
281  parameter( zero = 0.0d0, one = 1.0d0 )
282 * ..
283 * .. Local Scalars ..
284  LOGICAL alleig, indeig, lower, test, valeig, wantz
285  CHARACTER order
286  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
287  $ indisp, indiwo, indwrk, iscale, itmp1, j, jj,
288  $ nsplit
289  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
290  $ sigma, smlnum, tmp1, vll, vuu
291 * ..
292 * .. External Functions ..
293  LOGICAL lsame
294  DOUBLE PRECISION dlamch, dlansb
295  EXTERNAL lsame, dlamch, dlansb
296 * ..
297 * .. External Subroutines ..
298  EXTERNAL dcopy, dgemv, dlacpy, dlascl, dsbtrd, dscal,
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC max, min, sqrt
303 * ..
304 * .. Executable Statements ..
305 *
306 * Test the input parameters.
307 *
308  wantz = lsame( jobz, 'V' )
309  alleig = lsame( range, 'A' )
310  valeig = lsame( range, 'V' )
311  indeig = lsame( range, 'I' )
312  lower = lsame( uplo, 'L' )
313 *
314  info = 0
315  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
316  info = -1
317  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
318  info = -2
319  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
320  info = -3
321  ELSE IF( n.LT.0 ) THEN
322  info = -4
323  ELSE IF( kd.LT.0 ) THEN
324  info = -5
325  ELSE IF( ldab.LT.kd+1 ) THEN
326  info = -7
327  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
328  info = -9
329  ELSE
330  IF( valeig ) THEN
331  IF( n.GT.0 .AND. vu.LE.vl )
332  $ info = -11
333  ELSE IF( indeig ) THEN
334  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
335  info = -12
336  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
337  info = -13
338  END IF
339  END IF
340  END IF
341  IF( info.EQ.0 ) THEN
342  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
343  $ info = -18
344  END IF
345 *
346  IF( info.NE.0 ) THEN
347  CALL xerbla( 'DSBEVX', -info )
348  RETURN
349  END IF
350 *
351 * Quick return if possible
352 *
353  m = 0
354  IF( n.EQ.0 )
355  $ RETURN
356 *
357  IF( n.EQ.1 ) THEN
358  m = 1
359  IF( lower ) THEN
360  tmp1 = ab( 1, 1 )
361  ELSE
362  tmp1 = ab( kd+1, 1 )
363  END IF
364  IF( valeig ) THEN
365  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
366  $ m = 0
367  END IF
368  IF( m.EQ.1 ) THEN
369  w( 1 ) = tmp1
370  IF( wantz )
371  $ z( 1, 1 ) = one
372  END IF
373  RETURN
374  END IF
375 *
376 * Get machine constants.
377 *
378  safmin = dlamch( 'Safe minimum' )
379  eps = dlamch( 'Precision' )
380  smlnum = safmin / eps
381  bignum = one / smlnum
382  rmin = sqrt( smlnum )
383  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
384 *
385 * Scale matrix to allowable range, if necessary.
386 *
387  iscale = 0
388  abstll = abstol
389  IF( valeig ) THEN
390  vll = vl
391  vuu = vu
392  ELSE
393  vll = zero
394  vuu = zero
395  END IF
396  anrm = dlansb( 'M', uplo, n, kd, ab, ldab, work )
397  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
398  iscale = 1
399  sigma = rmin / anrm
400  ELSE IF( anrm.GT.rmax ) THEN
401  iscale = 1
402  sigma = rmax / anrm
403  END IF
404  IF( iscale.EQ.1 ) THEN
405  IF( lower ) THEN
406  CALL dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
407  ELSE
408  CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
409  END IF
410  IF( abstol.GT.0 )
411  $ abstll = abstol*sigma
412  IF( valeig ) THEN
413  vll = vl*sigma
414  vuu = vu*sigma
415  END IF
416  END IF
417 *
418 * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
419 *
420  indd = 1
421  inde = indd + n
422  indwrk = inde + n
423  CALL dsbtrd( jobz, uplo, n, kd, ab, ldab, work( indd ),
424  $ work( inde ), q, ldq, work( indwrk ), iinfo )
425 *
426 * If all eigenvalues are desired and ABSTOL is less than or equal
427 * to zero, then call DSTERF or SSTEQR. If this fails for some
428 * eigenvalue, then try DSTEBZ.
429 *
430  test = .false.
431  IF (indeig) THEN
432  IF (il.EQ.1 .AND. iu.EQ.n) THEN
433  test = .true.
434  END IF
435  END IF
436  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
437  CALL dcopy( n, work( indd ), 1, w, 1 )
438  indee = indwrk + 2*n
439  IF( .NOT.wantz ) THEN
440  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
441  CALL dsterf( n, w, work( indee ), info )
442  ELSE
443  CALL dlacpy( 'A', n, n, q, ldq, z, ldz )
444  CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
445  CALL dsteqr( jobz, n, w, work( indee ), z, ldz,
446  $ work( indwrk ), info )
447  IF( info.EQ.0 ) THEN
448  DO 10 i = 1, n
449  ifail( i ) = 0
450  10 CONTINUE
451  END IF
452  END IF
453  IF( info.EQ.0 ) THEN
454  m = n
455  GO TO 30
456  END IF
457  info = 0
458  END IF
459 *
460 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
461 *
462  IF( wantz ) THEN
463  order = 'B'
464  ELSE
465  order = 'E'
466  END IF
467  indibl = 1
468  indisp = indibl + n
469  indiwo = indisp + n
470  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
471  $ work( indd ), work( inde ), m, nsplit, w,
472  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
473  $ iwork( indiwo ), info )
474 *
475  IF( wantz ) THEN
476  CALL dstein( n, work( indd ), work( inde ), m, w,
477  $ iwork( indibl ), iwork( indisp ), z, ldz,
478  $ work( indwrk ), iwork( indiwo ), ifail, info )
479 *
480 * Apply orthogonal matrix used in reduction to tridiagonal
481 * form to eigenvectors returned by DSTEIN.
482 *
483  DO 20 j = 1, m
484  CALL dcopy( n, z( 1, j ), 1, work( 1 ), 1 )
485  CALL dgemv( 'N', n, n, one, q, ldq, work, 1, zero,
486  $ z( 1, j ), 1 )
487  20 CONTINUE
488  END IF
489 *
490 * If matrix was scaled, then rescale eigenvalues appropriately.
491 *
492  30 CONTINUE
493  IF( iscale.EQ.1 ) THEN
494  IF( info.EQ.0 ) THEN
495  imax = m
496  ELSE
497  imax = info - 1
498  END IF
499  CALL dscal( imax, one / sigma, w, 1 )
500  END IF
501 *
502 * If eigenvalues are not in order, then sort them, along with
503 * eigenvectors.
504 *
505  IF( wantz ) THEN
506  DO 50 j = 1, m - 1
507  i = 0
508  tmp1 = w( j )
509  DO 40 jj = j + 1, m
510  IF( w( jj ).LT.tmp1 ) THEN
511  i = jj
512  tmp1 = w( jj )
513  END IF
514  40 CONTINUE
515 *
516  IF( i.NE.0 ) THEN
517  itmp1 = iwork( indibl+i-1 )
518  w( i ) = w( j )
519  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
520  w( j ) = tmp1
521  iwork( indibl+j-1 ) = itmp1
522  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
523  IF( info.NE.0 ) THEN
524  itmp1 = ifail( i )
525  ifail( i ) = ifail( j )
526  ifail( j ) = itmp1
527  END IF
528  END IF
529  50 CONTINUE
530  END IF
531 *
532  RETURN
533 *
534 * End of DSBEVX
535 *
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 dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine dsbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
DSBTRD
Definition: dsbtrd.f:165
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:141
double precision function dlansb(NORM, UPLO, N, K, AB, LDAB, WORK)
DLANSB 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 band matrix.
Definition: dlansb.f:131
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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
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
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158

Here is the call graph for this function:

Here is the caller graph for this function: