LAPACK: SRC/dsyev.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
00002 *
00003 *  -- LAPACK driver routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          JOBZ, UPLO
00010       INTEGER            INFO, LDA, LWORK, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DSYEV computes all eigenvalues and, optionally, eigenvectors of a
00020 *  real symmetric matrix A.
00021 *
00022 *  Arguments
00023 *  =========
00024 *
00025 *  JOBZ    (input) CHARACTER*1
00026 *          = 'N':  Compute eigenvalues only;
00027 *          = 'V':  Compute eigenvalues and eigenvectors.
00028 *
00029 *  UPLO    (input) CHARACTER*1
00030 *          = 'U':  Upper triangle of A is stored;
00031 *          = 'L':  Lower triangle of A is stored.
00032 *
00033 *  N       (input) INTEGER
00034 *          The order of the matrix A.  N >= 0.
00035 *
00036 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
00037 *          On entry, the symmetric matrix A.  If UPLO = 'U', the
00038 *          leading N-by-N upper triangular part of A contains the
00039 *          upper triangular part of the matrix A.  If UPLO = 'L',
00040 *          the leading N-by-N lower triangular part of A contains
00041 *          the lower triangular part of the matrix A.
00042 *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
00043 *          orthonormal eigenvectors of the matrix A.
00044 *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
00045 *          or the upper triangle (if UPLO='U') of A, including the
00046 *          diagonal, is destroyed.
00047 *
00048 *  LDA     (input) INTEGER
00049 *          The leading dimension of the array A.  LDA >= max(1,N).
00050 *
00051 *  W       (output) DOUBLE PRECISION array, dimension (N)
00052 *          If INFO = 0, the eigenvalues in ascending order.
00053 *
00054 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00055 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00056 *
00057 *  LWORK   (input) INTEGER
00058 *          The length of the array WORK.  LWORK >= max(1,3*N-1).
00059 *          For optimal efficiency, LWORK >= (NB+2)*N,
00060 *          where NB is the blocksize for DSYTRD returned by ILAENV.
00061 *
00062 *          If LWORK = -1, then a workspace query is assumed; the routine
00063 *          only calculates the optimal size of the WORK array, returns
00064 *          this value as the first entry of the WORK array, and no error
00065 *          message related to LWORK is issued by XERBLA.
00066 *
00067 *  INFO    (output) INTEGER
00068 *          = 0:  successful exit
00069 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00070 *          > 0:  if INFO = i, the algorithm failed to converge; i
00071 *                off-diagonal elements of an intermediate tridiagonal
00072 *                form did not converge to zero.
00073 *
00074 *  =====================================================================
00075 *
00076 *     .. Parameters ..
00077       DOUBLE PRECISION   ZERO, ONE
00078       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00079 *     ..
00080 *     .. Local Scalars ..
00081       LOGICAL            LOWER, LQUERY, WANTZ
00082       INTEGER            IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
00083      $                   LLWORK, LWKOPT, NB
00084       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
00085      $                   SMLNUM
00086 *     ..
00087 *     .. External Functions ..
00088       LOGICAL            LSAME
00089       INTEGER            ILAENV
00090       DOUBLE PRECISION   DLAMCH, DLANSY
00091       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANSY
00092 *     ..
00093 *     .. External Subroutines ..
00094       EXTERNAL           DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD,
00095      $                   XERBLA
00096 *     ..
00097 *     .. Intrinsic Functions ..
00098       INTRINSIC          MAX, SQRT
00099 *     ..
00100 *     .. Executable Statements ..
00101 *
00102 *     Test the input parameters.
00103 *
00104       WANTZ = LSAME( JOBZ, 'V' )
00105       LOWER = LSAME( UPLO, 'L' )
00106       LQUERY = ( LWORK.EQ.-1 )
00107 *
00108       INFO = 0
00109       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00110          INFO = -1
00111       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
00112          INFO = -2
00113       ELSE IF( N.LT.0 ) THEN
00114          INFO = -3
00115       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00116          INFO = -5
00117       END IF
00118 *
00119       IF( INFO.EQ.0 ) THEN
00120          NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
00121          LWKOPT = MAX( 1, ( NB+2 )*N )
00122          WORK( 1 ) = LWKOPT
00123 *
00124          IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
00125      $      INFO = -8
00126       END IF
00127 *
00128       IF( INFO.NE.0 ) THEN
00129          CALL XERBLA( 'DSYEV ', -INFO )
00130          RETURN
00131       ELSE IF( LQUERY ) THEN
00132          RETURN
00133       END IF
00134 *
00135 *     Quick return if possible
00136 *
00137       IF( N.EQ.0 ) THEN
00138          RETURN
00139       END IF
00140 *
00141       IF( N.EQ.1 ) THEN
00142          W( 1 ) = A( 1, 1 )
00143          WORK( 1 ) = 2
00144          IF( WANTZ )
00145      $      A( 1, 1 ) = ONE
00146          RETURN
00147       END IF
00148 *
00149 *     Get machine constants.
00150 *
00151       SAFMIN = DLAMCH( 'Safe minimum' )
00152       EPS = DLAMCH( 'Precision' )
00153       SMLNUM = SAFMIN / EPS
00154       BIGNUM = ONE / SMLNUM
00155       RMIN = SQRT( SMLNUM )
00156       RMAX = SQRT( BIGNUM )
00157 *
00158 *     Scale matrix to allowable range, if necessary.
00159 *
00160       ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
00161       ISCALE = 0
00162       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
00163          ISCALE = 1
00164          SIGMA = RMIN / ANRM
00165       ELSE IF( ANRM.GT.RMAX ) THEN
00166          ISCALE = 1
00167          SIGMA = RMAX / ANRM
00168       END IF
00169       IF( ISCALE.EQ.1 )
00170      $   CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
00171 *
00172 *     Call DSYTRD to reduce symmetric matrix to tridiagonal form.
00173 *
00174       INDE = 1
00175       INDTAU = INDE + N
00176       INDWRK = INDTAU + N
00177       LLWORK = LWORK - INDWRK + 1
00178       CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
00179      $             WORK( INDWRK ), LLWORK, IINFO )
00180 *
00181 *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
00182 *     DORGTR to generate the orthogonal matrix, then call DSTEQR.
00183 *
00184       IF( .NOT.WANTZ ) THEN
00185          CALL DSTERF( N, W, WORK( INDE ), INFO )
00186       ELSE
00187          CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
00188      $                LLWORK, IINFO )
00189          CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
00190      $                INFO )
00191       END IF
00192 *
00193 *     If matrix was scaled, then rescale eigenvalues appropriately.
00194 *
00195       IF( ISCALE.EQ.1 ) THEN
00196          IF( INFO.EQ.0 ) THEN
00197             IMAX = N
00198          ELSE
00199             IMAX = INFO - 1
00200          END IF
00201          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
00202       END IF
00203 *
00204 *     Set WORK(1) to optimal workspace size.
00205 *
00206       WORK( 1 ) = LWKOPT
00207 *
00208       RETURN
00209 *
00210 *     End of DSYEV
00211 *
00212       END