![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> SSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSYEVD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssyevd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssyevd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssyevd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, 00022 * LIWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBZ, UPLO 00026 * INTEGER INFO, LDA, LIWORK, LWORK, N 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * REAL A( LDA, * ), W( * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SSYEVD computes all eigenvalues and, optionally, eigenvectors of a 00040 *> real symmetric matrix A. If eigenvectors are desired, it uses a 00041 *> divide and conquer algorithm. 00042 *> 00043 *> The divide and conquer algorithm makes very mild assumptions about 00044 *> floating point arithmetic. It will work on machines with a guard 00045 *> digit in add/subtract, or on those binary machines without guard 00046 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00047 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines 00048 *> without guard digits, but we know of none. 00049 *> 00050 *> Because of large use of BLAS of level 3, SSYEVD needs N**2 more 00051 *> workspace than SSYEVX. 00052 *> \endverbatim 00053 * 00054 * Arguments: 00055 * ========== 00056 * 00057 *> \param[in] JOBZ 00058 *> \verbatim 00059 *> JOBZ is CHARACTER*1 00060 *> = 'N': Compute eigenvalues only; 00061 *> = 'V': Compute eigenvalues and eigenvectors. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] UPLO 00065 *> \verbatim 00066 *> UPLO is CHARACTER*1 00067 *> = 'U': Upper triangle of A is stored; 00068 *> = 'L': Lower triangle of A is stored. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] N 00072 *> \verbatim 00073 *> N is INTEGER 00074 *> The order of the matrix A. N >= 0. 00075 *> \endverbatim 00076 *> 00077 *> \param[in,out] A 00078 *> \verbatim 00079 *> A is REAL array, dimension (LDA, N) 00080 *> On entry, the symmetric matrix A. If UPLO = 'U', the 00081 *> leading N-by-N upper triangular part of A contains the 00082 *> upper triangular part of the matrix A. If UPLO = 'L', 00083 *> the leading N-by-N lower triangular part of A contains 00084 *> the lower triangular part of the matrix A. 00085 *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00086 *> orthonormal eigenvectors of the matrix A. 00087 *> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 00088 *> or the upper triangle (if UPLO='U') of A, including the 00089 *> diagonal, is destroyed. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] LDA 00093 *> \verbatim 00094 *> LDA is INTEGER 00095 *> The leading dimension of the array A. LDA >= max(1,N). 00096 *> \endverbatim 00097 *> 00098 *> \param[out] W 00099 *> \verbatim 00100 *> W is REAL array, dimension (N) 00101 *> If INFO = 0, the eigenvalues in ascending order. 00102 *> \endverbatim 00103 *> 00104 *> \param[out] WORK 00105 *> \verbatim 00106 *> WORK is REAL array, 00107 *> dimension (LWORK) 00108 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] LWORK 00112 *> \verbatim 00113 *> LWORK is INTEGER 00114 *> The dimension of the array WORK. 00115 *> If N <= 1, LWORK must be at least 1. 00116 *> If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1. 00117 *> If JOBZ = 'V' and N > 1, LWORK must be at least 00118 *> 1 + 6*N + 2*N**2. 00119 *> 00120 *> If LWORK = -1, then a workspace query is assumed; the routine 00121 *> only calculates the optimal sizes of the WORK and IWORK 00122 *> arrays, returns these values as the first entries of the WORK 00123 *> and IWORK arrays, and no error message related to LWORK or 00124 *> LIWORK is issued by XERBLA. 00125 *> \endverbatim 00126 *> 00127 *> \param[out] IWORK 00128 *> \verbatim 00129 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 00130 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00131 *> \endverbatim 00132 *> 00133 *> \param[in] LIWORK 00134 *> \verbatim 00135 *> LIWORK is INTEGER 00136 *> The dimension of the array IWORK. 00137 *> If N <= 1, LIWORK must be at least 1. 00138 *> If JOBZ = 'N' and N > 1, LIWORK must be at least 1. 00139 *> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 00140 *> 00141 *> If LIWORK = -1, then a workspace query is assumed; the 00142 *> routine only calculates the optimal sizes of the WORK and 00143 *> IWORK arrays, returns these values as the first entries of 00144 *> the WORK and IWORK arrays, and no error message related to 00145 *> LWORK or LIWORK is issued by XERBLA. 00146 *> \endverbatim 00147 *> 00148 *> \param[out] INFO 00149 *> \verbatim 00150 *> INFO is INTEGER 00151 *> = 0: successful exit 00152 *> < 0: if INFO = -i, the i-th argument had an illegal value 00153 *> > 0: if INFO = i and JOBZ = 'N', then the algorithm failed 00154 *> to converge; i off-diagonal elements of an intermediate 00155 *> tridiagonal form did not converge to zero; 00156 *> if INFO = i and JOBZ = 'V', then the algorithm failed 00157 *> to compute an eigenvalue while working on the submatrix 00158 *> lying in rows and columns INFO/(N+1) through 00159 *> mod(INFO,N+1). 00160 *> \endverbatim 00161 * 00162 * Authors: 00163 * ======== 00164 * 00165 *> \author Univ. of Tennessee 00166 *> \author Univ. of California Berkeley 00167 *> \author Univ. of Colorado Denver 00168 *> \author NAG Ltd. 00169 * 00170 *> \date November 2011 00171 * 00172 *> \ingroup realSYeigen 00173 * 00174 *> \par Contributors: 00175 * ================== 00176 *> 00177 *> Jeff Rutter, Computer Science Division, University of California 00178 *> at Berkeley, USA \n 00179 *> Modified by Francoise Tisseur, University of Tennessee \n 00180 *> Modified description of INFO. Sven, 16 Feb 05. \n 00181 *> 00182 * ===================================================================== 00183 SUBROUTINE SSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, 00184 $ LIWORK, INFO ) 00185 * 00186 * -- LAPACK driver routine (version 3.4.0) -- 00187 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00188 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00189 * November 2011 00190 * 00191 * .. Scalar Arguments .. 00192 CHARACTER JOBZ, UPLO 00193 INTEGER INFO, LDA, LIWORK, LWORK, N 00194 * .. 00195 * .. Array Arguments .. 00196 INTEGER IWORK( * ) 00197 REAL A( LDA, * ), W( * ), WORK( * ) 00198 * .. 00199 * 00200 * ===================================================================== 00201 * 00202 * .. Parameters .. 00203 REAL ZERO, ONE 00204 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00205 * .. 00206 * .. Local Scalars .. 00207 * 00208 LOGICAL LOWER, LQUERY, WANTZ 00209 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE, 00210 $ LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN 00211 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00212 $ SMLNUM 00213 * .. 00214 * .. External Functions .. 00215 LOGICAL LSAME 00216 INTEGER ILAENV 00217 REAL SLAMCH, SLANSY 00218 EXTERNAL ILAENV, LSAME, SLAMCH, SLANSY 00219 * .. 00220 * .. External Subroutines .. 00221 EXTERNAL SLACPY, SLASCL, SORMTR, SSCAL, SSTEDC, SSTERF, 00222 $ SSYTRD, XERBLA 00223 * .. 00224 * .. Intrinsic Functions .. 00225 INTRINSIC MAX, SQRT 00226 * .. 00227 * .. Executable Statements .. 00228 * 00229 * Test the input parameters. 00230 * 00231 WANTZ = LSAME( JOBZ, 'V' ) 00232 LOWER = LSAME( UPLO, 'L' ) 00233 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00234 * 00235 INFO = 0 00236 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00237 INFO = -1 00238 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00239 INFO = -2 00240 ELSE IF( N.LT.0 ) THEN 00241 INFO = -3 00242 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00243 INFO = -5 00244 END IF 00245 * 00246 IF( INFO.EQ.0 ) THEN 00247 IF( N.LE.1 ) THEN 00248 LIWMIN = 1 00249 LWMIN = 1 00250 LOPT = LWMIN 00251 LIOPT = LIWMIN 00252 ELSE 00253 IF( WANTZ ) THEN 00254 LIWMIN = 3 + 5*N 00255 LWMIN = 1 + 6*N + 2*N**2 00256 ELSE 00257 LIWMIN = 1 00258 LWMIN = 2*N + 1 00259 END IF 00260 LOPT = MAX( LWMIN, 2*N + 00261 $ ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 ) ) 00262 LIOPT = LIWMIN 00263 END IF 00264 WORK( 1 ) = LOPT 00265 IWORK( 1 ) = LIOPT 00266 * 00267 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00268 INFO = -8 00269 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00270 INFO = -10 00271 END IF 00272 END IF 00273 * 00274 IF( INFO.NE.0 ) THEN 00275 CALL XERBLA( 'SSYEVD', -INFO ) 00276 RETURN 00277 ELSE IF( LQUERY ) THEN 00278 RETURN 00279 END IF 00280 * 00281 * Quick return if possible 00282 * 00283 IF( N.EQ.0 ) 00284 $ RETURN 00285 * 00286 IF( N.EQ.1 ) THEN 00287 W( 1 ) = A( 1, 1 ) 00288 IF( WANTZ ) 00289 $ A( 1, 1 ) = ONE 00290 RETURN 00291 END IF 00292 * 00293 * Get machine constants. 00294 * 00295 SAFMIN = SLAMCH( 'Safe minimum' ) 00296 EPS = SLAMCH( 'Precision' ) 00297 SMLNUM = SAFMIN / EPS 00298 BIGNUM = ONE / SMLNUM 00299 RMIN = SQRT( SMLNUM ) 00300 RMAX = SQRT( BIGNUM ) 00301 * 00302 * Scale matrix to allowable range, if necessary. 00303 * 00304 ANRM = SLANSY( 'M', UPLO, N, A, LDA, WORK ) 00305 ISCALE = 0 00306 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00307 ISCALE = 1 00308 SIGMA = RMIN / ANRM 00309 ELSE IF( ANRM.GT.RMAX ) THEN 00310 ISCALE = 1 00311 SIGMA = RMAX / ANRM 00312 END IF 00313 IF( ISCALE.EQ.1 ) 00314 $ CALL SLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 00315 * 00316 * Call SSYTRD to reduce symmetric matrix to tridiagonal form. 00317 * 00318 INDE = 1 00319 INDTAU = INDE + N 00320 INDWRK = INDTAU + N 00321 LLWORK = LWORK - INDWRK + 1 00322 INDWK2 = INDWRK + N*N 00323 LLWRK2 = LWORK - INDWK2 + 1 00324 * 00325 CALL SSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ), 00326 $ WORK( INDWRK ), LLWORK, IINFO ) 00327 * 00328 * For eigenvalues only, call SSTERF. For eigenvectors, first call 00329 * SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the 00330 * tridiagonal matrix, then call SORMTR to multiply it by the 00331 * Householder transformations stored in A. 00332 * 00333 IF( .NOT.WANTZ ) THEN 00334 CALL SSTERF( N, W, WORK( INDE ), INFO ) 00335 ELSE 00336 CALL SSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N, 00337 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO ) 00338 CALL SORMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ), 00339 $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO ) 00340 CALL SLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA ) 00341 END IF 00342 * 00343 * If matrix was scaled, then rescale eigenvalues appropriately. 00344 * 00345 IF( ISCALE.EQ.1 ) 00346 $ CALL SSCAL( N, ONE / SIGMA, W, 1 ) 00347 * 00348 WORK( 1 ) = LOPT 00349 IWORK( 1 ) = LIOPT 00350 * 00351 RETURN 00352 * 00353 * End of SSYEVD 00354 * 00355 END