|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE SSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, 00002 $ LIWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBZ, UPLO 00011 INTEGER INFO, LDA, LIWORK, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 REAL A( LDA, * ), W( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * SSYEVD computes all eigenvalues and, optionally, eigenvectors of a 00022 * real symmetric matrix A. If eigenvectors are desired, it uses a 00023 * divide and conquer algorithm. 00024 * 00025 * The divide and conquer algorithm makes very mild assumptions about 00026 * floating point arithmetic. It will work on machines with a guard 00027 * digit in add/subtract, or on those binary machines without guard 00028 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00029 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00030 * without guard digits, but we know of none. 00031 * 00032 * Because of large use of BLAS of level 3, SSYEVD needs N**2 more 00033 * workspace than SSYEVX. 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * JOBZ (input) CHARACTER*1 00039 * = 'N': Compute eigenvalues only; 00040 * = 'V': Compute eigenvalues and eigenvectors. 00041 * 00042 * UPLO (input) CHARACTER*1 00043 * = 'U': Upper triangle of A is stored; 00044 * = 'L': Lower triangle of A is stored. 00045 * 00046 * N (input) INTEGER 00047 * The order of the matrix A. N >= 0. 00048 * 00049 * A (input/output) REAL array, dimension (LDA, N) 00050 * On entry, the symmetric matrix A. If UPLO = 'U', the 00051 * leading N-by-N upper triangular part of A contains the 00052 * upper triangular part of the matrix A. If UPLO = 'L', 00053 * the leading N-by-N lower triangular part of A contains 00054 * the lower triangular part of the matrix A. 00055 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00056 * orthonormal eigenvectors of the matrix A. 00057 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 00058 * or the upper triangle (if UPLO='U') of A, including the 00059 * diagonal, is destroyed. 00060 * 00061 * LDA (input) INTEGER 00062 * The leading dimension of the array A. LDA >= max(1,N). 00063 * 00064 * W (output) REAL array, dimension (N) 00065 * If INFO = 0, the eigenvalues in ascending order. 00066 * 00067 * WORK (workspace/output) REAL array, 00068 * dimension (LWORK) 00069 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00070 * 00071 * LWORK (input) INTEGER 00072 * The dimension of the array WORK. 00073 * If N <= 1, LWORK must be at least 1. 00074 * If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1. 00075 * If JOBZ = 'V' and N > 1, LWORK must be at least 00076 * 1 + 6*N + 2*N**2. 00077 * 00078 * If LWORK = -1, then a workspace query is assumed; the routine 00079 * only calculates the optimal sizes of the WORK and IWORK 00080 * arrays, returns these values as the first entries of the WORK 00081 * and IWORK arrays, and no error message related to LWORK or 00082 * LIWORK is issued by XERBLA. 00083 * 00084 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00085 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00086 * 00087 * LIWORK (input) INTEGER 00088 * The dimension of the array IWORK. 00089 * If N <= 1, LIWORK must be at least 1. 00090 * If JOBZ = 'N' and N > 1, LIWORK must be at least 1. 00091 * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 00092 * 00093 * If LIWORK = -1, then a workspace query is assumed; the 00094 * routine only calculates the optimal sizes of the WORK and 00095 * IWORK arrays, returns these values as the first entries of 00096 * the WORK and IWORK arrays, and no error message related to 00097 * LWORK or LIWORK is issued by XERBLA. 00098 * 00099 * INFO (output) INTEGER 00100 * = 0: successful exit 00101 * < 0: if INFO = -i, the i-th argument had an illegal value 00102 * > 0: if INFO = i and JOBZ = 'N', then the algorithm failed 00103 * to converge; i off-diagonal elements of an intermediate 00104 * tridiagonal form did not converge to zero; 00105 * if INFO = i and JOBZ = 'V', then the algorithm failed 00106 * to compute an eigenvalue while working on the submatrix 00107 * lying in rows and columns INFO/(N+1) through 00108 * mod(INFO,N+1). 00109 * 00110 * Further Details 00111 * =============== 00112 * 00113 * Based on contributions by 00114 * Jeff Rutter, Computer Science Division, University of California 00115 * at Berkeley, USA 00116 * Modified by Francoise Tisseur, University of Tennessee. 00117 * 00118 * Modified description of INFO. Sven, 16 Feb 05. 00119 * ===================================================================== 00120 * 00121 * .. Parameters .. 00122 REAL ZERO, ONE 00123 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00124 * .. 00125 * .. Local Scalars .. 00126 * 00127 LOGICAL LOWER, LQUERY, WANTZ 00128 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE, 00129 $ LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN 00130 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00131 $ SMLNUM 00132 * .. 00133 * .. External Functions .. 00134 LOGICAL LSAME 00135 INTEGER ILAENV 00136 REAL SLAMCH, SLANSY 00137 EXTERNAL ILAENV, LSAME, SLAMCH, SLANSY 00138 * .. 00139 * .. External Subroutines .. 00140 EXTERNAL SLACPY, SLASCL, SORMTR, SSCAL, SSTEDC, SSTERF, 00141 $ SSYTRD, XERBLA 00142 * .. 00143 * .. Intrinsic Functions .. 00144 INTRINSIC MAX, SQRT 00145 * .. 00146 * .. Executable Statements .. 00147 * 00148 * Test the input parameters. 00149 * 00150 WANTZ = LSAME( JOBZ, 'V' ) 00151 LOWER = LSAME( UPLO, 'L' ) 00152 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00153 * 00154 INFO = 0 00155 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00156 INFO = -1 00157 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00158 INFO = -2 00159 ELSE IF( N.LT.0 ) THEN 00160 INFO = -3 00161 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00162 INFO = -5 00163 END IF 00164 * 00165 IF( INFO.EQ.0 ) THEN 00166 IF( N.LE.1 ) THEN 00167 LIWMIN = 1 00168 LWMIN = 1 00169 LOPT = LWMIN 00170 LIOPT = LIWMIN 00171 ELSE 00172 IF( WANTZ ) THEN 00173 LIWMIN = 3 + 5*N 00174 LWMIN = 1 + 6*N + 2*N**2 00175 ELSE 00176 LIWMIN = 1 00177 LWMIN = 2*N + 1 00178 END IF 00179 LOPT = MAX( LWMIN, 2*N + 00180 $ ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 ) ) 00181 LIOPT = LIWMIN 00182 END IF 00183 WORK( 1 ) = LOPT 00184 IWORK( 1 ) = LIOPT 00185 * 00186 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00187 INFO = -8 00188 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00189 INFO = -10 00190 END IF 00191 END IF 00192 * 00193 IF( INFO.NE.0 ) THEN 00194 CALL XERBLA( 'SSYEVD', -INFO ) 00195 RETURN 00196 ELSE IF( LQUERY ) THEN 00197 RETURN 00198 END IF 00199 * 00200 * Quick return if possible 00201 * 00202 IF( N.EQ.0 ) 00203 $ RETURN 00204 * 00205 IF( N.EQ.1 ) THEN 00206 W( 1 ) = A( 1, 1 ) 00207 IF( WANTZ ) 00208 $ A( 1, 1 ) = ONE 00209 RETURN 00210 END IF 00211 * 00212 * Get machine constants. 00213 * 00214 SAFMIN = SLAMCH( 'Safe minimum' ) 00215 EPS = SLAMCH( 'Precision' ) 00216 SMLNUM = SAFMIN / EPS 00217 BIGNUM = ONE / SMLNUM 00218 RMIN = SQRT( SMLNUM ) 00219 RMAX = SQRT( BIGNUM ) 00220 * 00221 * Scale matrix to allowable range, if necessary. 00222 * 00223 ANRM = SLANSY( 'M', UPLO, N, A, LDA, WORK ) 00224 ISCALE = 0 00225 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00226 ISCALE = 1 00227 SIGMA = RMIN / ANRM 00228 ELSE IF( ANRM.GT.RMAX ) THEN 00229 ISCALE = 1 00230 SIGMA = RMAX / ANRM 00231 END IF 00232 IF( ISCALE.EQ.1 ) 00233 $ CALL SLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 00234 * 00235 * Call SSYTRD to reduce symmetric matrix to tridiagonal form. 00236 * 00237 INDE = 1 00238 INDTAU = INDE + N 00239 INDWRK = INDTAU + N 00240 LLWORK = LWORK - INDWRK + 1 00241 INDWK2 = INDWRK + N*N 00242 LLWRK2 = LWORK - INDWK2 + 1 00243 * 00244 CALL SSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ), 00245 $ WORK( INDWRK ), LLWORK, IINFO ) 00246 * 00247 * For eigenvalues only, call SSTERF. For eigenvectors, first call 00248 * SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the 00249 * tridiagonal matrix, then call SORMTR to multiply it by the 00250 * Householder transformations stored in A. 00251 * 00252 IF( .NOT.WANTZ ) THEN 00253 CALL SSTERF( N, W, WORK( INDE ), INFO ) 00254 ELSE 00255 CALL SSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N, 00256 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO ) 00257 CALL SORMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ), 00258 $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO ) 00259 CALL SLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA ) 00260 END IF 00261 * 00262 * If matrix was scaled, then rescale eigenvalues appropriately. 00263 * 00264 IF( ISCALE.EQ.1 ) 00265 $ CALL SSCAL( N, ONE / SIGMA, W, 1 ) 00266 * 00267 WORK( 1 ) = LOPT 00268 IWORK( 1 ) = LIOPT 00269 * 00270 RETURN 00271 * 00272 * End of SSYEVD 00273 * 00274 END
1.7.3 
