 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) 00002 * 00003 * -- LAPACK driver routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER JOBZ, UPLO 00010 INTEGER INFO, LDA, LWORK, N 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DSYEV computes all eigenvalues and, optionally, eigenvectors of a 00020 * real symmetric matrix A. 00021 * 00022 * Arguments 00023 * ========= 00024 * 00025 * JOBZ (input) CHARACTER*1 00026 * = 'N': Compute eigenvalues only; 00027 * = 'V': Compute eigenvalues and eigenvectors. 00028 * 00029 * UPLO (input) CHARACTER*1 00030 * = 'U': Upper triangle of A is stored; 00031 * = 'L': Lower triangle of A is stored. 00032 * 00033 * N (input) INTEGER 00034 * The order of the matrix A. N >= 0. 00035 * 00036 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00037 * On entry, the symmetric matrix A. If UPLO = 'U', the 00038 * leading N-by-N upper triangular part of A contains the 00039 * upper triangular part of the matrix A. If UPLO = 'L', 00040 * the leading N-by-N lower triangular part of A contains 00041 * the lower triangular part of the matrix A. 00042 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00043 * orthonormal eigenvectors of the matrix A. 00044 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 00045 * or the upper triangle (if UPLO='U') of A, including the 00046 * diagonal, is destroyed. 00047 * 00048 * LDA (input) INTEGER 00049 * The leading dimension of the array A. LDA >= max(1,N). 00050 * 00051 * W (output) DOUBLE PRECISION array, dimension (N) 00052 * If INFO = 0, the eigenvalues in ascending order. 00053 * 00054 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00055 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00056 * 00057 * LWORK (input) INTEGER 00058 * The length of the array WORK. LWORK >= max(1,3*N-1). 00059 * For optimal efficiency, LWORK >= (NB+2)*N, 00060 * where NB is the blocksize for DSYTRD returned by ILAENV. 00061 * 00062 * If LWORK = -1, then a workspace query is assumed; the routine 00063 * only calculates the optimal size of the WORK array, returns 00064 * this value as the first entry of the WORK array, and no error 00065 * message related to LWORK is issued by XERBLA. 00066 * 00067 * INFO (output) INTEGER 00068 * = 0: successful exit 00069 * < 0: if INFO = -i, the i-th argument had an illegal value 00070 * > 0: if INFO = i, the algorithm failed to converge; i 00071 * off-diagonal elements of an intermediate tridiagonal 00072 * form did not converge to zero. 00073 * 00074 * ===================================================================== 00075 * 00076 * .. Parameters .. 00077 DOUBLE PRECISION ZERO, ONE 00078 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00079 * .. 00080 * .. Local Scalars .. 00081 LOGICAL LOWER, LQUERY, WANTZ 00082 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, 00083 $ LLWORK, LWKOPT, NB 00084 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00085 $ SMLNUM 00086 * .. 00087 * .. External Functions .. 00088 LOGICAL LSAME 00089 INTEGER ILAENV 00090 DOUBLE PRECISION DLAMCH, DLANSY 00091 EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY 00092 * .. 00093 * .. External Subroutines .. 00094 EXTERNAL DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD, 00095 $ XERBLA 00096 * .. 00097 * .. Intrinsic Functions .. 00098 INTRINSIC MAX, SQRT 00099 * .. 00100 * .. Executable Statements .. 00101 * 00102 * Test the input parameters. 00103 * 00104 WANTZ = LSAME( JOBZ, 'V' ) 00105 LOWER = LSAME( UPLO, 'L' ) 00106 LQUERY = ( LWORK.EQ.-1 ) 00107 * 00108 INFO = 0 00109 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00110 INFO = -1 00111 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00112 INFO = -2 00113 ELSE IF( N.LT.0 ) THEN 00114 INFO = -3 00115 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00116 INFO = -5 00117 END IF 00118 * 00119 IF( INFO.EQ.0 ) THEN 00120 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) 00121 LWKOPT = MAX( 1, ( NB+2 )*N ) 00122 WORK( 1 ) = LWKOPT 00123 * 00124 IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY ) 00125 $ INFO = -8 00126 END IF 00127 * 00128 IF( INFO.NE.0 ) THEN 00129 CALL XERBLA( 'DSYEV ', -INFO ) 00130 RETURN 00131 ELSE IF( LQUERY ) THEN 00132 RETURN 00133 END IF 00134 * 00135 * Quick return if possible 00136 * 00137 IF( N.EQ.0 ) THEN 00138 RETURN 00139 END IF 00140 * 00141 IF( N.EQ.1 ) THEN 00142 W( 1 ) = A( 1, 1 ) 00143 WORK( 1 ) = 2 00144 IF( WANTZ ) 00145 $ A( 1, 1 ) = ONE 00146 RETURN 00147 END IF 00148 * 00149 * Get machine constants. 00150 * 00151 SAFMIN = DLAMCH( 'Safe minimum' ) 00152 EPS = DLAMCH( 'Precision' ) 00153 SMLNUM = SAFMIN / EPS 00154 BIGNUM = ONE / SMLNUM 00155 RMIN = SQRT( SMLNUM ) 00156 RMAX = SQRT( BIGNUM ) 00157 * 00158 * Scale matrix to allowable range, if necessary. 00159 * 00160 ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK ) 00161 ISCALE = 0 00162 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00163 ISCALE = 1 00164 SIGMA = RMIN / ANRM 00165 ELSE IF( ANRM.GT.RMAX ) THEN 00166 ISCALE = 1 00167 SIGMA = RMAX / ANRM 00168 END IF 00169 IF( ISCALE.EQ.1 ) 00170 $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 00171 * 00172 * Call DSYTRD to reduce symmetric matrix to tridiagonal form. 00173 * 00174 INDE = 1 00175 INDTAU = INDE + N 00176 INDWRK = INDTAU + N 00177 LLWORK = LWORK - INDWRK + 1 00178 CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ), 00179 $ WORK( INDWRK ), LLWORK, IINFO ) 00180 * 00181 * For eigenvalues only, call DSTERF. For eigenvectors, first call 00182 * DORGTR to generate the orthogonal matrix, then call DSTEQR. 00183 * 00184 IF( .NOT.WANTZ ) THEN 00185 CALL DSTERF( N, W, WORK( INDE ), INFO ) 00186 ELSE 00187 CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ), 00188 $ LLWORK, IINFO ) 00189 CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ), 00190 $ INFO ) 00191 END IF 00192 * 00193 * If matrix was scaled, then rescale eigenvalues appropriately. 00194 * 00195 IF( ISCALE.EQ.1 ) THEN 00196 IF( INFO.EQ.0 ) THEN 00197 IMAX = N 00198 ELSE 00199 IMAX = INFO - 1 00200 END IF 00201 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00202 END IF 00203 * 00204 * Set WORK(1) to optimal workspace size. 00205 * 00206 WORK( 1 ) = LWKOPT 00207 * 00208 RETURN 00209 * 00210 * End of DSYEV 00211 * 00212 END
1.7.3 
