LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cheevd.f
Go to the documentation of this file.
00001 *> \brief <b> CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CHEEVD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheevd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheevd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheevd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
00022 *                          LRWORK, IWORK, LIWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBZ, UPLO
00026 *       INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IWORK( * )
00030 *       REAL               RWORK( * ), W( * )
00031 *       COMPLEX            A( LDA, * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> CHEEVD computes all eigenvalues and, optionally, eigenvectors of a
00041 *> complex Hermitian matrix A.  If eigenvectors are desired, it uses a
00042 *> divide and conquer algorithm.
00043 *>
00044 *> The divide and conquer algorithm makes very mild assumptions about
00045 *> floating point arithmetic. It will work on machines with a guard
00046 *> digit in add/subtract, or on those binary machines without guard
00047 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00048 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
00049 *> without guard digits, but we know of none.
00050 *> \endverbatim
00051 *
00052 *  Arguments:
00053 *  ==========
00054 *
00055 *> \param[in] JOBZ
00056 *> \verbatim
00057 *>          JOBZ is CHARACTER*1
00058 *>          = 'N':  Compute eigenvalues only;
00059 *>          = 'V':  Compute eigenvalues and eigenvectors.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] UPLO
00063 *> \verbatim
00064 *>          UPLO is CHARACTER*1
00065 *>          = 'U':  Upper triangle of A is stored;
00066 *>          = 'L':  Lower triangle of A is stored.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] N
00070 *> \verbatim
00071 *>          N is INTEGER
00072 *>          The order of the matrix A.  N >= 0.
00073 *> \endverbatim
00074 *>
00075 *> \param[in,out] A
00076 *> \verbatim
00077 *>          A is COMPLEX array, dimension (LDA, N)
00078 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
00079 *>          leading N-by-N upper triangular part of A contains the
00080 *>          upper triangular part of the matrix A.  If UPLO = 'L',
00081 *>          the leading N-by-N lower triangular part of A contains
00082 *>          the lower triangular part of the matrix A.
00083 *>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
00084 *>          orthonormal eigenvectors of the matrix A.
00085 *>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
00086 *>          or the upper triangle (if UPLO='U') of A, including the
00087 *>          diagonal, is destroyed.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] LDA
00091 *> \verbatim
00092 *>          LDA is INTEGER
00093 *>          The leading dimension of the array A.  LDA >= max(1,N).
00094 *> \endverbatim
00095 *>
00096 *> \param[out] W
00097 *> \verbatim
00098 *>          W is REAL array, dimension (N)
00099 *>          If INFO = 0, the eigenvalues in ascending order.
00100 *> \endverbatim
00101 *>
00102 *> \param[out] WORK
00103 *> \verbatim
00104 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
00105 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00106 *> \endverbatim
00107 *>
00108 *> \param[in] LWORK
00109 *> \verbatim
00110 *>          LWORK is INTEGER
00111 *>          The length of the array WORK.
00112 *>          If N <= 1,                LWORK must be at least 1.
00113 *>          If JOBZ  = 'N' and N > 1, LWORK must be at least N + 1.
00114 *>          If JOBZ  = 'V' and N > 1, LWORK must be at least 2*N + N**2.
00115 *>
00116 *>          If LWORK = -1, then a workspace query is assumed; the routine
00117 *>          only calculates the optimal sizes of the WORK, RWORK and
00118 *>          IWORK arrays, returns these values as the first entries of
00119 *>          the WORK, RWORK and IWORK arrays, and no error message
00120 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
00121 *> \endverbatim
00122 *>
00123 *> \param[out] RWORK
00124 *> \verbatim
00125 *>          RWORK is REAL array,
00126 *>                                         dimension (LRWORK)
00127 *>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
00128 *> \endverbatim
00129 *>
00130 *> \param[in] LRWORK
00131 *> \verbatim
00132 *>          LRWORK is INTEGER
00133 *>          The dimension of the array RWORK.
00134 *>          If N <= 1,                LRWORK must be at least 1.
00135 *>          If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
00136 *>          If JOBZ  = 'V' and N > 1, LRWORK must be at least
00137 *>                         1 + 5*N + 2*N**2.
00138 *>
00139 *>          If LRWORK = -1, then a workspace query is assumed; the
00140 *>          routine only calculates the optimal sizes of the WORK, RWORK
00141 *>          and IWORK arrays, returns these values as the first entries
00142 *>          of the WORK, RWORK and IWORK arrays, and no error message
00143 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
00144 *> \endverbatim
00145 *>
00146 *> \param[out] IWORK
00147 *> \verbatim
00148 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
00149 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00150 *> \endverbatim
00151 *>
00152 *> \param[in] LIWORK
00153 *> \verbatim
00154 *>          LIWORK is INTEGER
00155 *>          The dimension of the array IWORK.
00156 *>          If N <= 1,                LIWORK must be at least 1.
00157 *>          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
00158 *>          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
00159 *>
00160 *>          If LIWORK = -1, then a workspace query is assumed; the
00161 *>          routine only calculates the optimal sizes of the WORK, RWORK
00162 *>          and IWORK arrays, returns these values as the first entries
00163 *>          of the WORK, RWORK and IWORK arrays, and no error message
00164 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
00165 *> \endverbatim
00166 *>
00167 *> \param[out] INFO
00168 *> \verbatim
00169 *>          INFO is INTEGER
00170 *>          = 0:  successful exit
00171 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00172 *>          > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
00173 *>                to converge; i off-diagonal elements of an intermediate
00174 *>                tridiagonal form did not converge to zero;
00175 *>                if INFO = i and JOBZ = 'V', then the algorithm failed
00176 *>                to compute an eigenvalue while working on the submatrix
00177 *>                lying in rows and columns INFO/(N+1) through
00178 *>                mod(INFO,N+1).
00179 *> \endverbatim
00180 *
00181 *  Authors:
00182 *  ========
00183 *
00184 *> \author Univ. of Tennessee 
00185 *> \author Univ. of California Berkeley 
00186 *> \author Univ. of Colorado Denver 
00187 *> \author NAG Ltd. 
00188 *
00189 *> \date November 2011
00190 *
00191 *> \ingroup complexHEeigen
00192 *
00193 *> \par Further Details:
00194 *  =====================
00195 *>
00196 *>  Modified description of INFO. Sven, 16 Feb 05.
00197 *
00198 *> \par Contributors:
00199 *  ==================
00200 *>
00201 *> Jeff Rutter, Computer Science Division, University of California
00202 *> at Berkeley, USA
00203 *>
00204 *  =====================================================================
00205       SUBROUTINE CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
00206      $                   LRWORK, IWORK, LIWORK, INFO )
00207 *
00208 *  -- LAPACK driver routine (version 3.4.0) --
00209 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00210 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00211 *     November 2011
00212 *
00213 *     .. Scalar Arguments ..
00214       CHARACTER          JOBZ, UPLO
00215       INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
00216 *     ..
00217 *     .. Array Arguments ..
00218       INTEGER            IWORK( * )
00219       REAL               RWORK( * ), W( * )
00220       COMPLEX            A( LDA, * ), WORK( * )
00221 *     ..
00222 *
00223 *  =====================================================================
00224 *
00225 *     .. Parameters ..
00226       REAL               ZERO, ONE
00227       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00228       COMPLEX            CONE
00229       PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
00230 *     ..
00231 *     .. Local Scalars ..
00232       LOGICAL            LOWER, LQUERY, WANTZ
00233       INTEGER            IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
00234      $                   INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK,
00235      $                   LLWRK2, LOPT, LROPT, LRWMIN, LWMIN
00236       REAL               ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
00237      $                   SMLNUM
00238 *     ..
00239 *     .. External Functions ..
00240       LOGICAL            LSAME
00241       INTEGER            ILAENV
00242       REAL               CLANHE, SLAMCH
00243       EXTERNAL           ILAENV, LSAME, CLANHE, SLAMCH
00244 *     ..
00245 *     .. External Subroutines ..
00246       EXTERNAL           CHETRD, CLACPY, CLASCL, CSTEDC, CUNMTR, SSCAL,
00247      $                   SSTERF, XERBLA
00248 *     ..
00249 *     .. Intrinsic Functions ..
00250       INTRINSIC          MAX, SQRT
00251 *     ..
00252 *     .. Executable Statements ..
00253 *
00254 *     Test the input parameters.
00255 *
00256       WANTZ = LSAME( JOBZ, 'V' )
00257       LOWER = LSAME( UPLO, 'L' )
00258       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00259 *
00260       INFO = 0
00261       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00262          INFO = -1
00263       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
00264          INFO = -2
00265       ELSE IF( N.LT.0 ) THEN
00266          INFO = -3
00267       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00268          INFO = -5
00269       END IF
00270 *
00271       IF( INFO.EQ.0 ) THEN
00272          IF( N.LE.1 ) THEN
00273             LWMIN = 1
00274             LRWMIN = 1
00275             LIWMIN = 1
00276             LOPT = LWMIN
00277             LROPT = LRWMIN
00278             LIOPT = LIWMIN
00279          ELSE
00280             IF( WANTZ ) THEN
00281                LWMIN = 2*N + N*N
00282                LRWMIN = 1 + 5*N + 2*N**2
00283                LIWMIN = 3 + 5*N
00284             ELSE
00285                LWMIN = N + 1
00286                LRWMIN = N
00287                LIWMIN = 1
00288             END IF
00289             LOPT = MAX( LWMIN, N +
00290      $                  ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 ) )
00291             LROPT = LRWMIN
00292             LIOPT = LIWMIN
00293          END IF
00294          WORK( 1 ) = LOPT
00295          RWORK( 1 ) = LROPT
00296          IWORK( 1 ) = LIOPT
00297 *
00298          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00299             INFO = -8
00300          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
00301             INFO = -10
00302          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00303             INFO = -12
00304          END IF
00305       END IF
00306 *
00307       IF( INFO.NE.0 ) THEN
00308          CALL XERBLA( 'CHEEVD', -INFO )
00309          RETURN 
00310       ELSE IF( LQUERY ) THEN
00311          RETURN
00312       END IF
00313 *
00314 *     Quick return if possible
00315 *
00316       IF( N.EQ.0 )
00317      $   RETURN
00318 *
00319       IF( N.EQ.1 ) THEN
00320          W( 1 ) = A( 1, 1 )
00321          IF( WANTZ )
00322      $      A( 1, 1 ) = CONE
00323          RETURN
00324       END IF
00325 *
00326 *     Get machine constants.
00327 *
00328       SAFMIN = SLAMCH( 'Safe minimum' )
00329       EPS = SLAMCH( 'Precision' )
00330       SMLNUM = SAFMIN / EPS
00331       BIGNUM = ONE / SMLNUM
00332       RMIN = SQRT( SMLNUM )
00333       RMAX = SQRT( BIGNUM )
00334 *
00335 *     Scale matrix to allowable range, if necessary.
00336 *
00337       ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK )
00338       ISCALE = 0
00339       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
00340          ISCALE = 1
00341          SIGMA = RMIN / ANRM
00342       ELSE IF( ANRM.GT.RMAX ) THEN
00343          ISCALE = 1
00344          SIGMA = RMAX / ANRM
00345       END IF
00346       IF( ISCALE.EQ.1 )
00347      $   CALL CLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
00348 *
00349 *     Call CHETRD to reduce Hermitian matrix to tridiagonal form.
00350 *
00351       INDE = 1
00352       INDTAU = 1
00353       INDWRK = INDTAU + N
00354       INDRWK = INDE + N
00355       INDWK2 = INDWRK + N*N
00356       LLWORK = LWORK - INDWRK + 1
00357       LLWRK2 = LWORK - INDWK2 + 1
00358       LLRWK = LRWORK - INDRWK + 1
00359       CALL CHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
00360      $             WORK( INDWRK ), LLWORK, IINFO )
00361 *
00362 *     For eigenvalues only, call SSTERF.  For eigenvectors, first call
00363 *     CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
00364 *     tridiagonal matrix, then call CUNMTR to multiply it to the
00365 *     Householder transformations represented as Householder vectors in
00366 *     A.
00367 *
00368       IF( .NOT.WANTZ ) THEN
00369          CALL SSTERF( N, W, RWORK( INDE ), INFO )
00370       ELSE
00371          CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N,
00372      $                WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK,
00373      $                IWORK, LIWORK, INFO )
00374          CALL CUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
00375      $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
00376          CALL CLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
00377       END IF
00378 *
00379 *     If matrix was scaled, then rescale eigenvalues appropriately.
00380 *
00381       IF( ISCALE.EQ.1 ) THEN
00382          IF( INFO.EQ.0 ) THEN
00383             IMAX = N
00384          ELSE
00385             IMAX = INFO - 1
00386          END IF
00387          CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
00388       END IF
00389 *
00390       WORK( 1 ) = LOPT
00391       RWORK( 1 ) = LROPT
00392       IWORK( 1 ) = LIOPT
00393 *
00394       RETURN
00395 *
00396 *     End of CHEEVD
00397 *
00398       END
 All Files Functions