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

◆ dstevr()

subroutine dstevr ( character jobz,
character range,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
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,
integer, dimension( * ) isuppz,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

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

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

Purpose:
!>
!> DSTEVR computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric tridiagonal matrix T.  Eigenvalues and
!> eigenvectors can be selected by specifying either a range of values
!> or a range of indices for the desired eigenvalues.
!>
!> Whenever possible, DSTEVR calls DSTEMR to compute the
!> eigenspectrum using Relatively Robust Representations.  DSTEMR
!> computes eigenvalues by the dqds algorithm, while orthogonal
!> eigenvectors are computed from various  L D L^T representations
!> (also known as Relatively Robust Representations). Gram-Schmidt
!> orthogonalization is avoided as far as possible. More specifically,
!> the various steps of the algorithm are as follows. For the i-th
!> unreduced block of T,
!>    (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
!>         is a relatively robust representation,
!>    (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
!>        relative accuracy by the dqds algorithm,
!>    (c) If there is a cluster of close eigenvalues,  sigma_i
!>        close to the cluster, and go to step (a),
!>    (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
!>        compute the corresponding eigenvector by forming a
!>        rank-revealing twisted factorization.
!> The desired accuracy of the output can be specified by the input
!> parameter ABSTOL.
!>
!> For more details, see , by Inderjit Dhillon,
!> Computer Science Division Technical Report No. UCB//CSD-97-971,
!> UC Berkeley, May 1997.
!>
!>
!> Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of DSTEMR may create NaNs and infinities and
!> hence may abort due to a floating point exception in environments
!> which do not handle NaNs and infinities in the ieee standard default
!> manner.
!>
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.
!>          For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
!>          DSTEIN are called
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix.  N >= 0.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          On entry, the n diagonal elements of the tridiagonal matrix
!>          A.
!>          On exit, D may be multiplied by a constant factor chosen
!>          to avoid over/underflow in computing the eigenvalues.
!>
[in,out]E
!>          E is DOUBLE PRECISION array, dimension (max(1,N-1))
!>          On entry, the (n-1) subdiagonal elements of the tridiagonal
!>          matrix A in elements 1 to N-1 of E.
!>          On exit, E may be multiplied by a constant factor chosen
!>          to avoid over/underflow in computing the eigenvalues.
!>
[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 A to tridiagonal form.
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
!>          If high relative accuracy is important, set ABSTOL to
!>          DLAMCH( 'Safe minimum' ).  Doing so will guarantee that
!>          eigenvalues are computed to high relative accuracy when
!>          possible in future releases.  The current code does not
!>          make any guarantees about high relative accuracy, but
!>          future releases will. See J. Barlow and J. Demmel,
!>          , LAPACK Working Note #7, for a discussion
!>          of which matrices define their eigenvalues to high relative
!>          accuracy.
!>
[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).
!>          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]ISUPPZ
!>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
!>          The support of the eigenvectors in Z, i.e., the indices
!>          indicating the nonzero elements in Z. The i-th eigenvector
!>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
!>          ISUPPZ( 2*i ).
!>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal (and
!>          minimal) LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,20*N).
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal sizes of the WORK and IWORK
!>          arrays, returns these values as the first entries of the WORK
!>          and IWORK arrays, and no error message related to LWORK or
!>          LIWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal (and
!>          minimal) LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.  LIWORK >= max(1,10*N).
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal sizes of the WORK and
!>          IWORK arrays, returns these values as the first entries of
!>          the WORK and IWORK arrays, and no error message related to
!>          LWORK or LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  Internal error
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Inderjit Dhillon, IBM Almaden, USA
Osni Marques, LBNL/NERSC, USA
Ken Stanley, Computer Science Division, University of California at Berkeley, USA

Definition at line 299 of file dstevr.f.

303*
304* -- LAPACK driver routine --
305* -- LAPACK is a software package provided by Univ. of Tennessee, --
306* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
307*
308* .. Scalar Arguments ..
309 CHARACTER JOBZ, RANGE
310 INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
311 DOUBLE PRECISION ABSTOL, VL, VU
312* ..
313* .. Array Arguments ..
314 INTEGER ISUPPZ( * ), IWORK( * )
315 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
316* ..
317*
318* =====================================================================
319*
320* .. Parameters ..
321 DOUBLE PRECISION ZERO, ONE, TWO
322 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
323* ..
324* .. Local Scalars ..
325 LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
326 $ TRYRAC
327 CHARACTER ORDER
328 INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
329 $ INDIWO, ISCALE, ITMP1, J, JJ, LIWMIN, LWMIN,
330 $ NSPLIT
331 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
332 $ TMP1, TNRM, VLL, VUU
333* ..
334* .. External Functions ..
335 LOGICAL LSAME
336 INTEGER ILAENV
337 DOUBLE PRECISION DLAMCH, DLANST
338 EXTERNAL lsame, ilaenv, dlamch, dlanst
339* ..
340* .. External Subroutines ..
341 EXTERNAL dcopy, dscal, dstebz, dstemr, dstein,
342 $ dsterf,
343 $ dswap, xerbla
344* ..
345* .. Intrinsic Functions ..
346 INTRINSIC max, min, sqrt
347* ..
348* .. Executable Statements ..
349*
350*
351* Test the input parameters.
352*
353 ieeeok = ilaenv( 10, 'DSTEVR', 'N', 1, 2, 3, 4 )
354*
355 wantz = lsame( jobz, 'V' )
356 alleig = lsame( range, 'A' )
357 valeig = lsame( range, 'V' )
358 indeig = lsame( range, 'I' )
359*
360 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
361 lwmin = max( 1, 20*n )
362 liwmin = max( 1, 10*n )
363*
364*
365 info = 0
366 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
367 info = -1
368 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
369 info = -2
370 ELSE IF( n.LT.0 ) THEN
371 info = -3
372 ELSE
373 IF( valeig ) THEN
374 IF( n.GT.0 .AND. vu.LE.vl )
375 $ info = -7
376 ELSE IF( indeig ) THEN
377 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
378 info = -8
379 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
380 info = -9
381 END IF
382 END IF
383 END IF
384 IF( info.EQ.0 ) THEN
385 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
386 info = -14
387 END IF
388 END IF
389*
390 IF( info.EQ.0 ) THEN
391 work( 1 ) = lwmin
392 iwork( 1 ) = liwmin
393*
394 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
395 info = -17
396 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
397 info = -19
398 END IF
399 END IF
400*
401 IF( info.NE.0 ) THEN
402 CALL xerbla( 'DSTEVR', -info )
403 RETURN
404 ELSE IF( lquery ) THEN
405 RETURN
406 END IF
407*
408* Quick return if possible
409*
410 m = 0
411 IF( n.EQ.0 )
412 $ RETURN
413*
414 IF( n.EQ.1 ) THEN
415 IF( alleig .OR. indeig ) THEN
416 m = 1
417 w( 1 ) = d( 1 )
418 ELSE
419 IF( vl.LT.d( 1 ) .AND. vu.GE.d( 1 ) ) THEN
420 m = 1
421 w( 1 ) = d( 1 )
422 END IF
423 END IF
424 IF( wantz )
425 $ z( 1, 1 ) = one
426 RETURN
427 END IF
428*
429* Get machine constants.
430*
431 safmin = dlamch( 'Safe minimum' )
432 eps = dlamch( 'Precision' )
433 smlnum = safmin / eps
434 bignum = one / smlnum
435 rmin = sqrt( smlnum )
436 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
437*
438*
439* Scale matrix to allowable range, if necessary.
440*
441 iscale = 0
442 IF( valeig ) THEN
443 vll = vl
444 vuu = vu
445 END IF
446*
447 tnrm = dlanst( 'M', n, d, e )
448 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
449 iscale = 1
450 sigma = rmin / tnrm
451 ELSE IF( tnrm.GT.rmax ) THEN
452 iscale = 1
453 sigma = rmax / tnrm
454 END IF
455 IF( iscale.EQ.1 ) THEN
456 CALL dscal( n, sigma, d, 1 )
457 CALL dscal( n-1, sigma, e( 1 ), 1 )
458 IF( valeig ) THEN
459 vll = vl*sigma
460 vuu = vu*sigma
461 END IF
462 END IF
463
464* Initialize indices into workspaces. Note: These indices are used only
465* if DSTERF or DSTEMR fail.
466
467* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
468* stores the block indices of each of the M<=N eigenvalues.
469 indibl = 1
470* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
471* stores the starting and finishing indices of each block.
472 indisp = indibl + n
473* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
474* that corresponding to eigenvectors that fail to converge in
475* DSTEIN. This information is discarded; if any fail, the driver
476* returns INFO > 0.
477 indifl = indisp + n
478* INDIWO is the offset of the remaining integer workspace.
479 indiwo = indisp + n
480*
481* If all eigenvalues are desired, then
482* call DSTERF or DSTEMR. If this fails for some eigenvalue, then
483* try DSTEBZ.
484*
485*
486 test = .false.
487 IF( indeig ) THEN
488 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
489 test = .true.
490 END IF
491 END IF
492 IF( ( alleig .OR. test ) .AND. ieeeok.EQ.1 ) THEN
493 CALL dcopy( n-1, e( 1 ), 1, work( 1 ), 1 )
494 IF( .NOT.wantz ) THEN
495 CALL dcopy( n, d, 1, w, 1 )
496 CALL dsterf( n, w, work, info )
497 ELSE
498 CALL dcopy( n, d, 1, work( n+1 ), 1 )
499 IF (abstol .LE. two*n*eps) THEN
500 tryrac = .true.
501 ELSE
502 tryrac = .false.
503 END IF
504 CALL dstemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,
505 $ iu, m, w, z, ldz, n, isuppz, tryrac,
506 $ work( 2*n+1 ), lwork-2*n, iwork, liwork, info )
507*
508 END IF
509 IF( info.EQ.0 ) THEN
510 m = n
511 GO TO 10
512 END IF
513 info = 0
514 END IF
515*
516* Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
517*
518 IF( wantz ) THEN
519 order = 'B'
520 ELSE
521 order = 'E'
522 END IF
523
524 CALL dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e,
525 $ m,
526 $ nsplit, w, iwork( indibl ), iwork( indisp ), work,
527 $ iwork( indiwo ), info )
528*
529 IF( wantz ) THEN
530 CALL dstein( n, d, e, m, w, iwork( indibl ),
531 $ iwork( indisp ),
532 $ z, ldz, work, iwork( indiwo ), iwork( indifl ),
533 $ info )
534 END IF
535*
536* If matrix was scaled, then rescale eigenvalues appropriately.
537*
538 10 CONTINUE
539 IF( iscale.EQ.1 ) THEN
540 IF( info.EQ.0 ) THEN
541 imax = m
542 ELSE
543 imax = info - 1
544 END IF
545 CALL dscal( imax, one / sigma, w, 1 )
546 END IF
547*
548* If eigenvalues are not in order, then sort them, along with
549* eigenvectors.
550*
551 IF( wantz ) THEN
552 DO 30 j = 1, m - 1
553 i = 0
554 tmp1 = w( j )
555 DO 20 jj = j + 1, m
556 IF( w( jj ).LT.tmp1 ) THEN
557 i = jj
558 tmp1 = w( jj )
559 END IF
560 20 CONTINUE
561*
562 IF( i.NE.0 ) THEN
563 itmp1 = iwork( i )
564 w( i ) = w( j )
565 iwork( i ) = iwork( j )
566 w( j ) = tmp1
567 iwork( j ) = itmp1
568 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
569 END IF
570 30 CONTINUE
571 END IF
572*
573* Causes problems with tests 19 & 20:
574* IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002
575*
576*
577 work( 1 ) = lwmin
578 iwork( 1 ) = liwmin
579 RETURN
580*
581* End of DSTEVR
582*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanst
double precision function dlanst(norm, n, d, e)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlanst.f:98
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dstebz
subroutine dstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
DSTEBZ
Definition dstebz.f:272
dstein
subroutine dstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
DSTEIN
Definition dstein.f:172
dstemr
subroutine dstemr(jobz, range, n, d, e, vl, vu, il, iu, m, w, z, ldz, nzc, isuppz, tryrac, work, lwork, iwork, liwork, info)
DSTEMR
Definition dstemr.f:320
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:84
dswap
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: