1*> \brief <b> DSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY 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 DSYEVD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyevd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyevd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyevd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, 22* LIWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER JOBZ, UPLO 26* INTEGER INFO, LDA, LIWORK, LWORK, N 27* .. 28* .. Array Arguments .. 29* INTEGER IWORK( * ) 30* DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DSYEVD computes all eigenvalues and, optionally, eigenvectors of a 40*> real symmetric matrix A. If eigenvectors are desired, it uses a 41*> divide and conquer algorithm. 42*> 43*> The divide and conquer algorithm makes very mild assumptions about 44*> floating point arithmetic. It will work on machines with a guard 45*> digit in add/subtract, or on those binary machines without guard 46*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 47*> Cray-2. It could conceivably fail on hexadecimal or decimal machines 48*> without guard digits, but we know of none. 49*> 50*> Because of large use of BLAS of level 3, DSYEVD needs N**2 more 51*> workspace than DSYEVX. 52*> \endverbatim 53* 54* Arguments: 55* ========== 56* 57*> \param[in] JOBZ 58*> \verbatim 59*> JOBZ is CHARACTER*1 60*> = 'N': Compute eigenvalues only; 61*> = 'V': Compute eigenvalues and eigenvectors. 62*> \endverbatim 63*> 64*> \param[in] UPLO 65*> \verbatim 66*> UPLO is CHARACTER*1 67*> = 'U': Upper triangle of A is stored; 68*> = 'L': Lower triangle of A is stored. 69*> \endverbatim 70*> 71*> \param[in] N 72*> \verbatim 73*> N is INTEGER 74*> The order of the matrix A. N >= 0. 75*> \endverbatim 76*> 77*> \param[in,out] A 78*> \verbatim 79*> A is DOUBLE PRECISION array, dimension (LDA, N) 80*> On entry, the symmetric matrix A. If UPLO = 'U', the 81*> leading N-by-N upper triangular part of A contains the 82*> upper triangular part of the matrix A. If UPLO = 'L', 83*> the leading N-by-N lower triangular part of A contains 84*> the lower triangular part of the matrix A. 85*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the 86*> orthonormal eigenvectors of the matrix A. 87*> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 88*> or the upper triangle (if UPLO='U') of A, including the 89*> diagonal, is destroyed. 90*> \endverbatim 91*> 92*> \param[in] LDA 93*> \verbatim 94*> LDA is INTEGER 95*> The leading dimension of the array A. LDA >= max(1,N). 96*> \endverbatim 97*> 98*> \param[out] W 99*> \verbatim 100*> W is DOUBLE PRECISION array, dimension (N) 101*> If INFO = 0, the eigenvalues in ascending order. 102*> \endverbatim 103*> 104*> \param[out] WORK 105*> \verbatim 106*> WORK is DOUBLE PRECISION array, 107*> dimension (LWORK) 108*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 109*> \endverbatim 110*> 111*> \param[in] LWORK 112*> \verbatim 113*> LWORK is INTEGER 114*> The dimension of the array WORK. 115*> If N <= 1, LWORK must be at least 1. 116*> If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1. 117*> If JOBZ = 'V' and N > 1, LWORK must be at least 118*> 1 + 6*N + 2*N**2. 119*> 120*> If LWORK = -1, then a workspace query is assumed; the routine 121*> only calculates the optimal sizes of the WORK and IWORK 122*> arrays, returns these values as the first entries of the WORK 123*> and IWORK arrays, and no error message related to LWORK or 124*> LIWORK is issued by XERBLA. 125*> \endverbatim 126*> 127*> \param[out] IWORK 128*> \verbatim 129*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 130*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 131*> \endverbatim 132*> 133*> \param[in] LIWORK 134*> \verbatim 135*> LIWORK is INTEGER 136*> The dimension of the array IWORK. 137*> If N <= 1, LIWORK must be at least 1. 138*> If JOBZ = 'N' and N > 1, LIWORK must be at least 1. 139*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 140*> 141*> If LIWORK = -1, then a workspace query is assumed; the 142*> routine only calculates the optimal sizes of the WORK and 143*> IWORK arrays, returns these values as the first entries of 144*> the WORK and IWORK arrays, and no error message related to 145*> LWORK or LIWORK is issued by XERBLA. 146*> \endverbatim 147*> 148*> \param[out] INFO 149*> \verbatim 150*> INFO is INTEGER 151*> = 0: successful exit 152*> < 0: if INFO = -i, the i-th argument had an illegal value 153*> > 0: if INFO = i and JOBZ = 'N', then the algorithm failed 154*> to converge; i off-diagonal elements of an intermediate 155*> tridiagonal form did not converge to zero; 156*> if INFO = i and JOBZ = 'V', then the algorithm failed 157*> to compute an eigenvalue while working on the submatrix 158*> lying in rows and columns INFO/(N+1) through 159*> mod(INFO,N+1). 160*> \endverbatim 161* 162* Authors: 163* ======== 164* 165*> \author Univ. of Tennessee 166*> \author Univ. of California Berkeley 167*> \author Univ. of Colorado Denver 168*> \author NAG Ltd. 169* 170*> \date September 2012 171* 172*> \ingroup doubleSYeigen 173* 174*> \par Contributors: 175* ================== 176*> 177*> Jeff Rutter, Computer Science Division, University of California 178*> at Berkeley, USA \n 179*> Modified by Francoise Tisseur, University of Tennessee \n 180*> Modified description of INFO. Sven, 16 Feb 05. \n 181 182 183*> 184* ===================================================================== 185 SUBROUTINE DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, 186 $ LIWORK, INFO ) 187* 188* -- LAPACK driver routine (version 3.4.2) -- 189* -- LAPACK is a software package provided by Univ. of Tennessee, -- 190* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 191* September 2012 192* 193* .. Scalar Arguments .. 194 CHARACTER JOBZ, UPLO 195 INTEGER INFO, LDA, LIWORK, LWORK, N 196* .. 197* .. Array Arguments .. 198 INTEGER IWORK( * ) 199 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ) 200* .. 201* 202* ===================================================================== 203* 204* .. Parameters .. 205 DOUBLE PRECISION ZERO, ONE 206 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 207* .. 208* .. Local Scalars .. 209* 210 LOGICAL LOWER, LQUERY, WANTZ 211 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE, 212 $ LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN 213 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 214 $ SMLNUM 215* .. 216* .. External Functions .. 217 LOGICAL LSAME 218 INTEGER ILAENV 219 DOUBLE PRECISION DLAMCH, DLANSY 220 EXTERNAL LSAME, DLAMCH, DLANSY, ILAENV 221* .. 222* .. External Subroutines .. 223 EXTERNAL DLACPY, DLASCL, DORMTR, DSCAL, DSTEDC, DSTERF, 224 $ DSYTRD, XERBLA 225* .. 226* .. Intrinsic Functions .. 227 INTRINSIC MAX, SQRT 228* .. 229* .. Executable Statements .. 230* 231* Test the input parameters. 232* 233 WANTZ = LSAME( JOBZ, 'V' ) 234 LOWER = LSAME( UPLO, 'L' ) 235 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 236* 237 INFO = 0 238 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 239 INFO = -1 240 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 241 INFO = -2 242 ELSE IF( N.LT.0 ) THEN 243 INFO = -3 244 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 245 INFO = -5 246 END IF 247* 248 IF( INFO.EQ.0 ) THEN 249 IF( N.LE.1 ) THEN 250 LIWMIN = 1 251 LWMIN = 1 252 LOPT = LWMIN 253 LIOPT = LIWMIN 254 ELSE 255 IF( WANTZ ) THEN 256 LIWMIN = 3 + 5*N 257 LWMIN = 1 + 6*N + 2*N**2 258 ELSE 259 LIWMIN = 1 260 LWMIN = 2*N + 1 261 END IF 262 LOPT = MAX( LWMIN, 2*N + 263 $ ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) ) 264 LIOPT = LIWMIN 265 END IF 266 WORK( 1 ) = LOPT 267 IWORK( 1 ) = LIOPT 268* 269 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 270 INFO = -8 271 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 272 INFO = -10 273 END IF 274 END IF 275* 276 IF( INFO.NE.0 ) THEN 277 CALL XERBLA( 'DSYEVD', -INFO ) 278 RETURN 279 ELSE IF( LQUERY ) THEN 280 RETURN 281 END IF 282* 283* Quick return if possible 284* 285 IF( N.EQ.0 ) 286 $ RETURN 287* 288 IF( N.EQ.1 ) THEN 289 W( 1 ) = A( 1, 1 ) 290 IF( WANTZ ) 291 $ A( 1, 1 ) = ONE 292 RETURN 293 END IF 294* 295* Get machine constants. 296* 297 SAFMIN = DLAMCH( 'Safe minimum' ) 298 EPS = DLAMCH( 'Precision' ) 299 SMLNUM = SAFMIN / EPS 300 BIGNUM = ONE / SMLNUM 301 RMIN = SQRT( SMLNUM ) 302 RMAX = SQRT( BIGNUM ) 303* 304* Scale matrix to allowable range, if necessary. 305* 306 ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK ) 307 ISCALE = 0 308 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 309 ISCALE = 1 310 SIGMA = RMIN / ANRM 311 ELSE IF( ANRM.GT.RMAX ) THEN 312 ISCALE = 1 313 SIGMA = RMAX / ANRM 314 END IF 315 IF( ISCALE.EQ.1 ) 316 $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 317* 318* Call DSYTRD to reduce symmetric matrix to tridiagonal form. 319* 320 INDE = 1 321 INDTAU = INDE + N 322 INDWRK = INDTAU + N 323 LLWORK = LWORK - INDWRK + 1 324 INDWK2 = INDWRK + N*N 325 LLWRK2 = LWORK - INDWK2 + 1 326* 327 CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ), 328 $ WORK( INDWRK ), LLWORK, IINFO ) 329* 330* For eigenvalues only, call DSTERF. For eigenvectors, first call 331* DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the 332* tridiagonal matrix, then call DORMTR to multiply it by the 333* Householder transformations stored in A. 334* 335 IF( .NOT.WANTZ ) THEN 336 CALL DSTERF( N, W, WORK( INDE ), INFO ) 337 ELSE 338 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N, 339 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO ) 340 CALL DORMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ), 341 $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO ) 342 CALL DLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA ) 343 END IF 344* 345* If matrix was scaled, then rescale eigenvalues appropriately. 346* 347 IF( ISCALE.EQ.1 ) 348 $ CALL DSCAL( N, ONE / SIGMA, W, 1 ) 349* 350 WORK( 1 ) = LOPT 351 IWORK( 1 ) = LIOPT 352* 353 RETURN 354* 355* End of DSYEVD 356* 357 END 358