1 SUBROUTINE SSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) 2* 3* -- LAPACK driver routine (version 3.0) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* June 30, 1999 7* 8* .. Scalar Arguments .. 9 CHARACTER JOBZ, UPLO 10 INTEGER INFO, LDA, LWORK, N 11* .. 12* .. Array Arguments .. 13 REAL A( LDA, * ), W( * ), WORK( * ) 14* .. 15* 16* Purpose 17* ======= 18* 19* SSYEV computes all eigenvalues and, optionally, eigenvectors of a 20* real symmetric matrix A. 21* 22* Arguments 23* ========= 24* 25* JOBZ (input) CHARACTER*1 26* = 'N': Compute eigenvalues only; 27* = 'V': Compute eigenvalues and eigenvectors. 28* 29* UPLO (input) CHARACTER*1 30* = 'U': Upper triangle of A is stored; 31* = 'L': Lower triangle of A is stored. 32* 33* N (input) INTEGER 34* The order of the matrix A. N >= 0. 35* 36* A (input/output) REAL array, dimension (LDA, N) 37* On entry, the symmetric matrix A. If UPLO = 'U', the 38* leading N-by-N upper triangular part of A contains the 39* upper triangular part of the matrix A. If UPLO = 'L', 40* the leading N-by-N lower triangular part of A contains 41* the lower triangular part of the matrix A. 42* On exit, if JOBZ = 'V', then if INFO = 0, A contains the 43* orthonormal eigenvectors of the matrix A. 44* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 45* or the upper triangle (if UPLO='U') of A, including the 46* diagonal, is destroyed. 47* 48* LDA (input) INTEGER 49* The leading dimension of the array A. LDA >= max(1,N). 50* 51* W (output) REAL array, dimension (N) 52* If INFO = 0, the eigenvalues in ascending order. 53* 54* WORK (workspace/output) REAL array, dimension (LWORK) 55* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 56* 57* LWORK (input) INTEGER 58* The length of the array WORK. LWORK >= max(1,3*N-1). 59* For optimal efficiency, LWORK >= (NB+2)*N, 60* where NB is the blocksize for SSYTRD returned by ILAENV. 61* 62* If LWORK = -1, then a workspace query is assumed; the routine 63* only calculates the optimal size of the WORK array, returns 64* this value as the first entry of the WORK array, and no error 65* message related to LWORK is issued by XERBLA. 66* 67* INFO (output) INTEGER 68* = 0: successful exit 69* < 0: if INFO = -i, the i-th argument had an illegal value 70* > 0: if INFO = i, the algorithm failed to converge; i 71* off-diagonal elements of an intermediate tridiagonal 72* form did not converge to zero. 73* 74* ===================================================================== 75* 76* .. Parameters .. 77 REAL ZERO, ONE 78 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 79* .. 80* .. Local Scalars .. 81 LOGICAL LOWER, LQUERY, WANTZ 82 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, 83 $ LLWORK, LOPT, LWKOPT, NB 84 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 85 $ SMLNUM 86* .. 87* .. External Functions .. 88 LOGICAL LSAME 89 INTEGER ILAENV 90 REAL SLAMCH, SLANSY 91 EXTERNAL ILAENV, LSAME, SLAMCH, SLANSY 92* .. 93* .. External Subroutines .. 94 EXTERNAL SLASCL, SORGTR, SSCAL, SSTEQR, SSTERF, SSYTRD, 95 $ XERBLA 96* .. 97* .. Intrinsic Functions .. 98 INTRINSIC MAX, SQRT 99* .. 100* .. Executable Statements .. 101* 102* Test the input parameters. 103* 104 WANTZ = LSAME( JOBZ, 'V' ) 105 LOWER = LSAME( UPLO, 'L' ) 106 LQUERY = ( LWORK.EQ.-1 ) 107* 108 INFO = 0 109 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 110 INFO = -1 111 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 112 INFO = -2 113 ELSE IF( N.LT.0 ) THEN 114 INFO = -3 115 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 116 INFO = -5 117 ELSE IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY ) THEN 118 INFO = -8 119 END IF 120* 121 IF( INFO.EQ.0 ) THEN 122 NB = ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 ) 123 LWKOPT = MAX( 1, ( NB+2 )*N ) 124 WORK( 1 ) = LWKOPT 125 END IF 126* 127 IF( INFO.NE.0 ) THEN 128 CALL XERBLA( 'SSYEV ', -INFO ) 129 RETURN 130 ELSE IF( LQUERY ) THEN 131 RETURN 132 END IF 133* 134* Quick return if possible 135* 136 IF( N.EQ.0 ) THEN 137 WORK( 1 ) = 1 138 RETURN 139 END IF 140* 141 IF( N.EQ.1 ) THEN 142 W( 1 ) = A( 1, 1 ) 143 WORK( 1 ) = 3 144 IF( WANTZ ) 145 $ A( 1, 1 ) = ONE 146 RETURN 147 END IF 148* 149* Get machine constants. 150* 151 SAFMIN = SLAMCH( 'Safe minimum' ) 152 EPS = SLAMCH( 'Precision' ) 153 SMLNUM = SAFMIN / EPS 154 BIGNUM = ONE / SMLNUM 155 RMIN = SQRT( SMLNUM ) 156 RMAX = SQRT( BIGNUM ) 157* 158* Scale matrix to allowable range, if necessary. 159* 160 ANRM = SLANSY( 'M', UPLO, N, A, LDA, WORK ) 161 ISCALE = 0 162 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 163 ISCALE = 1 164 SIGMA = RMIN / ANRM 165 ELSE IF( ANRM.GT.RMAX ) THEN 166 ISCALE = 1 167 SIGMA = RMAX / ANRM 168 END IF 169 IF( ISCALE.EQ.1 ) 170 $ CALL SLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 171* 172* Call SSYTRD to reduce symmetric matrix to tridiagonal form. 173* 174 INDE = 1 175 INDTAU = INDE + N 176 INDWRK = INDTAU + N 177 LLWORK = LWORK - INDWRK + 1 178 CALL SSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ), 179 $ WORK( INDWRK ), LLWORK, IINFO ) 180 LOPT = 2*N + WORK( INDWRK ) 181* 182* For eigenvalues only, call SSTERF. For eigenvectors, first call 183* SORGTR to generate the orthogonal matrix, then call SSTEQR. 184* 185 IF( .NOT.WANTZ ) THEN 186 CALL SSTERF( N, W, WORK( INDE ), INFO ) 187 ELSE 188 CALL SORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ), 189 $ LLWORK, IINFO ) 190 CALL SSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ), 191 $ INFO ) 192 END IF 193* 194* If matrix was scaled, then rescale eigenvalues appropriately. 195* 196 IF( ISCALE.EQ.1 ) THEN 197 IF( INFO.EQ.0 ) THEN 198 IMAX = N 199 ELSE 200 IMAX = INFO - 1 201 END IF 202 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 203 END IF 204* 205* Set WORK(1) to optimal workspace size. 206* 207 WORK( 1 ) = LWKOPT 208* 209 RETURN 210* 211* End of SSYEV 212* 213 END 214