|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBZ, UPLO 00011 INTEGER INFO, LDA, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION RWORK( * ), W( * ) 00015 COMPLEX*16 A( LDA, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZHEEV computes all eigenvalues and, optionally, eigenvectors of a 00022 * complex Hermitian matrix A. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * JOBZ (input) CHARACTER*1 00028 * = 'N': Compute eigenvalues only; 00029 * = 'V': Compute eigenvalues and eigenvectors. 00030 * 00031 * UPLO (input) CHARACTER*1 00032 * = 'U': Upper triangle of A is stored; 00033 * = 'L': Lower triangle of A is stored. 00034 * 00035 * N (input) INTEGER 00036 * The order of the matrix A. N >= 0. 00037 * 00038 * A (input/output) COMPLEX*16 array, dimension (LDA, N) 00039 * On entry, the Hermitian matrix A. If UPLO = 'U', the 00040 * leading N-by-N upper triangular part of A contains the 00041 * upper triangular part of the matrix A. If UPLO = 'L', 00042 * the leading N-by-N lower triangular part of A contains 00043 * the lower triangular part of the matrix A. 00044 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00045 * orthonormal eigenvectors of the matrix A. 00046 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 00047 * or the upper triangle (if UPLO='U') of A, including the 00048 * diagonal, is destroyed. 00049 * 00050 * LDA (input) INTEGER 00051 * The leading dimension of the array A. LDA >= max(1,N). 00052 * 00053 * W (output) DOUBLE PRECISION array, dimension (N) 00054 * If INFO = 0, the eigenvalues in ascending order. 00055 * 00056 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00057 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00058 * 00059 * LWORK (input) INTEGER 00060 * The length of the array WORK. LWORK >= max(1,2*N-1). 00061 * For optimal efficiency, LWORK >= (NB+1)*N, 00062 * where NB is the blocksize for ZHETRD returned by ILAENV. 00063 * 00064 * If LWORK = -1, then a workspace query is assumed; the routine 00065 * only calculates the optimal size of the WORK array, returns 00066 * this value as the first entry of the WORK array, and no error 00067 * message related to LWORK is issued by XERBLA. 00068 * 00069 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2)) 00070 * 00071 * INFO (output) INTEGER 00072 * = 0: successful exit 00073 * < 0: if INFO = -i, the i-th argument had an illegal value 00074 * > 0: if INFO = i, the algorithm failed to converge; i 00075 * off-diagonal elements of an intermediate tridiagonal 00076 * form did not converge to zero. 00077 * 00078 * ===================================================================== 00079 * 00080 * .. Parameters .. 00081 DOUBLE PRECISION ZERO, ONE 00082 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00083 COMPLEX*16 CONE 00084 PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) 00085 * .. 00086 * .. Local Scalars .. 00087 LOGICAL LOWER, LQUERY, WANTZ 00088 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, 00089 $ LLWORK, LWKOPT, NB 00090 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00091 $ SMLNUM 00092 * .. 00093 * .. External Functions .. 00094 LOGICAL LSAME 00095 INTEGER ILAENV 00096 DOUBLE PRECISION DLAMCH, ZLANHE 00097 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE 00098 * .. 00099 * .. External Subroutines .. 00100 EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLASCL, ZSTEQR, 00101 $ ZUNGTR 00102 * .. 00103 * .. Intrinsic Functions .. 00104 INTRINSIC MAX, SQRT 00105 * .. 00106 * .. Executable Statements .. 00107 * 00108 * Test the input parameters. 00109 * 00110 WANTZ = LSAME( JOBZ, 'V' ) 00111 LOWER = LSAME( UPLO, 'L' ) 00112 LQUERY = ( LWORK.EQ.-1 ) 00113 * 00114 INFO = 0 00115 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00116 INFO = -1 00117 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00118 INFO = -2 00119 ELSE IF( N.LT.0 ) THEN 00120 INFO = -3 00121 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00122 INFO = -5 00123 END IF 00124 * 00125 IF( INFO.EQ.0 ) THEN 00126 NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) 00127 LWKOPT = MAX( 1, ( NB+1 )*N ) 00128 WORK( 1 ) = LWKOPT 00129 * 00130 IF( LWORK.LT.MAX( 1, 2*N-1 ) .AND. .NOT.LQUERY ) 00131 $ INFO = -8 00132 END IF 00133 * 00134 IF( INFO.NE.0 ) THEN 00135 CALL XERBLA( 'ZHEEV ', -INFO ) 00136 RETURN 00137 ELSE IF( LQUERY ) THEN 00138 RETURN 00139 END IF 00140 * 00141 * Quick return if possible 00142 * 00143 IF( N.EQ.0 ) THEN 00144 RETURN 00145 END IF 00146 * 00147 IF( N.EQ.1 ) THEN 00148 W( 1 ) = A( 1, 1 ) 00149 WORK( 1 ) = 1 00150 IF( WANTZ ) 00151 $ A( 1, 1 ) = CONE 00152 RETURN 00153 END IF 00154 * 00155 * Get machine constants. 00156 * 00157 SAFMIN = DLAMCH( 'Safe minimum' ) 00158 EPS = DLAMCH( 'Precision' ) 00159 SMLNUM = SAFMIN / EPS 00160 BIGNUM = ONE / SMLNUM 00161 RMIN = SQRT( SMLNUM ) 00162 RMAX = SQRT( BIGNUM ) 00163 * 00164 * Scale matrix to allowable range, if necessary. 00165 * 00166 ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK ) 00167 ISCALE = 0 00168 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00169 ISCALE = 1 00170 SIGMA = RMIN / ANRM 00171 ELSE IF( ANRM.GT.RMAX ) THEN 00172 ISCALE = 1 00173 SIGMA = RMAX / ANRM 00174 END IF 00175 IF( ISCALE.EQ.1 ) 00176 $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 00177 * 00178 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form. 00179 * 00180 INDE = 1 00181 INDTAU = 1 00182 INDWRK = INDTAU + N 00183 LLWORK = LWORK - INDWRK + 1 00184 CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ), 00185 $ WORK( INDWRK ), LLWORK, IINFO ) 00186 * 00187 * For eigenvalues only, call DSTERF. For eigenvectors, first call 00188 * ZUNGTR to generate the unitary matrix, then call ZSTEQR. 00189 * 00190 IF( .NOT.WANTZ ) THEN 00191 CALL DSTERF( N, W, RWORK( INDE ), INFO ) 00192 ELSE 00193 CALL ZUNGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ), 00194 $ LLWORK, IINFO ) 00195 INDWRK = INDE + N 00196 CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), A, LDA, 00197 $ RWORK( INDWRK ), INFO ) 00198 END IF 00199 * 00200 * If matrix was scaled, then rescale eigenvalues appropriately. 00201 * 00202 IF( ISCALE.EQ.1 ) THEN 00203 IF( INFO.EQ.0 ) THEN 00204 IMAX = N 00205 ELSE 00206 IMAX = INFO - 1 00207 END IF 00208 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00209 END IF 00210 * 00211 * Set WORK(1) to optimal complex workspace size. 00212 * 00213 WORK( 1 ) = LWKOPT 00214 * 00215 RETURN 00216 * 00217 * End of ZHEEV 00218 * 00219 END
1.7.3 
