 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CHPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, 00002 $ RWORK, LRWORK, IWORK, LIWORK, 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, LDZ, LIWORK, LRWORK, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 REAL RWORK( * ), W( * ) 00016 COMPLEX AP( * ), WORK( * ), Z( LDZ, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CHPEVD computes all the eigenvalues and, optionally, eigenvectors of 00023 * a complex Hermitian matrix A in packed storage. If eigenvectors are 00024 * desired, it uses a divide and conquer algorithm. 00025 * 00026 * The divide and conquer algorithm makes very mild assumptions about 00027 * floating point arithmetic. It will work on machines with a guard 00028 * digit in add/subtract, or on those binary machines without guard 00029 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00030 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00031 * without guard digits, but we know of none. 00032 * 00033 * Arguments 00034 * ========= 00035 * 00036 * JOBZ (input) CHARACTER*1 00037 * = 'N': Compute eigenvalues only; 00038 * = 'V': Compute eigenvalues and eigenvectors. 00039 * 00040 * UPLO (input) CHARACTER*1 00041 * = 'U': Upper triangle of A is stored; 00042 * = 'L': Lower triangle of A is stored. 00043 * 00044 * N (input) INTEGER 00045 * The order of the matrix A. N >= 0. 00046 * 00047 * AP (input/output) COMPLEX array, dimension (N*(N+1)/2) 00048 * On entry, the upper or lower triangle of the Hermitian matrix 00049 * A, packed columnwise in a linear array. The j-th column of A 00050 * is stored in the array AP as follows: 00051 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00052 * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 00053 * 00054 * On exit, AP is overwritten by values generated during the 00055 * reduction to tridiagonal form. If UPLO = 'U', the diagonal 00056 * and first superdiagonal of the tridiagonal matrix T overwrite 00057 * the corresponding elements of A, and if UPLO = 'L', the 00058 * diagonal and first subdiagonal of T overwrite the 00059 * corresponding elements of A. 00060 * 00061 * W (output) REAL array, dimension (N) 00062 * If INFO = 0, the eigenvalues in ascending order. 00063 * 00064 * Z (output) COMPLEX array, dimension (LDZ, N) 00065 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 00066 * eigenvectors of the matrix A, with the i-th column of Z 00067 * holding the eigenvector associated with W(i). 00068 * If JOBZ = 'N', then Z is not referenced. 00069 * 00070 * LDZ (input) INTEGER 00071 * The leading dimension of the array Z. LDZ >= 1, and if 00072 * JOBZ = 'V', LDZ >= max(1,N). 00073 * 00074 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00075 * On exit, if INFO = 0, WORK(1) returns the required LWORK. 00076 * 00077 * LWORK (input) INTEGER 00078 * The dimension of array WORK. 00079 * If N <= 1, LWORK must be at least 1. 00080 * If JOBZ = 'N' and N > 1, LWORK must be at least N. 00081 * If JOBZ = 'V' and N > 1, LWORK must be at least 2*N. 00082 * 00083 * If LWORK = -1, then a workspace query is assumed; the routine 00084 * only calculates the required sizes of the WORK, RWORK and 00085 * IWORK arrays, returns these values as the first entries of 00086 * the WORK, RWORK and IWORK arrays, and no error message 00087 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00088 * 00089 * RWORK (workspace/output) REAL array, dimension (MAX(1,LRWORK)) 00090 * On exit, if INFO = 0, RWORK(1) returns the required LRWORK. 00091 * 00092 * LRWORK (input) INTEGER 00093 * The dimension of array RWORK. 00094 * If N <= 1, LRWORK must be at least 1. 00095 * If JOBZ = 'N' and N > 1, LRWORK must be at least N. 00096 * If JOBZ = 'V' and N > 1, LRWORK must be at least 00097 * 1 + 5*N + 2*N**2. 00098 * 00099 * If LRWORK = -1, then a workspace query is assumed; the 00100 * routine only calculates the required sizes of the WORK, RWORK 00101 * and IWORK arrays, returns these values as the first entries 00102 * of the WORK, RWORK and IWORK arrays, and no error message 00103 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00104 * 00105 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00106 * On exit, if INFO = 0, IWORK(1) returns the required LIWORK. 00107 * 00108 * LIWORK (input) INTEGER 00109 * The dimension of array IWORK. 00110 * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. 00111 * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 00112 * 00113 * If LIWORK = -1, then a workspace query is assumed; the 00114 * routine only calculates the required sizes of the WORK, RWORK 00115 * and IWORK arrays, returns these values as the first entries 00116 * of the WORK, RWORK and IWORK arrays, and no error message 00117 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00118 * 00119 * INFO (output) INTEGER 00120 * = 0: successful exit 00121 * < 0: if INFO = -i, the i-th argument had an illegal value. 00122 * > 0: if INFO = i, the algorithm failed to converge; i 00123 * off-diagonal elements of an intermediate tridiagonal 00124 * form did not converge to zero. 00125 * 00126 * ===================================================================== 00127 * 00128 * .. Parameters .. 00129 REAL ZERO, ONE 00130 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00131 COMPLEX CONE 00132 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00133 * .. 00134 * .. Local Scalars .. 00135 LOGICAL LQUERY, WANTZ 00136 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWRK, 00137 $ ISCALE, LIWMIN, LLRWK, LLWRK, LRWMIN, LWMIN 00138 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00139 $ SMLNUM 00140 * .. 00141 * .. External Functions .. 00142 LOGICAL LSAME 00143 REAL CLANHP, SLAMCH 00144 EXTERNAL LSAME, CLANHP, SLAMCH 00145 * .. 00146 * .. External Subroutines .. 00147 EXTERNAL CHPTRD, CSSCAL, CSTEDC, CUPMTR, SSCAL, SSTERF, 00148 $ XERBLA 00149 * .. 00150 * .. Intrinsic Functions .. 00151 INTRINSIC SQRT 00152 * .. 00153 * .. Executable Statements .. 00154 * 00155 * Test the input parameters. 00156 * 00157 WANTZ = LSAME( JOBZ, 'V' ) 00158 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00159 * 00160 INFO = 0 00161 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00162 INFO = -1 00163 ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) ) 00164 $ THEN 00165 INFO = -2 00166 ELSE IF( N.LT.0 ) THEN 00167 INFO = -3 00168 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00169 INFO = -7 00170 END IF 00171 * 00172 IF( INFO.EQ.0 ) THEN 00173 IF( N.LE.1 ) THEN 00174 LWMIN = 1 00175 LIWMIN = 1 00176 LRWMIN = 1 00177 ELSE 00178 IF( WANTZ ) THEN 00179 LWMIN = 2*N 00180 LRWMIN = 1 + 5*N + 2*N**2 00181 LIWMIN = 3 + 5*N 00182 ELSE 00183 LWMIN = N 00184 LRWMIN = N 00185 LIWMIN = 1 00186 END IF 00187 END IF 00188 WORK( 1 ) = LWMIN 00189 RWORK( 1 ) = LRWMIN 00190 IWORK( 1 ) = LIWMIN 00191 * 00192 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00193 INFO = -9 00194 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00195 INFO = -11 00196 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00197 INFO = -13 00198 END IF 00199 END IF 00200 * 00201 IF( INFO.NE.0 ) THEN 00202 CALL XERBLA( 'CHPEVD', -INFO ) 00203 RETURN 00204 ELSE IF( LQUERY ) THEN 00205 RETURN 00206 END IF 00207 * 00208 * Quick return if possible 00209 * 00210 IF( N.EQ.0 ) 00211 $ RETURN 00212 * 00213 IF( N.EQ.1 ) THEN 00214 W( 1 ) = AP( 1 ) 00215 IF( WANTZ ) 00216 $ Z( 1, 1 ) = CONE 00217 RETURN 00218 END IF 00219 * 00220 * Get machine constants. 00221 * 00222 SAFMIN = SLAMCH( 'Safe minimum' ) 00223 EPS = SLAMCH( 'Precision' ) 00224 SMLNUM = SAFMIN / EPS 00225 BIGNUM = ONE / SMLNUM 00226 RMIN = SQRT( SMLNUM ) 00227 RMAX = SQRT( BIGNUM ) 00228 * 00229 * Scale matrix to allowable range, if necessary. 00230 * 00231 ANRM = CLANHP( 'M', UPLO, N, AP, RWORK ) 00232 ISCALE = 0 00233 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00234 ISCALE = 1 00235 SIGMA = RMIN / ANRM 00236 ELSE IF( ANRM.GT.RMAX ) THEN 00237 ISCALE = 1 00238 SIGMA = RMAX / ANRM 00239 END IF 00240 IF( ISCALE.EQ.1 ) THEN 00241 CALL CSSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 ) 00242 END IF 00243 * 00244 * Call CHPTRD to reduce Hermitian packed matrix to tridiagonal form. 00245 * 00246 INDE = 1 00247 INDTAU = 1 00248 INDRWK = INDE + N 00249 INDWRK = INDTAU + N 00250 LLWRK = LWORK - INDWRK + 1 00251 LLRWK = LRWORK - INDRWK + 1 00252 CALL CHPTRD( UPLO, N, AP, W, RWORK( INDE ), WORK( INDTAU ), 00253 $ IINFO ) 00254 * 00255 * For eigenvalues only, call SSTERF. For eigenvectors, first call 00256 * CUPGTR to generate the orthogonal matrix, then call CSTEDC. 00257 * 00258 IF( .NOT.WANTZ ) THEN 00259 CALL SSTERF( N, W, RWORK( INDE ), INFO ) 00260 ELSE 00261 CALL CSTEDC( 'I', N, W, RWORK( INDE ), Z, LDZ, WORK( INDWRK ), 00262 $ LLWRK, RWORK( INDRWK ), LLRWK, IWORK, LIWORK, 00263 $ INFO ) 00264 CALL CUPMTR( 'L', UPLO, 'N', N, N, AP, WORK( INDTAU ), Z, LDZ, 00265 $ WORK( INDWRK ), IINFO ) 00266 END IF 00267 * 00268 * If matrix was scaled, then rescale eigenvalues appropriately. 00269 * 00270 IF( ISCALE.EQ.1 ) THEN 00271 IF( INFO.EQ.0 ) THEN 00272 IMAX = N 00273 ELSE 00274 IMAX = INFO - 1 00275 END IF 00276 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 00277 END IF 00278 * 00279 WORK( 1 ) = LWMIN 00280 RWORK( 1 ) = LRWMIN 00281 IWORK( 1 ) = LIWMIN 00282 RETURN 00283 * 00284 * End of CHPEVD 00285 * 00286 END
1.7.3 
