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