SRC/chpevd.f

   1 *> \brief <b> CHPEVD 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 CHPEVD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chpevd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chpevd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chpevd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CHPEVD( 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 *       REAL               RWORK( * ), W( * )
  31 *       COMPLEX            AP( * ), WORK( * ), Z( LDZ, * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> CHPEVD 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 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 REAL 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 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 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 REAL array, dimension (MAX(1,LRWORK))
 138 *>          On exit, if INFO = 0, RWORK(1) returns the required LRWORK.
 139 *> \endverbatim
 140 *>
 141 *> \param[in] LRWORK
 142 *> \verbatim
 143 *>          LRWORK is INTEGER
 144 *>          The dimension of array RWORK.
 145 *>          If N <= 1,               LRWORK must be at least 1.
 146 *>          If JOBZ = 'N' and N > 1, LRWORK must be at least N.
 147 *>          If JOBZ = 'V' and N > 1, LRWORK must be at least
 148 *>                    1 + 5*N + 2*N**2.
 149 *>
 150 *>          If LRWORK = -1, then a workspace query is assumed; the
 151 *>          routine only calculates the required sizes of the WORK, RWORK
 152 *>          and IWORK arrays, returns these values as the first entries
 153 *>          of the WORK, RWORK and IWORK arrays, and no error message
 154 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 155 *> \endverbatim
 156 *>
 157 *> \param[out] IWORK
 158 *> \verbatim
 159 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 160 *>          On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
 161 *> \endverbatim
 162 *>
 163 *> \param[in] LIWORK
 164 *> \verbatim
 165 *>          LIWORK is INTEGER
 166 *>          The dimension of array IWORK.
 167 *>          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
 168 *>          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
 169 *>
 170 *>          If LIWORK = -1, then a workspace query is assumed; the
 171 *>          routine only calculates the required sizes of the WORK, RWORK
 172 *>          and IWORK arrays, returns these values as the first entries
 173 *>          of the WORK, RWORK and IWORK arrays, and no error message
 174 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 175 *> \endverbatim
 176 *>
 177 *> \param[out] INFO
 178 *> \verbatim
 179 *>          INFO is INTEGER
 180 *>          = 0:  successful exit
 181 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 182 *>          > 0:  if INFO = i, the algorithm failed to converge; i
 183 *>                off-diagonal elements of an intermediate tridiagonal
 184 *>                form did not converge to zero.
 185 *> \endverbatim
 186 *
 187 *  Authors:
 188 *  ========
 189 *
 190 *> \author Univ. of Tennessee
 191 *> \author Univ. of California Berkeley
 192 *> \author Univ. of Colorado Denver
 193 *> \author NAG Ltd.
 194 *
 195 *> \date November 2011
 196 *
 197 *> \ingroup complexOTHEReigen
 198 *
 199 *  =====================================================================
 200       SUBROUTINE CHPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
 201      $                   RWORK, LRWORK, IWORK, LIWORK, INFO )
 202 *
 203 *  -- LAPACK driver routine (version 3.4.0) --
 204 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 205 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 206 *     November 2011
 207 *
 208 *     .. Scalar Arguments ..
 209       CHARACTER          JOBZ, UPLO
 210       INTEGER            INFO, LDZ, LIWORK, LRWORK, LWORK, N
 211 *     ..
 212 *     .. Array Arguments ..
 213       INTEGER            IWORK( * )
 214       REAL               RWORK( * ), W( * )
 215       COMPLEX            AP( * ), WORK( * ), Z( LDZ, * )
 216 *     ..
 217 *
 218 *  =====================================================================
 219 *
 220 *     .. Parameters ..
 221       REAL               ZERO, ONE
 222       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 223       COMPLEX            CONE
 224       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
 225 *     ..
 226 *     .. Local Scalars ..
 227       LOGICAL            LQUERY, WANTZ
 228       INTEGER            IINFO, IMAX, INDE, INDRWK, INDTAU, INDWRK,
 229      $                   ISCALE, LIWMIN, LLRWK, LLWRK, LRWMIN, LWMIN
 230       REAL               ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
 231      $                   SMLNUM
 232 *     ..
 233 *     .. External Functions ..
 234       LOGICAL            LSAME
 235       REAL               CLANHP, SLAMCH
 236       EXTERNAL           LSAME, CLANHP, SLAMCH
 237 *     ..
 238 *     .. External Subroutines ..
 239       EXTERNAL           CHPTRD, CSSCAL, CSTEDC, CUPMTR, SSCAL, SSTERF,
 240      $                   XERBLA
 241 *     ..
 242 *     .. Intrinsic Functions ..
 243       INTRINSIC          SQRT
 244 *     ..
 245 *     .. Executable Statements ..
 246 *
 247 *     Test the input parameters.
 248 *
 249       WANTZ = LSAME( JOBZ, 'V' )
 250       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
 251 *
 252       INFO = 0
 253       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 254          INFO = -1
 255       ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) )
 256      $          THEN
 257          INFO = -2
 258       ELSE IF( N.LT.0 ) THEN
 259          INFO = -3
 260       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
 261          INFO = -7
 262       END IF
 263 *
 264       IF( INFO.EQ.0 ) THEN
 265          IF( N.LE.1 ) THEN
 266             LWMIN = 1
 267             LIWMIN = 1
 268             LRWMIN = 1
 269          ELSE
 270             IF( WANTZ ) THEN
 271                LWMIN = 2*N
 272                LRWMIN = 1 + 5*N + 2*N**2
 273                LIWMIN = 3 + 5*N
 274             ELSE
 275                LWMIN = N
 276                LRWMIN = N
 277                LIWMIN = 1
 278             END IF
 279          END IF
 280          WORK( 1 ) = LWMIN
 281          RWORK( 1 ) = LRWMIN
 282          IWORK( 1 ) = LIWMIN
 283 *
 284          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 285             INFO = -9
 286          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
 287             INFO = -11
 288          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
 289             INFO = -13
 290          END IF
 291       END IF
 292 *
 293       IF( INFO.NE.0 ) THEN
 294          CALL XERBLA( 'CHPEVD', -INFO )
 295          RETURN
 296       ELSE IF( LQUERY ) THEN
 297          RETURN
 298       END IF
 299 *
 300 *     Quick return if possible
 301 *
 302       IF( N.EQ.0 )
 303      $   RETURN
 304 *
 305       IF( N.EQ.1 ) THEN
 306          W( 1 ) = AP( 1 )
 307          IF( WANTZ )
 308      $      Z( 1, 1 ) = CONE
 309          RETURN
 310       END IF
 311 *
 312 *     Get machine constants.
 313 *
 314       SAFMIN = SLAMCH( 'Safe minimum' )
 315       EPS = SLAMCH( 'Precision' )
 316       SMLNUM = SAFMIN / EPS
 317       BIGNUM = ONE / SMLNUM
 318       RMIN = SQRT( SMLNUM )
 319       RMAX = SQRT( BIGNUM )
 320 *
 321 *     Scale matrix to allowable range, if necessary.
 322 *
 323       ANRM = CLANHP( 'M', UPLO, N, AP, RWORK )
 324       ISCALE = 0
 325       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
 326          ISCALE = 1
 327          SIGMA = RMIN / ANRM
 328       ELSE IF( ANRM.GT.RMAX ) THEN
 329          ISCALE = 1
 330          SIGMA = RMAX / ANRM
 331       END IF
 332       IF( ISCALE.EQ.1 ) THEN
 333          CALL CSSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
 334       END IF
 335 *
 336 *     Call CHPTRD to reduce Hermitian packed matrix to tridiagonal form.
 337 *
 338       INDE = 1
 339       INDTAU = 1
 340       INDRWK = INDE + N
 341       INDWRK = INDTAU + N
 342       LLWRK = LWORK - INDWRK + 1
 343       LLRWK = LRWORK - INDRWK + 1
 344       CALL CHPTRD( UPLO, N, AP, W, RWORK( INDE ), WORK( INDTAU ),
 345      $             IINFO )
 346 *
 347 *     For eigenvalues only, call SSTERF.  For eigenvectors, first call
 348 *     CUPGTR to generate the orthogonal matrix, then call CSTEDC.
 349 *
 350       IF( .NOT.WANTZ ) THEN
 351          CALL SSTERF( N, W, RWORK( INDE ), INFO )
 352       ELSE
 353          CALL CSTEDC( 'I', N, W, RWORK( INDE ), Z, LDZ, WORK( INDWRK ),
 354      $                LLWRK, RWORK( INDRWK ), LLRWK, IWORK, LIWORK,
 355      $                INFO )
 356          CALL CUPMTR( 'L', UPLO, 'N', N, N, AP, WORK( INDTAU ), Z, LDZ,
 357      $                WORK( INDWRK ), IINFO )
 358       END IF
 359 *
 360 *     If matrix was scaled, then rescale eigenvalues appropriately.
 361 *
 362       IF( ISCALE.EQ.1 ) THEN
 363          IF( INFO.EQ.0 ) THEN
 364             IMAX = N
 365          ELSE
 366             IMAX = INFO - 1
 367          END IF
 368          CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
 369       END IF
 370 *
 371       WORK( 1 ) = LWMIN
 372       RWORK( 1 ) = LRWMIN
 373       IWORK( 1 ) = LIWMIN
 374       RETURN
 375 *
 376 *     End of CHPEVD
 377 *
 378       END