SRC/zheevd.f

   1 *> \brief <b> ZHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHEEVD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheevd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheevd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheevd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
  22 *                          LRWORK, IWORK, LIWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ, UPLO
  26 *       INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IWORK( * )
  30 *       DOUBLE PRECISION   RWORK( * ), W( * )
  31 *       COMPLEX*16         A( LDA, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZHEEVD computes all eigenvalues and, optionally, eigenvectors of a
  41 *> complex Hermitian matrix A.  If eigenvectors are desired, it uses a
  42 *> divide and conquer algorithm.
  43 *>
  44 *> The divide and conquer algorithm makes very mild assumptions about
  45 *> floating point arithmetic. It will work on machines with a guard
  46 *> digit in add/subtract, or on those binary machines without guard
  47 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
  48 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
  49 *> without guard digits, but we know of none.
  50 *> \endverbatim
  51 *
  52 *  Arguments:
  53 *  ==========
  54 *
  55 *> \param[in] JOBZ
  56 *> \verbatim
  57 *>          JOBZ is CHARACTER*1
  58 *>          = 'N':  Compute eigenvalues only;
  59 *>          = 'V':  Compute eigenvalues and eigenvectors.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] UPLO
  63 *> \verbatim
  64 *>          UPLO is CHARACTER*1
  65 *>          = 'U':  Upper triangle of A is stored;
  66 *>          = 'L':  Lower triangle of A is stored.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] N
  70 *> \verbatim
  71 *>          N is INTEGER
  72 *>          The order of the matrix A.  N >= 0.
  73 *> \endverbatim
  74 *>
  75 *> \param[in,out] A
  76 *> \verbatim
  77 *>          A is COMPLEX*16 array, dimension (LDA, N)
  78 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
  79 *>          leading N-by-N upper triangular part of A contains the
  80 *>          upper triangular part of the matrix A.  If UPLO = 'L',
  81 *>          the leading N-by-N lower triangular part of A contains
  82 *>          the lower triangular part of the matrix A.
  83 *>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  84 *>          orthonormal eigenvectors of the matrix A.
  85 *>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
  86 *>          or the upper triangle (if UPLO='U') of A, including the
  87 *>          diagonal, is destroyed.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] LDA
  91 *> \verbatim
  92 *>          LDA is INTEGER
  93 *>          The leading dimension of the array A.  LDA >= max(1,N).
  94 *> \endverbatim
  95 *>
  96 *> \param[out] W
  97 *> \verbatim
  98 *>          W is DOUBLE PRECISION array, dimension (N)
  99 *>          If INFO = 0, the eigenvalues in ascending order.
 100 *> \endverbatim
 101 *>
 102 *> \param[out] WORK
 103 *> \verbatim
 104 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 105 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] LWORK
 109 *> \verbatim
 110 *>          LWORK is INTEGER
 111 *>          The length of the array WORK.
 112 *>          If N <= 1,                LWORK must be at least 1.
 113 *>          If JOBZ  = 'N' and N > 1, LWORK must be at least N + 1.
 114 *>          If JOBZ  = 'V' and N > 1, LWORK must be at least 2*N + N**2.
 115 *>
 116 *>          If LWORK = -1, then a workspace query is assumed; the routine
 117 *>          only calculates the optimal sizes of the WORK, RWORK and
 118 *>          IWORK arrays, returns these values as the first entries of
 119 *>          the WORK, RWORK and IWORK arrays, and no error message
 120 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 121 *> \endverbatim
 122 *>
 123 *> \param[out] RWORK
 124 *> \verbatim
 125 *>          RWORK is DOUBLE PRECISION array,
 126 *>                                         dimension (LRWORK)
 127 *>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
 128 *> \endverbatim
 129 *>
 130 *> \param[in] LRWORK
 131 *> \verbatim
 132 *>          LRWORK is INTEGER
 133 *>          The dimension of the array RWORK.
 134 *>          If N <= 1,                LRWORK must be at least 1.
 135 *>          If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
 136 *>          If JOBZ  = 'V' and N > 1, LRWORK must be at least
 137 *>                         1 + 5*N + 2*N**2.
 138 *>
 139 *>          If LRWORK = -1, then a workspace query is assumed; the
 140 *>          routine only calculates the optimal sizes of the WORK, RWORK
 141 *>          and IWORK arrays, returns these values as the first entries
 142 *>          of the WORK, RWORK and IWORK arrays, and no error message
 143 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 144 *> \endverbatim
 145 *>
 146 *> \param[out] IWORK
 147 *> \verbatim
 148 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 149 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 150 *> \endverbatim
 151 *>
 152 *> \param[in] LIWORK
 153 *> \verbatim
 154 *>          LIWORK is INTEGER
 155 *>          The dimension of the array IWORK.
 156 *>          If N <= 1,                LIWORK must be at least 1.
 157 *>          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
 158 *>          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
 159 *>
 160 *>          If LIWORK = -1, then a workspace query is assumed; the
 161 *>          routine only calculates the optimal sizes of the WORK, RWORK
 162 *>          and IWORK arrays, returns these values as the first entries
 163 *>          of the WORK, RWORK and IWORK arrays, and no error message
 164 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 165 *> \endverbatim
 166 *>
 167 *> \param[out] INFO
 168 *> \verbatim
 169 *>          INFO is INTEGER
 170 *>          = 0:  successful exit
 171 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 172 *>          > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
 173 *>                to converge; i off-diagonal elements of an intermediate
 174 *>                tridiagonal form did not converge to zero;
 175 *>                if INFO = i and JOBZ = 'V', then the algorithm failed
 176 *>                to compute an eigenvalue while working on the submatrix
 177 *>                lying in rows and columns INFO/(N+1) through
 178 *>                mod(INFO,N+1).
 179 *> \endverbatim
 180 *
 181 *  Authors:
 182 *  ========
 183 *
 184 *> \author Univ. of Tennessee
 185 *> \author Univ. of California Berkeley
 186 *> \author Univ. of Colorado Denver
 187 *> \author NAG Ltd.
 188 *
 189 *> \date November 2011
 190 *
 191 *> \ingroup complex16HEeigen
 192 *
 193 *> \par Further Details:
 194 *  =====================
 195 *>
 196 *>  Modified description of INFO. Sven, 16 Feb 05.
 197 *
 198 *> \par Contributors:
 199 *  ==================
 200 *>
 201 *> Jeff Rutter, Computer Science Division, University of California
 202 *> at Berkeley, USA
 203 *>
 204 *  =====================================================================
 205       SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
 206      $                   LRWORK, IWORK, LIWORK, INFO )
 207 *
 208 *  -- LAPACK driver routine (version 3.4.0) --
 209 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 210 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 211 *     November 2011
 212 *
 213 *     .. Scalar Arguments ..
 214       CHARACTER          JOBZ, UPLO
 215       INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
 216 *     ..
 217 *     .. Array Arguments ..
 218       INTEGER            IWORK( * )
 219       DOUBLE PRECISION   RWORK( * ), W( * )
 220       COMPLEX*16         A( LDA, * ), WORK( * )
 221 *     ..
 222 *
 223 *  =====================================================================
 224 *
 225 *     .. Parameters ..
 226       DOUBLE PRECISION   ZERO, ONE
 227       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 228       COMPLEX*16         CONE
 229       PARAMETER          ( CONE = ( 1.0D0, 0.0D0 ) )
 230 *     ..
 231 *     .. Local Scalars ..
 232       LOGICAL            LOWER, LQUERY, WANTZ
 233       INTEGER            IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
 234      $                   INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK,
 235      $                   LLWRK2, LOPT, LROPT, LRWMIN, LWMIN
 236       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
 237      $                   SMLNUM
 238 *     ..
 239 *     .. External Functions ..
 240       LOGICAL            LSAME
 241       INTEGER            ILAENV
 242       DOUBLE PRECISION   DLAMCH, ZLANHE
 243       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANHE
 244 *     ..
 245 *     .. External Subroutines ..
 246       EXTERNAL           DSCAL, DSTERF, XERBLA, ZHETRD, ZLACPY, ZLASCL,
 247      $                   ZSTEDC, ZUNMTR
 248 *     ..
 249 *     .. Intrinsic Functions ..
 250       INTRINSIC          MAX, SQRT
 251 *     ..
 252 *     .. Executable Statements ..
 253 *
 254 *     Test the input parameters.
 255 *
 256       WANTZ = LSAME( JOBZ, 'V' )
 257       LOWER = LSAME( UPLO, 'L' )
 258       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
 259 *
 260       INFO = 0
 261       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 262          INFO = -1
 263       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
 264          INFO = -2
 265       ELSE IF( N.LT.0 ) THEN
 266          INFO = -3
 267       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 268          INFO = -5
 269       END IF
 270 *
 271       IF( INFO.EQ.0 ) THEN
 272          IF( N.LE.1 ) THEN
 273             LWMIN = 1
 274             LRWMIN = 1
 275             LIWMIN = 1
 276             LOPT = LWMIN
 277             LROPT = LRWMIN
 278             LIOPT = LIWMIN
 279          ELSE
 280             IF( WANTZ ) THEN
 281                LWMIN = 2*N + N*N
 282                LRWMIN = 1 + 5*N + 2*N**2
 283                LIWMIN = 3 + 5*N
 284             ELSE
 285                LWMIN = N + 1
 286                LRWMIN = N
 287                LIWMIN = 1
 288             END IF
 289             LOPT = MAX( LWMIN, N +
 290      $                  ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
 291             LROPT = LRWMIN
 292             LIOPT = LIWMIN
 293          END IF
 294          WORK( 1 ) = LOPT
 295          RWORK( 1 ) = LROPT
 296          IWORK( 1 ) = LIOPT
 297 *
 298          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 299             INFO = -8
 300          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
 301             INFO = -10
 302          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
 303             INFO = -12
 304          END IF
 305       END IF
 306 *
 307       IF( INFO.NE.0 ) THEN
 308          CALL XERBLA( 'ZHEEVD', -INFO )
 309          RETURN
 310       ELSE IF( LQUERY ) THEN
 311          RETURN
 312       END IF
 313 *
 314 *     Quick return if possible
 315 *
 316       IF( N.EQ.0 )
 317      $   RETURN
 318 *
 319       IF( N.EQ.1 ) THEN
 320          W( 1 ) = A( 1, 1 )
 321          IF( WANTZ )
 322      $      A( 1, 1 ) = CONE
 323          RETURN
 324       END IF
 325 *
 326 *     Get machine constants.
 327 *
 328       SAFMIN = DLAMCH( 'Safe minimum' )
 329       EPS = DLAMCH( 'Precision' )
 330       SMLNUM = SAFMIN / EPS
 331       BIGNUM = ONE / SMLNUM
 332       RMIN = SQRT( SMLNUM )
 333       RMAX = SQRT( BIGNUM )
 334 *
 335 *     Scale matrix to allowable range, if necessary.
 336 *
 337       ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK )
 338       ISCALE = 0
 339       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
 340          ISCALE = 1
 341          SIGMA = RMIN / ANRM
 342       ELSE IF( ANRM.GT.RMAX ) THEN
 343          ISCALE = 1
 344          SIGMA = RMAX / ANRM
 345       END IF
 346       IF( ISCALE.EQ.1 )
 347      $   CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
 348 *
 349 *     Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
 350 *
 351       INDE = 1
 352       INDTAU = 1
 353       INDWRK = INDTAU + N
 354       INDRWK = INDE + N
 355       INDWK2 = INDWRK + N*N
 356       LLWORK = LWORK - INDWRK + 1
 357       LLWRK2 = LWORK - INDWK2 + 1
 358       LLRWK = LRWORK - INDRWK + 1
 359       CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
 360      $             WORK( INDWRK ), LLWORK, IINFO )
 361 *
 362 *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
 363 *     ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
 364 *     tridiagonal matrix, then call ZUNMTR to multiply it to the
 365 *     Householder transformations represented as Householder vectors in
 366 *     A.
 367 *
 368       IF( .NOT.WANTZ ) THEN
 369          CALL DSTERF( N, W, RWORK( INDE ), INFO )
 370       ELSE
 371          CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N,
 372      $                WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK,
 373      $                IWORK, LIWORK, INFO )
 374          CALL ZUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
 375      $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
 376          CALL ZLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
 377       END IF
 378 *
 379 *     If matrix was scaled, then rescale eigenvalues appropriately.
 380 *
 381       IF( ISCALE.EQ.1 ) THEN
 382          IF( INFO.EQ.0 ) THEN
 383             IMAX = N
 384          ELSE
 385             IMAX = INFO - 1
 386          END IF
 387          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
 388       END IF
 389 *
 390       WORK( 1 ) = LOPT
 391       RWORK( 1 ) = LROPT
 392       IWORK( 1 ) = LIOPT
 393 *
 394       RETURN
 395 *
 396 *     End of ZHEEVD
 397 *
 398       END