dsyevr.f −
Functions/Subroutines
subroutine dsyevr (JOBZ,
RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z,
LDZ, ISUPPZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSYEVR computes the eigenvalues and, optionally, the left
and/or right eigenvectors for SY matrices
subroutine
dsyevr (character JOBZ, character RANGE, character UPLO,
integer N, double precision, dimension( lda, * ) A, integer
LDA, 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)
Purpose:
DSYEVR computes
selected eigenvalues and, optionally, eigenvectors
of a real symmetric matrix A. Eigenvalues and eigenvectors
can be
selected by specifying either a range of values or a range
of
indices for the desired eigenvalues.
DSYEVR first
reduces the matrix A to tridiagonal form T with a call
to DSYTRD. Then, whenever possible, DSYEVR calls DSTEMR to
compute
the eigenspectrum using Relatively Robust Representations.
DSTEMR
computes eigenvalues by the dqds algorithm, while orthogonal
eigenvectors are computed from various "good" 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 each
unreduced block (submatrix) of T,
(a) Compute T - sigma I = L D L^T, so that L and D
define all the wanted eigenvalues to high relative accuracy.
This means that small relative changes in the entries of D
and L
cause only small relative changes in the eigenvalues and
eigenvectors. The standard (unfactored) representation of
the
tridiagonal matrix T does not have this property in general.
(b) Compute the eigenvalues to suitable accuracy.
If the eigenvectors are desired, the algorithm attains full
accuracy of the computed eigenvalues only right before
the corresponding vectors have to be computed, see steps c)
and d).
(c) For each cluster of close eigenvalues, select a new
shift close to the cluster, find a new factorization, and
refine
the shifted eigenvalues to suitable accuracy.
(d) For each eigenvalue with a large enough relative
separation compute
the corresponding eigenvector by forming a rank revealing
twisted
factorization. Go back to (c) for any clusters that
remain.
The desired
accuracy of the output can be specified by the input
parameter ABSTOL.
For more
details, see DSTEMR’s documentation and:
- Inderjit S. Dhillon and Beresford N. Parlett:
"Multiple representations
to compute orthogonal eigenvectors of symmetric tridiagonal
matrices,"
Linear Algebra and its Applications, 387(1), pp. 1-28,
August 2004.
- Inderjit Dhillon and Beresford Parlett: "Orthogonal
Eigenvectors and
Relative Gaps," SIAM Journal on Matrix Analysis and
Applications, Vol. 25,
2004. Also LAPACK Working Note 154.
- Inderjit Dhillon: "A new O(n^2) algorithm for the
symmetric
tridiagonal eigenvalue/eigenvector problem",
Computer Science Division Technical Report No.
UCB/CSD-97-971,
UC Berkeley, May 1997.
Note 1 : DSYEVR
calls DSTEMR when the full spectrum is requested
on machines which conform to the ieee-754 floating point
standard.
DSYEVR calls DSTEBZ and SSTEIN 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:
JOBZ
JOBZ is
CHARACTER*1
= ’N’: Compute eigenvalues only;
= ’V’: Compute eigenvalues and eigenvectors.
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
UPLO
UPLO is
CHARACTER*1
= ’U’: Upper triangle of A is stored;
= ’L’: Lower triangle of A is stored.
N
N is INTEGER
The order of the matrix A. N >= 0.
A
A is DOUBLE
PRECISION array, dimension (LDA, N)
On entry, the symmetric matrix A. If UPLO = ’U’,
the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO =
’L’,
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO=’L’) or the
upper
triangle (if UPLO=’U’) of A, including the
diagonal, is
destroyed.
LDA
LDA is INTEGER
The leading dimension of the array A. LDA >=
max(1,N).
VL
VL is DOUBLE PRECISION
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’.
IL
IL is INTEGER
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’.
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
"Computing Small Singular Values of Bidiagonal Matrices
with Guaranteed High Relative Accuracy," 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,
"Computing Accurate Eigensystems of Scaled Diagonally
Dominant Matrices", LAPACK Working Note #7, for a
discussion
of which matrices define their eigenvalues to high relative
accuracy.
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.
W
W is DOUBLE
PRECISION array, dimension (N)
The first M elements contain the selected eigenvalues in
ascending order.
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 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.
Supplying N columns is always safe.
LDZ
LDZ is INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = ’V’, LDZ >= max(1,N).
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
WORK
WORK is DOUBLE
PRECISION array, dimension (MAX(1,LWORK))
On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
LWORK
LWORK is
INTEGER
The dimension of the array WORK. LWORK >= max(1,26*N).
For optimal efficiency, LWORK >= (NB+6)*N,
where NB is the max of the blocksize for DSYTRD and DORMTR
returned by ILAENV.
If LWORK = -1,
then a workspace query is assumed; the routine
only calculates the optimal size of the WORK array, returns
this value as the first entry of the WORK array, and no
error
message related to LWORK is issued by XERBLA.
IWORK
IWORK is
INTEGER array, dimension (MAX(1,LIWORK))
On exit, if INFO = 0, IWORK(1) returns the optimal
LWORK.
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 size of the IWORK array,
returns this value as the first entry of the IWORK array,
and
no error message related to LIWORK is issued by XERBLA.
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.
Date:
September 2012
Contributors:
Inderjit Dhillon, IBM Almaden,
USA
Osni Marques, LBNL/NERSC, USA
Ken Stanley, Computer Science Division, University of
California at Berkeley, USA
Jason Riedy, Computer Science Division, University of
California at Berkeley, USA
Definition at line 327 of file dsyevr.f.
327 *
328 * -- LAPACK driver routine (version 3.4.2) --
329 * -- LAPACK is a software package provided by Univ. of
Tennessee, --
330 * -- Univ. of California Berkeley, Univ. of Colorado
Denver and NAG Ltd..--
331 * September 2012
332 *
333 * .. Scalar Arguments ..
334 CHARACTER jobz, range, uplo
335 INTEGER il, info, iu, lda, ldz, liwork, lwork, m, n
336 DOUBLE PRECISION abstol, vl, vu
337 * ..
338 * .. Array Arguments ..
339 INTEGER isuppz( * ), iwork( * )
340 DOUBLE PRECISION a( lda, * ), w( * ), work( * ), z( ldz,
* )
341 * ..
342 *
343 *
=====================================================================
344 *
345 * .. Parameters ..
346 DOUBLE PRECISION zero, one, two
347 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
348 * ..
349 * .. Local Scalars ..
350 LOGICAL alleig, indeig, lower, lquery, valeig, wantz,
351 $ tryrac
352 CHARACTER order
353 INTEGER i, ieeeok, iinfo, imax, indd, inddd, inde,
354 $ indee, indibl, indifl, indisp, indiwo, indtau,
355 $ indwk, indwkn, iscale, j, jj, liwmin,
356 $ llwork, llwrkn, lwkopt, lwmin, nb, nsplit
357 DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin,
safmin,
358 $ sigma, smlnum, tmp1, vll, vuu
359 * ..
360 * .. External Functions ..
361 LOGICAL lsame
362 INTEGER ilaenv
363 DOUBLE PRECISION dlamch, dlansy
364 EXTERNAL lsame, ilaenv, dlamch, dlansy
365 * ..
366 * .. External Subroutines ..
367 EXTERNAL dcopy, dormtr, dscal, dstebz, dstemr, dstein,
368 $ dsterf, dswap, dsytrd, xerbla
369 * ..
370 * .. Intrinsic Functions ..
371 INTRINSIC max, min, sqrt
372 * ..
373 * .. Executable Statements ..
374 *
375 * Test the input parameters.
376 *
377 ieeeok = ilaenv( 10, ’DSYEVR’,
’N’, 1, 2, 3, 4 )
378 *
379 lower = lsame( uplo, ’L’ )
380 wantz = lsame( jobz, ’V’ )
381 alleig = lsame( range, ’A’ )
382 valeig = lsame( range, ’V’ )
383 indeig = lsame( range, ’I’ )
384 *
385 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
386 *
387 lwmin = max( 1, 26*n )
388 liwmin = max( 1, 10*n )
389 *
390 info = 0
391 IF( .NOT.( wantz .OR. lsame( jobz, ’N’ ) ) )
THEN
392 info = -1
393 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
394 info = -2
395 ELSE IF( .NOT.( lower .OR. lsame( uplo, ’U’
) ) ) THEN
396 info = -3
397 ELSE IF( n.LT.0 ) THEN
398 info = -4
399 ELSE IF( lda.LT.max( 1, n ) ) THEN
400 info = -6
401 ELSE
402 IF( valeig ) THEN
403 IF( n.GT.0 .AND. vu.LE.vl )
404 $ info = -8
405 ELSE IF( indeig ) THEN
406 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
407 info = -9
408 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
409 info = -10
410 END IF
411 END IF
412 END IF
413 IF( info.EQ.0 ) THEN
414 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
415 info = -15
416 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
417 info = -18
418 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
419 info = -20
420 END IF
421 END IF
422 *
423 IF( info.EQ.0 ) THEN
424 nb = ilaenv( 1, ’DSYTRD’, uplo, n, -1, -1,
-1 )
425 nb = max( nb, ilaenv( 1, ’DORMTR’, uplo, n,
-1, -1, -1 ) )
426 lwkopt = max( ( nb+1 )*n, lwmin )
427 work( 1 ) = lwkopt
428 iwork( 1 ) = liwmin
429 END IF
430 *
431 IF( info.NE.0 ) THEN
432 CALL xerbla( ’DSYEVR’, -info )
433 RETURN
434 ELSE IF( lquery ) THEN
435 RETURN
436 END IF
437 *
438 * Quick return if possible
439 *
440 m = 0
441 IF( n.EQ.0 ) THEN
442 work( 1 ) = 1
443 RETURN
444 END IF
445 *
446 IF( n.EQ.1 ) THEN
447 work( 1 ) = 7
448 IF( alleig .OR. indeig ) THEN
449 m = 1
450 w( 1 ) = a( 1, 1 )
451 ELSE
452 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
453 m = 1
454 w( 1 ) = a( 1, 1 )
455 END IF
456 END IF
457 IF( wantz ) THEN
458 z( 1, 1 ) = one
459 isuppz( 1 ) = 1
460 isuppz( 2 ) = 1
461 END IF
462 RETURN
463 END IF
464 *
465 * Get machine constants.
466 *
467 safmin = dlamch( ’Safe minimum’ )
468 eps = dlamch( ’Precision’ )
469 smlnum = safmin / eps
470 bignum = one / smlnum
471 rmin = sqrt( smlnum )
472 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) )
)
473 *
474 * Scale matrix to allowable range, if necessary.
475 *
476 iscale = 0
477 abstll = abstol
478 IF (valeig) THEN
479 vll = vl
480 vuu = vu
481 END IF
482 anrm = dlansy( ’M’, uplo, n, a, lda, work )
483 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
484 iscale = 1
485 sigma = rmin / anrm
486 ELSE IF( anrm.GT.rmax ) THEN
487 iscale = 1
488 sigma = rmax / anrm
489 END IF
490 IF( iscale.EQ.1 ) THEN
491 IF( lower ) THEN
492 DO 10 j = 1, n
493 CALL dscal( n-j+1, sigma, a( j, j ), 1 )
494 10 CONTINUE
495 ELSE
496 DO 20 j = 1, n
497 CALL dscal( j, sigma, a( 1, j ), 1 )
498 20 CONTINUE
499 END IF
500 IF( abstol.GT.0 )
501 $ abstll = abstol*sigma
502 IF( valeig ) THEN
503 vll = vl*sigma
504 vuu = vu*sigma
505 END IF
506 END IF
507
508 * Initialize indices into workspaces. Note: The IWORK
indices are
509 * used only if DSTERF or DSTEMR fail.
510
511 * WORK(INDTAU:INDTAU+N-1) stores the scalar factors of
the
512 * elementary reflectors used in DSYTRD.
513 indtau = 1
514 * WORK(INDD:INDD+N-1) stores the tridiagonal’s
diagonal entries.
515 indd = indtau + n
516 * WORK(INDE:INDE+N-1) stores the off-diagonal entries of
the
517 * tridiagonal matrix from DSYTRD.
518 inde = indd + n
519 * WORK(INDDD:INDDD+N-1) is a copy of the diagonal
entries over
520 * -written by DSTEMR (the DSTERF path copies the
diagonal to W).
521 inddd = inde + n
522 * WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal
entries over
523 * -written while computing the eigenvalues in DSTERF and
DSTEMR.
524 indee = inddd + n
525 * INDWK is the starting offset of the left-over
workspace, and
526 * LLWORK is the remaining workspace size.
527 indwk = indee + n
528 llwork = lwork - indwk + 1
529
530 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in
DSTEBZ and
531 * stores the block indices of each of the M<=N
eigenvalues.
532 indibl = 1
533 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in
DSTEBZ and
534 * stores the starting and finishing indices of each
block.
535 indisp = indibl + n
536 * IWORK(INDIFL:INDIFL+N-1) stores the indices of
eigenvectors
537 * that corresponding to eigenvectors that fail to
converge in
538 * DSTEIN. This information is discarded; if any fail,
the driver
539 * returns INFO > 0.
540 indifl = indisp + n
541 * INDIWO is the offset of the remaining integer
workspace.
542 indiwo = indifl + n
543
544 *
545 * Call DSYTRD to reduce symmetric matrix to tridiagonal
form.
546 *
547 CALL dsytrd( uplo, n, a, lda, work( indd ), work( inde
),
548 $ work( indtau ), work( indwk ), llwork, iinfo )
549 *
550 * If all eigenvalues are desired
551 * then call DSTERF or DSTEMR and DORMTR.
552 *
553 IF( ( alleig .OR. ( indeig .AND. il.EQ.1 .AND. iu.EQ.n )
) .AND.
554 $ ieeeok.EQ.1 ) THEN
555 IF( .NOT.wantz ) THEN
556 CALL dcopy( n, work( indd ), 1, w, 1 )
557 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
558 CALL dsterf( n, w, work( indee ), info )
559 ELSE
560 CALL dcopy( n-1, work( inde ), 1, work( indee ), 1 )
561 CALL dcopy( n, work( indd ), 1, work( inddd ), 1 )
562 *
563 IF (abstol .LE. two*n*eps) THEN
564 tryrac = .true.
565 ELSE
566 tryrac = .false.
567 END IF
568 CALL dstemr( jobz, ’A’, n, work( inddd ),
work( indee ),
569 $ vl, vu, il, iu, m, w, z, ldz, n, isuppz,
570 $ tryrac, work( indwk ), lwork, iwork, liwork,
571 $ info )
572 *
573 *
574 *
575 * Apply orthogonal matrix used in reduction to
tridiagonal
576 * form to eigenvectors returned by DSTEIN.
577 *
578 IF( wantz .AND. info.EQ.0 ) THEN
579 indwkn = inde
580 llwrkn = lwork - indwkn + 1
581 CALL dormtr( ’L’, uplo, ’N’, n,
m, a, lda,
582 $ work( indtau ), z, ldz, work( indwkn ),
583 $ llwrkn, iinfo )
584 END IF
585 END IF
586 *
587 *
588 IF( info.EQ.0 ) THEN
589 * Everything worked. Skip DSTEBZ/DSTEIN. IWORK(:) are
590 * undefined.
591 m = n
592 GO TO 30
593 END IF
594 info = 0
595 END IF
596 *
597 * Otherwise, call DSTEBZ and, if eigenvectors are
desired, DSTEIN.
598 * Also call DSTEBZ and DSTEIN if DSTEMR fails.
599 *
600 IF( wantz ) THEN
601 order = ’B’
602 ELSE
603 order = ’E’
604 END IF
605
606 CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
607 $ work( indd ), work( inde ), m, nsplit, w,
608 $ iwork( indibl ), iwork( indisp ), work( indwk ),
609 $ iwork( indiwo ), info )
610 *
611 IF( wantz ) THEN
612 CALL dstein( n, work( indd ), work( inde ), m, w,
613 $ iwork( indibl ), iwork( indisp ), z, ldz,
614 $ work( indwk ), iwork( indiwo ), iwork( indifl ),
615 $ info )
616 *
617 * Apply orthogonal matrix used in reduction to
tridiagonal
618 * form to eigenvectors returned by DSTEIN.
619 *
620 indwkn = inde
621 llwrkn = lwork - indwkn + 1
622 CALL dormtr( ’L’, uplo, ’N’, n,
m, a, lda, work( indtau ), z,
623 $ ldz, work( indwkn ), llwrkn, iinfo )
624 END IF
625 *
626 * If matrix was scaled, then rescale eigenvalues
appropriately.
627 *
628 * Jump here if DSTEMR/DSTEIN succeeded.
629 30 CONTINUE
630 IF( iscale.EQ.1 ) THEN
631 IF( info.EQ.0 ) THEN
632 imax = m
633 ELSE
634 imax = info - 1
635 END IF
636 CALL dscal( imax, one / sigma, w, 1 )
637 END IF
638 *
639 * If eigenvalues are not in order, then sort them, along
with
640 * eigenvectors. Note: We do not sort the IFAIL portion
of IWORK.
641 * It may not be initialized (if DSTEMR/DSTEIN
succeeded), and we do
642 * not return this detailed information to the user.
643 *
644 IF( wantz ) THEN
645 DO 50 j = 1, m - 1
646 i = 0
647 tmp1 = w( j )
648 DO 40 jj = j + 1, m
649 IF( w( jj ).LT.tmp1 ) THEN
650 i = jj
651 tmp1 = w( jj )
652 END IF
653 40 CONTINUE
654 *
655 IF( i.NE.0 ) THEN
656 w( i ) = w( j )
657 w( j ) = tmp1
658 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
659 END IF
660 50 CONTINUE
661 END IF
662 *
663 * Set WORK(1) to optimal workspace size.
664 *
665 work( 1 ) = lwkopt
666 iwork( 1 ) = liwmin
667 *
668 RETURN
669 *
670 * End of DSYEVR
671 *
Generated automatically by Doxygen for LAPACK from the source code.