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