SRC/zhpevd.f

   1 *> \brief <b> ZHPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER 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 ZHPEVD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpevd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpevd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpevd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
  22 *                          RWORK, LRWORK, IWORK, LIWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ, UPLO
  26 *       INTEGER            INFO, LDZ, LIWORK, LRWORK, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IWORK( * )
  30 *       DOUBLE PRECISION   RWORK( * ), W( * )
  31 *       COMPLEX*16         AP( * ), WORK( * ), Z( LDZ, * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZHPEVD computes all the eigenvalues and, optionally, eigenvectors of
  41 *> a complex Hermitian matrix A in packed storage.  If eigenvectors are
  42 *> desired, it uses a 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] AP
  76 *> \verbatim
  77 *>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
  78 *>          On entry, the upper or lower triangle of the Hermitian matrix
  79 *>          A, packed columnwise in a linear array.  The j-th column of A
  80 *>          is stored in the array AP as follows:
  81 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  82 *>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
  83 *>
  84 *>          On exit, AP is overwritten by values generated during the
  85 *>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
  86 *>          and first superdiagonal of the tridiagonal matrix T overwrite
  87 *>          the corresponding elements of A, and if UPLO = 'L', the
  88 *>          diagonal and first subdiagonal of T overwrite the
  89 *>          corresponding elements of A.
  90 *> \endverbatim
  91 *>
  92 *> \param[out] W
  93 *> \verbatim
  94 *>          W is DOUBLE PRECISION array, dimension (N)
  95 *>          If INFO = 0, the eigenvalues in ascending order.
  96 *> \endverbatim
  97 *>
  98 *> \param[out] Z
  99 *> \verbatim
 100 *>          Z is COMPLEX*16 array, dimension (LDZ, N)
 101 *>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
 102 *>          eigenvectors of the matrix A, with the i-th column of Z
 103 *>          holding the eigenvector associated with W(i).
 104 *>          If JOBZ = 'N', then Z is not referenced.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] LDZ
 108 *> \verbatim
 109 *>          LDZ is INTEGER
 110 *>          The leading dimension of the array Z.  LDZ >= 1, and if
 111 *>          JOBZ = 'V', LDZ >= max(1,N).
 112 *> \endverbatim
 113 *>
 114 *> \param[out] WORK
 115 *> \verbatim
 116 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 117 *>          On exit, if INFO = 0, WORK(1) returns the required LWORK.
 118 *> \endverbatim
 119 *>
 120 *> \param[in] LWORK
 121 *> \verbatim
 122 *>          LWORK is INTEGER
 123 *>          The dimension of array WORK.
 124 *>          If N <= 1,               LWORK must be at least 1.
 125 *>          If JOBZ = 'N' and N > 1, LWORK must be at least N.
 126 *>          If JOBZ = 'V' and N > 1, LWORK must be at least 2*N.
 127 *>
 128 *>          If LWORK = -1, then a workspace query is assumed; the routine
 129 *>          only calculates the required sizes of the WORK, RWORK and
 130 *>          IWORK arrays, returns these values as the first entries of
 131 *>          the WORK, RWORK and IWORK arrays, and no error message
 132 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 133 *> \endverbatim
 134 *>
 135 *> \param[out] RWORK
 136 *> \verbatim
 137 *>          RWORK is DOUBLE PRECISION array,
 138 *>                                         dimension (LRWORK)
 139 *>          On exit, if INFO = 0, RWORK(1) returns the required LRWORK.
 140 *> \endverbatim
 141 *>
 142 *> \param[in] LRWORK
 143 *> \verbatim
 144 *>          LRWORK is INTEGER
 145 *>          The dimension of array RWORK.
 146 *>          If N <= 1,               LRWORK must be at least 1.
 147 *>          If JOBZ = 'N' and N > 1, LRWORK must be at least N.
 148 *>          If JOBZ = 'V' and N > 1, LRWORK must be at least
 149 *>                    1 + 5*N + 2*N**2.
 150 *>
 151 *>          If LRWORK = -1, then a workspace query is assumed; the
 152 *>          routine only calculates the required sizes of the WORK, RWORK
 153 *>          and IWORK arrays, returns these values as the first entries
 154 *>          of the WORK, RWORK and IWORK arrays, and no error message
 155 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 156 *> \endverbatim
 157 *>
 158 *> \param[out] IWORK
 159 *> \verbatim
 160 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 161 *>          On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
 162 *> \endverbatim
 163 *>
 164 *> \param[in] LIWORK
 165 *> \verbatim
 166 *>          LIWORK is INTEGER
 167 *>          The dimension of array IWORK.
 168 *>          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
 169 *>          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
 170 *>
 171 *>          If LIWORK = -1, then a workspace query is assumed; the
 172 *>          routine only calculates the required sizes of the WORK, RWORK
 173 *>          and IWORK arrays, returns these values as the first entries
 174 *>          of the WORK, RWORK and IWORK arrays, and no error message
 175 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 176 *> \endverbatim
 177 *>
 178 *> \param[out] INFO
 179 *> \verbatim
 180 *>          INFO is INTEGER
 181 *>          = 0:  successful exit
 182 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 183 *>          > 0:  if INFO = i, the algorithm failed to converge; i
 184 *>                off-diagonal elements of an intermediate tridiagonal
 185 *>                form did not converge to zero.
 186 *> \endverbatim
 187 *
 188 *  Authors:
 189 *  ========
 190 *
 191 *> \author Univ. of Tennessee
 192 *> \author Univ. of California Berkeley
 193 *> \author Univ. of Colorado Denver
 194 *> \author NAG Ltd.
 195 *
 196 *> \date November 2011
 197 *
 198 *> \ingroup complex16OTHEReigen
 199 *
 200 *  =====================================================================
 201       SUBROUTINE ZHPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
 202      $                   RWORK, LRWORK, IWORK, LIWORK, INFO )
 203 *
 204 *  -- LAPACK driver routine (version 3.4.0) --
 205 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 206 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 207 *     November 2011
 208 *
 209 *     .. Scalar Arguments ..
 210       CHARACTER          JOBZ, UPLO
 211       INTEGER            INFO, LDZ, LIWORK, LRWORK, LWORK, N
 212 *     ..
 213 *     .. Array Arguments ..
 214       INTEGER            IWORK( * )
 215       DOUBLE PRECISION   RWORK( * ), W( * )
 216       COMPLEX*16         AP( * ), WORK( * ), Z( LDZ, * )
 217 *     ..
 218 *
 219 *  =====================================================================
 220 *
 221 *     .. Parameters ..
 222       DOUBLE PRECISION   ZERO, ONE
 223       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 224       COMPLEX*16         CONE
 225       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
 226 *     ..
 227 *     .. Local Scalars ..
 228       LOGICAL            LQUERY, WANTZ
 229       INTEGER            IINFO, IMAX, INDE, INDRWK, INDTAU, INDWRK,
 230      $                   ISCALE, LIWMIN, LLRWK, LLWRK, LRWMIN, LWMIN
 231       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
 232      $                   SMLNUM
 233 *     ..
 234 *     .. External Functions ..
 235       LOGICAL            LSAME
 236       DOUBLE PRECISION   DLAMCH, ZLANHP
 237       EXTERNAL           LSAME, DLAMCH, ZLANHP
 238 *     ..
 239 *     .. External Subroutines ..
 240       EXTERNAL           DSCAL, DSTERF, XERBLA, ZDSCAL, ZHPTRD, ZSTEDC,
 241      $                   ZUPMTR
 242 *     ..
 243 *     .. Intrinsic Functions ..
 244       INTRINSIC          SQRT
 245 *     ..
 246 *     .. Executable Statements ..
 247 *
 248 *     Test the input parameters.
 249 *
 250       WANTZ = LSAME( JOBZ, 'V' )
 251       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
 252 *
 253       INFO = 0
 254       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 255          INFO = -1
 256       ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) )
 257      $          THEN
 258          INFO = -2
 259       ELSE IF( N.LT.0 ) THEN
 260          INFO = -3
 261       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
 262          INFO = -7
 263       END IF
 264 *
 265       IF( INFO.EQ.0 ) THEN
 266          IF( N.LE.1 ) THEN
 267             LWMIN = 1
 268             LIWMIN = 1
 269             LRWMIN = 1
 270          ELSE
 271             IF( WANTZ ) THEN
 272                LWMIN = 2*N
 273                LRWMIN = 1 + 5*N + 2*N**2
 274                LIWMIN = 3 + 5*N
 275             ELSE
 276                LWMIN = N
 277                LRWMIN = N
 278                LIWMIN = 1
 279             END IF
 280          END IF
 281          WORK( 1 ) = LWMIN
 282          RWORK( 1 ) = LRWMIN
 283          IWORK( 1 ) = LIWMIN
 284 *
 285          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 286             INFO = -9
 287          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
 288             INFO = -11
 289          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
 290             INFO = -13
 291          END IF
 292       END IF
 293 *
 294       IF( INFO.NE.0 ) THEN
 295          CALL XERBLA( 'ZHPEVD', -INFO )
 296          RETURN
 297       ELSE IF( LQUERY ) THEN
 298          RETURN
 299       END IF
 300 *
 301 *     Quick return if possible
 302 *
 303       IF( N.EQ.0 )
 304      $   RETURN
 305 *
 306       IF( N.EQ.1 ) THEN
 307          W( 1 ) = AP( 1 )
 308          IF( WANTZ )
 309      $      Z( 1, 1 ) = CONE
 310          RETURN
 311       END IF
 312 *
 313 *     Get machine constants.
 314 *
 315       SAFMIN = DLAMCH( 'Safe minimum' )
 316       EPS = DLAMCH( 'Precision' )
 317       SMLNUM = SAFMIN / EPS
 318       BIGNUM = ONE / SMLNUM
 319       RMIN = SQRT( SMLNUM )
 320       RMAX = SQRT( BIGNUM )
 321 *
 322 *     Scale matrix to allowable range, if necessary.
 323 *
 324       ANRM = ZLANHP( 'M', UPLO, N, AP, RWORK )
 325       ISCALE = 0
 326       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
 327          ISCALE = 1
 328          SIGMA = RMIN / ANRM
 329       ELSE IF( ANRM.GT.RMAX ) THEN
 330          ISCALE = 1
 331          SIGMA = RMAX / ANRM
 332       END IF
 333       IF( ISCALE.EQ.1 ) THEN
 334          CALL ZDSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
 335       END IF
 336 *
 337 *     Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form.
 338 *
 339       INDE = 1
 340       INDTAU = 1
 341       INDRWK = INDE + N
 342       INDWRK = INDTAU + N
 343       LLWRK = LWORK - INDWRK + 1
 344       LLRWK = LRWORK - INDRWK + 1
 345       CALL ZHPTRD( UPLO, N, AP, W, RWORK( INDE ), WORK( INDTAU ),
 346      $             IINFO )
 347 *
 348 *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
 349 *     ZUPGTR to generate the orthogonal matrix, then call ZSTEDC.
 350 *
 351       IF( .NOT.WANTZ ) THEN
 352          CALL DSTERF( N, W, RWORK( INDE ), INFO )
 353       ELSE
 354          CALL ZSTEDC( 'I', N, W, RWORK( INDE ), Z, LDZ, WORK( INDWRK ),
 355      $                LLWRK, RWORK( INDRWK ), LLRWK, IWORK, LIWORK,
 356      $                INFO )
 357          CALL ZUPMTR( 'L', UPLO, 'N', N, N, AP, WORK( INDTAU ), Z, LDZ,
 358      $                WORK( INDWRK ), IINFO )
 359       END IF
 360 *
 361 *     If matrix was scaled, then rescale eigenvalues appropriately.
 362 *
 363       IF( ISCALE.EQ.1 ) THEN
 364          IF( INFO.EQ.0 ) THEN
 365             IMAX = N
 366          ELSE
 367             IMAX = INFO - 1
 368          END IF
 369          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
 370       END IF
 371 *
 372       WORK( 1 ) = LWMIN
 373       RWORK( 1 ) = LRWMIN
 374       IWORK( 1 ) = LIWMIN
 375       RETURN
 376 *
 377 *     End of ZHPEVD
 378 *
 379       END