src/LAPACK/dsyev.f

   1       SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
   2 *
   3 *  -- LAPACK driver routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       CHARACTER          JOBZ, UPLO
   9       INTEGER            INFO, LDA, LWORK, N
  10 *     ..
  11 *     .. Array Arguments ..
  12       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
  13 *     ..
  14 *
  15 *  Purpose
  16 *  =======
  17 *
  18 *  DSYEV computes all eigenvalues and, optionally, eigenvectors of a
  19 *  real symmetric matrix A.
  20 *
  21 *  Arguments
  22 *  =========
  23 *
  24 *  JOBZ    (input) CHARACTER*1
  25 *          = 'N':  Compute eigenvalues only;
  26 *          = 'V':  Compute eigenvalues and eigenvectors.
  27 *
  28 *  UPLO    (input) CHARACTER*1
  29 *          = 'U':  Upper triangle of A is stored;
  30 *          = 'L':  Lower triangle of A is stored.
  31 *
  32 *  N       (input) INTEGER
  33 *          The order of the matrix A.  N >= 0.
  34 *
  35 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
  36 *          On entry, the symmetric matrix A.  If UPLO = 'U', the
  37 *          leading N-by-N upper triangular part of A contains the
  38 *          upper triangular part of the matrix A.  If UPLO = 'L',
  39 *          the leading N-by-N lower triangular part of A contains
  40 *          the lower triangular part of the matrix A.
  41 *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  42 *          orthonormal eigenvectors of the matrix A.
  43 *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
  44 *          or the upper triangle (if UPLO='U') of A, including the
  45 *          diagonal, is destroyed.
  46 *
  47 *  LDA     (input) INTEGER
  48 *          The leading dimension of the array A.  LDA >= max(1,N).
  49 *
  50 *  W       (output) DOUBLE PRECISION array, dimension (N)
  51 *          If INFO = 0, the eigenvalues in ascending order.
  52 *
  53 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  54 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  55 *
  56 *  LWORK   (input) INTEGER
  57 *          The length of the array WORK.  LWORK >= max(1,3*N-1).
  58 *          For optimal efficiency, LWORK >= (NB+2)*N,
  59 *          where NB is the blocksize for DSYTRD returned by ILAENV.
  60 *
  61 *          If LWORK = -1, then a workspace query is assumed; the routine
  62 *          only calculates the optimal size of the WORK array, returns
  63 *          this value as the first entry of the WORK array, and no error
  64 *          message related to LWORK is issued by XERBLA.
  65 *
  66 *  INFO    (output) INTEGER
  67 *          = 0:  successful exit
  68 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  69 *          > 0:  if INFO = i, the algorithm failed to converge; i
  70 *                off-diagonal elements of an intermediate tridiagonal
  71 *                form did not converge to zero.
  72 *
  73 *  =====================================================================
  74 *
  75 *     .. Parameters ..
  76       DOUBLE PRECISION   ZERO, ONE
  77       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  78 *     ..
  79 *     .. Local Scalars ..
  80       LOGICAL            LOWER, LQUERY, WANTZ
  81       INTEGER            IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
  82      $                   LLWORK, LWKOPT, NB
  83       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
  84      $                   SMLNUM
  85 *     ..
  86 *     .. External Functions ..
  87       LOGICAL            LSAME
  88       INTEGER            ILAENV
  89       DOUBLE PRECISION   DLAMCH, DLANSY
  90       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANSY
  91 *     ..
  92 *     .. External Subroutines ..
  93       EXTERNAL           DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD,
  94      $                   XERBLA
  95 *     ..
  96 *     .. Intrinsic Functions ..
  97       INTRINSIC          MAX, SQRT
  98 *     ..
  99 *     .. Executable Statements ..
 100 *
 101 *     Test the input parameters.
 102 *
 103       WANTZ = LSAME( JOBZ, 'V' )
 104       LOWER = LSAME( UPLO, 'L' )
 105       LQUERY = ( LWORK.EQ.-1 )
 106 *
 107       INFO = 0
 108       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 109          INFO = -1
 110       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
 111          INFO = -2
 112       ELSE IF( N.LT.0 ) THEN
 113          INFO = -3
 114       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 115          INFO = -5
 116       END IF
 117 *
 118       IF( INFO.EQ.0 ) THEN
 119          NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
 120          LWKOPT = MAX( 1, ( NB+2 )*N )
 121          WORK( 1 ) = LWKOPT
 122 *
 123          IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
 124      $      INFO = -8
 125       END IF
 126 *
 127       IF( INFO.NE.0 ) THEN
 128          CALL XERBLA( 'DSYEV ', -INFO )
 129          RETURN
 130       ELSE IF( LQUERY ) THEN
 131          RETURN
 132       END IF
 133 *
 134 *     Quick return if possible
 135 *
 136       IF( N.EQ.0 ) THEN
 137          RETURN
 138       END IF
 139 *
 140       IF( N.EQ.1 ) THEN
 141          W( 1 ) = A( 1, 1 )
 142          WORK( 1 ) = 2
 143          IF( WANTZ )
 144      $      A( 1, 1 ) = ONE
 145          RETURN
 146       END IF
 147 *
 148 *     Get machine constants.
 149 *
 150       SAFMIN = DLAMCH( 'Safe minimum' )
 151       EPS = DLAMCH( 'Precision' )
 152       SMLNUM = SAFMIN / EPS
 153       BIGNUM = ONE / SMLNUM
 154       RMIN = SQRT( SMLNUM )
 155       RMAX = SQRT( BIGNUM )
 156 *
 157 *     Scale matrix to allowable range, if necessary.
 158 *
 159       ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
 160       ISCALE = 0
 161       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
 162          ISCALE = 1
 163          SIGMA = RMIN / ANRM
 164       ELSE IF( ANRM.GT.RMAX ) THEN
 165          ISCALE = 1
 166          SIGMA = RMAX / ANRM
 167       END IF
 168       IF( ISCALE.EQ.1 )
 169      $   CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
 170 *
 171 *     Call DSYTRD to reduce symmetric matrix to tridiagonal form.
 172 *
 173       INDE = 1
 174       INDTAU = INDE + N
 175       INDWRK = INDTAU + N
 176       LLWORK = LWORK - INDWRK + 1
 177       CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
 178      $             WORK( INDWRK ), LLWORK, IINFO )
 179 *
 180 *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
 181 *     DORGTR to generate the orthogonal matrix, then call DSTEQR.
 182 *
 183       IF( .NOT.WANTZ ) THEN
 184          CALL DSTERF( N, W, WORK( INDE ), INFO )
 185       ELSE
 186          CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
 187      $                LLWORK, IINFO )
 188          CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
 189      $                INFO )
 190       END IF
 191 *
 192 *     If matrix was scaled, then rescale eigenvalues appropriately.
 193 *
 194       IF( ISCALE.EQ.1 ) THEN
 195          IF( INFO.EQ.0 ) THEN
 196             IMAX = N
 197          ELSE
 198             IMAX = INFO - 1
 199          END IF
 200          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
 201       END IF
 202 *
 203 *     Set WORK(1) to optimal workspace size.
 204 *
 205       WORK( 1 ) = LWKOPT
 206 *
 207       RETURN
 208 *
 209 *     End of DSYEV
 210 *
 211       END