1 SUBROUTINE SSPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, 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* March 31, 1993 7* 8* .. Scalar Arguments .. 9 CHARACTER JOBZ, UPLO 10 INTEGER INFO, LDZ, N 11* .. 12* .. Array Arguments .. 13 REAL AP( * ), W( * ), WORK( * ), Z( LDZ, * ) 14* .. 15* 16* Purpose 17* ======= 18* 19* SSPEV computes all the eigenvalues and, optionally, eigenvectors of a 20* real symmetric matrix A in packed storage. 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* AP (input/output) REAL array, dimension (N*(N+1)/2) 37* On entry, the upper or lower triangle of the symmetric matrix 38* A, packed columnwise in a linear array. The j-th column of A 39* is stored in the array AP as follows: 40* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 41* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 42* 43* On exit, AP is overwritten by values generated during the 44* reduction to tridiagonal form. If UPLO = 'U', the diagonal 45* and first superdiagonal of the tridiagonal matrix T overwrite 46* the corresponding elements of A, and if UPLO = 'L', the 47* diagonal and first subdiagonal of T overwrite the 48* corresponding elements of A. 49* 50* W (output) REAL array, dimension (N) 51* If INFO = 0, the eigenvalues in ascending order. 52* 53* Z (output) REAL array, dimension (LDZ, N) 54* If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 55* eigenvectors of the matrix A, with the i-th column of Z 56* holding the eigenvector associated with W(i). 57* If JOBZ = 'N', then Z is not referenced. 58* 59* LDZ (input) INTEGER 60* The leading dimension of the array Z. LDZ >= 1, and if 61* JOBZ = 'V', LDZ >= max(1,N). 62* 63* WORK (workspace) REAL array, dimension (3*N) 64* 65* INFO (output) INTEGER 66* = 0: successful exit. 67* < 0: if INFO = -i, the i-th argument had an illegal value. 68* > 0: if INFO = i, the algorithm failed to converge; i 69* off-diagonal elements of an intermediate tridiagonal 70* form did not converge to zero. 71* 72* ===================================================================== 73* 74* .. Parameters .. 75 REAL ZERO, ONE 76 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 77* .. 78* .. Local Scalars .. 79 LOGICAL WANTZ 80 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE 81 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 82 $ SMLNUM 83* .. 84* .. External Functions .. 85 LOGICAL LSAME 86 REAL SLAMCH, SLANSP 87 EXTERNAL LSAME, SLAMCH, SLANSP 88* .. 89* .. External Subroutines .. 90 EXTERNAL SOPGTR, SSCAL, SSPTRD, SSTEQR, SSTERF, XERBLA 91* .. 92* .. Intrinsic Functions .. 93 INTRINSIC SQRT 94* .. 95* .. Executable Statements .. 96* 97* Test the input parameters. 98* 99 WANTZ = LSAME( JOBZ, 'V' ) 100* 101 INFO = 0 102 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 103 INFO = -1 104 ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) 105 $ THEN 106 INFO = -2 107 ELSE IF( N.LT.0 ) THEN 108 INFO = -3 109 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 110 INFO = -7 111 END IF 112* 113 IF( INFO.NE.0 ) THEN 114 CALL XERBLA( 'SSPEV ', -INFO ) 115 RETURN 116 END IF 117* 118* Quick return if possible 119* 120 IF( N.EQ.0 ) 121 $ RETURN 122* 123 IF( N.EQ.1 ) THEN 124 W( 1 ) = AP( 1 ) 125 IF( WANTZ ) 126 $ Z( 1, 1 ) = ONE 127 RETURN 128 END IF 129* 130* Get machine constants. 131* 132 SAFMIN = SLAMCH( 'Safe minimum' ) 133 EPS = SLAMCH( 'Precision' ) 134 SMLNUM = SAFMIN / EPS 135 BIGNUM = ONE / SMLNUM 136 RMIN = SQRT( SMLNUM ) 137 RMAX = SQRT( BIGNUM ) 138* 139* Scale matrix to allowable range, if necessary. 140* 141 ANRM = SLANSP( 'M', UPLO, N, AP, WORK ) 142 ISCALE = 0 143 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 144 ISCALE = 1 145 SIGMA = RMIN / ANRM 146 ELSE IF( ANRM.GT.RMAX ) THEN 147 ISCALE = 1 148 SIGMA = RMAX / ANRM 149 END IF 150 IF( ISCALE.EQ.1 ) THEN 151 CALL SSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 ) 152 END IF 153* 154* Call SSPTRD to reduce symmetric packed matrix to tridiagonal form. 155* 156 INDE = 1 157 INDTAU = INDE + N 158 CALL SSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO ) 159* 160* For eigenvalues only, call SSTERF. For eigenvectors, first call 161* SOPGTR to generate the orthogonal matrix, then call SSTEQR. 162* 163 IF( .NOT.WANTZ ) THEN 164 CALL SSTERF( N, W, WORK( INDE ), INFO ) 165 ELSE 166 INDWRK = INDTAU + N 167 CALL SOPGTR( UPLO, N, AP, WORK( INDTAU ), Z, LDZ, 168 $ WORK( INDWRK ), IINFO ) 169 CALL SSTEQR( JOBZ, N, W, WORK( INDE ), Z, LDZ, WORK( INDTAU ), 170 $ INFO ) 171 END IF 172* 173* If matrix was scaled, then rescale eigenvalues appropriately. 174* 175 IF( ISCALE.EQ.1 ) THEN 176 IF( INFO.EQ.0 ) THEN 177 IMAX = N 178 ELSE 179 IMAX = INFO - 1 180 END IF 181 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 182 END IF 183* 184 RETURN 185* 186* End of SSPEV 187* 188 END 189