1*> \brief <b> DSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER 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 DSPEVD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dspevd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dspevd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dspevd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, 22* IWORK, LIWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER JOBZ, UPLO 26* INTEGER INFO, LDZ, LIWORK, LWORK, N 27* .. 28* .. Array Arguments .. 29* INTEGER IWORK( * ) 30* DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DSPEVD computes all the eigenvalues and, optionally, eigenvectors 40*> of a real symmetric matrix A in packed storage. If eigenvectors are 41*> desired, it uses a 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*> \endverbatim 50* 51* Arguments: 52* ========== 53* 54*> \param[in] JOBZ 55*> \verbatim 56*> JOBZ is CHARACTER*1 57*> = 'N': Compute eigenvalues only; 58*> = 'V': Compute eigenvalues and eigenvectors. 59*> \endverbatim 60*> 61*> \param[in] UPLO 62*> \verbatim 63*> UPLO is CHARACTER*1 64*> = 'U': Upper triangle of A is stored; 65*> = 'L': Lower triangle of A is stored. 66*> \endverbatim 67*> 68*> \param[in] N 69*> \verbatim 70*> N is INTEGER 71*> The order of the matrix A. N >= 0. 72*> \endverbatim 73*> 74*> \param[in,out] AP 75*> \verbatim 76*> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2) 77*> On entry, the upper or lower triangle of the symmetric matrix 78*> A, packed columnwise in a linear array. The j-th column of A 79*> is stored in the array AP as follows: 80*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 81*> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 82*> 83*> On exit, AP is overwritten by values generated during the 84*> reduction to tridiagonal form. If UPLO = 'U', the diagonal 85*> and first superdiagonal of the tridiagonal matrix T overwrite 86*> the corresponding elements of A, and if UPLO = 'L', the 87*> diagonal and first subdiagonal of T overwrite the 88*> corresponding elements of A. 89*> \endverbatim 90*> 91*> \param[out] W 92*> \verbatim 93*> W is DOUBLE PRECISION array, dimension (N) 94*> If INFO = 0, the eigenvalues in ascending order. 95*> \endverbatim 96*> 97*> \param[out] Z 98*> \verbatim 99*> Z is DOUBLE PRECISION array, dimension (LDZ, N) 100*> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 101*> eigenvectors of the matrix A, with the i-th column of Z 102*> holding the eigenvector associated with W(i). 103*> If JOBZ = 'N', then Z is not referenced. 104*> \endverbatim 105*> 106*> \param[in] LDZ 107*> \verbatim 108*> LDZ is INTEGER 109*> The leading dimension of the array Z. LDZ >= 1, and if 110*> JOBZ = 'V', LDZ >= max(1,N). 111*> \endverbatim 112*> 113*> \param[out] WORK 114*> \verbatim 115*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 116*> On exit, if INFO = 0, WORK(1) returns the required LWORK. 117*> \endverbatim 118*> 119*> \param[in] LWORK 120*> \verbatim 121*> LWORK is INTEGER 122*> The dimension of the array WORK. 123*> If N <= 1, LWORK must be at least 1. 124*> If JOBZ = 'N' and N > 1, LWORK must be at least 2*N. 125*> If JOBZ = 'V' and N > 1, LWORK must be at least 126*> 1 + 6*N + N**2. 127*> 128*> If LWORK = -1, then a workspace query is assumed; the routine 129*> only calculates the required sizes of the WORK and IWORK 130*> arrays, returns these values as the first entries of the WORK 131*> and IWORK arrays, and no error message related to LWORK or 132*> LIWORK is issued by XERBLA. 133*> \endverbatim 134*> 135*> \param[out] IWORK 136*> \verbatim 137*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 138*> On exit, if INFO = 0, IWORK(1) returns the required LIWORK. 139*> \endverbatim 140*> 141*> \param[in] LIWORK 142*> \verbatim 143*> LIWORK is INTEGER 144*> The dimension of the array IWORK. 145*> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. 146*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 147*> 148*> If LIWORK = -1, then a workspace query is assumed; the 149*> routine only calculates the required sizes of the WORK and 150*> IWORK arrays, returns these values as the first entries of 151*> the WORK and IWORK arrays, and no error message related to 152*> LWORK or LIWORK is issued by XERBLA. 153*> \endverbatim 154*> 155*> \param[out] INFO 156*> \verbatim 157*> INFO is INTEGER 158*> = 0: successful exit 159*> < 0: if INFO = -i, the i-th argument had an illegal value. 160*> > 0: if INFO = i, the algorithm failed to converge; i 161*> off-diagonal elements of an intermediate tridiagonal 162*> form did not converge to zero. 163*> \endverbatim 164* 165* Authors: 166* ======== 167* 168*> \author Univ. of Tennessee 169*> \author Univ. of California Berkeley 170*> \author Univ. of Colorado Denver 171*> \author NAG Ltd. 172* 173*> \date June 2017 174* 175*> \ingroup doubleOTHEReigen 176* 177* ===================================================================== 178 SUBROUTINE DSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, 179 $ IWORK, LIWORK, INFO ) 180* 181* -- LAPACK driver routine (version 3.7.1) -- 182* -- LAPACK is a software package provided by Univ. of Tennessee, -- 183* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 184* June 2017 185* 186* .. Scalar Arguments .. 187 CHARACTER JOBZ, UPLO 188 INTEGER INFO, LDZ, LIWORK, LWORK, N 189* .. 190* .. Array Arguments .. 191 INTEGER IWORK( * ) 192 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * ) 193* .. 194* 195* ===================================================================== 196* 197* .. Parameters .. 198 DOUBLE PRECISION ZERO, ONE 199 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 200* .. 201* .. Local Scalars .. 202 LOGICAL LQUERY, WANTZ 203 INTEGER IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN, 204 $ LLWORK, LWMIN 205 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 206 $ SMLNUM 207* .. 208* .. External Functions .. 209 LOGICAL LSAME 210 DOUBLE PRECISION DLAMCH, DLANSP 211 EXTERNAL LSAME, DLAMCH, DLANSP 212* .. 213* .. External Subroutines .. 214 EXTERNAL DOPMTR, DSCAL, DSPTRD, DSTEDC, DSTERF, XERBLA 215* .. 216* .. Intrinsic Functions .. 217 INTRINSIC SQRT 218* .. 219* .. Executable Statements .. 220* 221* Test the input parameters. 222* 223 WANTZ = LSAME( JOBZ, 'V' ) 224 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 225* 226 INFO = 0 227 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 228 INFO = -1 229 ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) 230 $ THEN 231 INFO = -2 232 ELSE IF( N.LT.0 ) THEN 233 INFO = -3 234 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 235 INFO = -7 236 END IF 237* 238 IF( INFO.EQ.0 ) THEN 239 IF( N.LE.1 ) THEN 240 LIWMIN = 1 241 LWMIN = 1 242 ELSE 243 IF( WANTZ ) THEN 244 LIWMIN = 3 + 5*N 245 LWMIN = 1 + 6*N + N**2 246 ELSE 247 LIWMIN = 1 248 LWMIN = 2*N 249 END IF 250 END IF 251 IWORK( 1 ) = LIWMIN 252 WORK( 1 ) = LWMIN 253* 254 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 255 INFO = -9 256 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 257 INFO = -11 258 END IF 259 END IF 260* 261 IF( INFO.NE.0 ) THEN 262 CALL XERBLA( 'DSPEVD', -INFO ) 263 RETURN 264 ELSE IF( LQUERY ) THEN 265 RETURN 266 END IF 267* 268* Quick return if possible 269* 270 IF( N.EQ.0 ) 271 $ RETURN 272* 273 IF( N.EQ.1 ) THEN 274 W( 1 ) = AP( 1 ) 275 IF( WANTZ ) 276 $ Z( 1, 1 ) = ONE 277 RETURN 278 END IF 279* 280* Get machine constants. 281* 282 SAFMIN = DLAMCH( 'Safe minimum' ) 283 EPS = DLAMCH( 'Precision' ) 284 SMLNUM = SAFMIN / EPS 285 BIGNUM = ONE / SMLNUM 286 RMIN = SQRT( SMLNUM ) 287 RMAX = SQRT( BIGNUM ) 288* 289* Scale matrix to allowable range, if necessary. 290* 291 ANRM = DLANSP( 'M', UPLO, N, AP, WORK ) 292 ISCALE = 0 293 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 294 ISCALE = 1 295 SIGMA = RMIN / ANRM 296 ELSE IF( ANRM.GT.RMAX ) THEN 297 ISCALE = 1 298 SIGMA = RMAX / ANRM 299 END IF 300 IF( ISCALE.EQ.1 ) THEN 301 CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 ) 302 END IF 303* 304* Call DSPTRD to reduce symmetric packed matrix to tridiagonal form. 305* 306 INDE = 1 307 INDTAU = INDE + N 308 CALL DSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO ) 309* 310* For eigenvalues only, call DSTERF. For eigenvectors, first call 311* DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the 312* tridiagonal matrix, then call DOPMTR to multiply it by the 313* Householder transformations represented in AP. 314* 315 IF( .NOT.WANTZ ) THEN 316 CALL DSTERF( N, W, WORK( INDE ), INFO ) 317 ELSE 318 INDWRK = INDTAU + N 319 LLWORK = LWORK - INDWRK + 1 320 CALL DSTEDC( 'I', N, W, WORK( INDE ), Z, LDZ, WORK( INDWRK ), 321 $ LLWORK, IWORK, LIWORK, INFO ) 322 CALL DOPMTR( 'L', UPLO, 'N', N, N, AP, WORK( INDTAU ), Z, LDZ, 323 $ WORK( INDWRK ), IINFO ) 324 END IF 325* 326* If matrix was scaled, then rescale eigenvalues appropriately. 327* 328 IF( ISCALE.EQ.1 ) 329 $ CALL DSCAL( N, ONE / SIGMA, W, 1 ) 330* 331 WORK( 1 ) = LWMIN 332 IWORK( 1 ) = LIWMIN 333 RETURN 334* 335* End of DSPEVD 336* 337 END 338