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

◆ dsbevx()

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
          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 DOUBLE PRECISION
          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 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.

Definition at line 262 of file dsbevx.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 DOUBLE PRECISION ABSTOL, VL, VU
274* ..
275* .. Array Arguments ..
276 INTEGER IFAIL( * ), IWORK( * )
277 DOUBLE PRECISION AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
278 $ Z( LDZ, * )
279* ..
280*
281* =====================================================================
282*
283* .. Parameters ..
284 DOUBLE PRECISION ZERO, ONE
285 parameter( zero = 0.0d0, one = 1.0d0 )
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 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
294 $ SIGMA, SMLNUM, TMP1, VLL, VUU
295* ..
296* .. External Functions ..
297 LOGICAL LSAME
298 DOUBLE PRECISION DLAMCH, DLANSB
299 EXTERNAL lsame, dlamch, dlansb
300* ..
301* .. External Subroutines ..
302 EXTERNAL dcopy, dgemv, dlacpy, dlascl, dsbtrd, dscal,
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( 'DSBEVX', -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 = dlamch( 'Safe minimum' )
383 eps = dlamch( '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 END IF
400 anrm = dlansb( '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 dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
411 ELSE
412 CALL dlascl( '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 DSBTRD to reduce symmetric band matrix to tridiagonal form.
423*
424 indd = 1
425 inde = indd + n
426 indwrk = inde + n
427 CALL dsbtrd( 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 DSTERF or SSTEQR. If this fails for some
432* eigenvalue, then try DSTEBZ.
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 dcopy( n, work( indd ), 1, w, 1 )
442 indee = indwrk + 2*n
443 IF( .NOT.wantz ) THEN
444 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
445 CALL dsterf( n, w, work( indee ), info )
446 ELSE
447 CALL dlacpy( 'A', n, n, q, ldq, z, ldz )
448 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
449 CALL dsteqr( 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 DSTEBZ 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 dstebz( 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 dstein( 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 DSTEIN.
486*
487 DO 20 j = 1, m
488 CALL dcopy( n, z( 1, j ), 1, work( 1 ), 1 )
489 CALL dgemv( '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 dscal( 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 dswap( 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 DSBEVX
539*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dsbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
DSBTRD
Definition dsbtrd.f:163
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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,...
Definition dlansb.f:129
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:143
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
DSTEBZ
Definition dstebz.f:273
subroutine dstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
DSTEIN
Definition dstein.f:174
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:131
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: