|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DSBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, 00002 $ LWORK, IWORK, 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, KD, LDAB, LDZ, LIWORK, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DSBEVD computes all the eigenvalues and, optionally, eigenvectors of 00022 * a real symmetric band matrix A. If eigenvectors are desired, it uses 00023 * a 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 * Arguments 00033 * ========= 00034 * 00035 * JOBZ (input) CHARACTER*1 00036 * = 'N': Compute eigenvalues only; 00037 * = 'V': Compute eigenvalues and eigenvectors. 00038 * 00039 * UPLO (input) CHARACTER*1 00040 * = 'U': Upper triangle of A is stored; 00041 * = 'L': Lower triangle of A is stored. 00042 * 00043 * N (input) INTEGER 00044 * The order of the matrix A. N >= 0. 00045 * 00046 * KD (input) INTEGER 00047 * The number of superdiagonals of the matrix A if UPLO = 'U', 00048 * or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00049 * 00050 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N) 00051 * On entry, the upper or lower triangle of the symmetric band 00052 * matrix A, stored in the first KD+1 rows of the array. The 00053 * j-th column of A is stored in the j-th column of the array AB 00054 * as follows: 00055 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00056 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00057 * 00058 * On exit, AB is overwritten by values generated during the 00059 * reduction to tridiagonal form. If UPLO = 'U', the first 00060 * superdiagonal and the diagonal of the tridiagonal matrix T 00061 * are returned in rows KD and KD+1 of AB, and if UPLO = 'L', 00062 * the diagonal and first subdiagonal of T are returned in the 00063 * first two rows of AB. 00064 * 00065 * LDAB (input) INTEGER 00066 * The leading dimension of the array AB. LDAB >= KD + 1. 00067 * 00068 * W (output) DOUBLE PRECISION array, dimension (N) 00069 * If INFO = 0, the eigenvalues in ascending order. 00070 * 00071 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N) 00072 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 00073 * eigenvectors of the matrix A, with the i-th column of Z 00074 * holding the eigenvector associated with W(i). 00075 * If JOBZ = 'N', then Z is not referenced. 00076 * 00077 * LDZ (input) INTEGER 00078 * The leading dimension of the array Z. LDZ >= 1, and if 00079 * JOBZ = 'V', LDZ >= max(1,N). 00080 * 00081 * WORK (workspace/output) DOUBLE PRECISION array, 00082 * dimension (LWORK) 00083 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00084 * 00085 * LWORK (input) INTEGER 00086 * The dimension of the array WORK. 00087 * IF N <= 1, LWORK must be at least 1. 00088 * If JOBZ = 'N' and N > 2, LWORK must be at least 2*N. 00089 * If JOBZ = 'V' and N > 2, LWORK must be at least 00090 * ( 1 + 5*N + 2*N**2 ). 00091 * 00092 * If LWORK = -1, then a workspace query is assumed; the routine 00093 * only calculates the optimal sizes of the WORK and IWORK 00094 * arrays, returns these values as the first entries of the WORK 00095 * and IWORK arrays, and no error message related to LWORK or 00096 * LIWORK is issued by XERBLA. 00097 * 00098 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00099 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00100 * 00101 * LIWORK (input) INTEGER 00102 * The dimension of the array LIWORK. 00103 * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. 00104 * If JOBZ = 'V' and N > 2, LIWORK must be at least 3 + 5*N. 00105 * 00106 * If LIWORK = -1, then a workspace query is assumed; the 00107 * routine only calculates the optimal sizes of the WORK and 00108 * IWORK arrays, returns these values as the first entries of 00109 * the WORK and IWORK arrays, and no error message related to 00110 * LWORK or LIWORK is issued by XERBLA. 00111 * 00112 * INFO (output) INTEGER 00113 * = 0: successful exit 00114 * < 0: if INFO = -i, the i-th argument had an illegal value 00115 * > 0: if INFO = i, the algorithm failed to converge; i 00116 * off-diagonal elements of an intermediate tridiagonal 00117 * form did not converge to zero. 00118 * 00119 * ===================================================================== 00120 * 00121 * .. Parameters .. 00122 DOUBLE PRECISION ZERO, ONE 00123 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00124 * .. 00125 * .. Local Scalars .. 00126 LOGICAL LOWER, LQUERY, WANTZ 00127 INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN, 00128 $ LLWRK2, LWMIN 00129 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00130 $ SMLNUM 00131 * .. 00132 * .. External Functions .. 00133 LOGICAL LSAME 00134 DOUBLE PRECISION DLAMCH, DLANSB 00135 EXTERNAL LSAME, DLAMCH, DLANSB 00136 * .. 00137 * .. External Subroutines .. 00138 EXTERNAL DGEMM, DLACPY, DLASCL, DSBTRD, DSCAL, DSTEDC, 00139 $ DSTERF, XERBLA 00140 * .. 00141 * .. Intrinsic Functions .. 00142 INTRINSIC SQRT 00143 * .. 00144 * .. Executable Statements .. 00145 * 00146 * Test the input parameters. 00147 * 00148 WANTZ = LSAME( JOBZ, 'V' ) 00149 LOWER = LSAME( UPLO, 'L' ) 00150 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00151 * 00152 INFO = 0 00153 IF( N.LE.1 ) THEN 00154 LIWMIN = 1 00155 LWMIN = 1 00156 ELSE 00157 IF( WANTZ ) THEN 00158 LIWMIN = 3 + 5*N 00159 LWMIN = 1 + 5*N + 2*N**2 00160 ELSE 00161 LIWMIN = 1 00162 LWMIN = 2*N 00163 END IF 00164 END IF 00165 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00166 INFO = -1 00167 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00168 INFO = -2 00169 ELSE IF( N.LT.0 ) THEN 00170 INFO = -3 00171 ELSE IF( KD.LT.0 ) THEN 00172 INFO = -4 00173 ELSE IF( LDAB.LT.KD+1 ) THEN 00174 INFO = -6 00175 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00176 INFO = -9 00177 END IF 00178 * 00179 IF( INFO.EQ.0 ) THEN 00180 WORK( 1 ) = LWMIN 00181 IWORK( 1 ) = LIWMIN 00182 * 00183 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00184 INFO = -11 00185 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00186 INFO = -13 00187 END IF 00188 END IF 00189 * 00190 IF( INFO.NE.0 ) THEN 00191 CALL XERBLA( 'DSBEVD', -INFO ) 00192 RETURN 00193 ELSE IF( LQUERY ) THEN 00194 RETURN 00195 END IF 00196 * 00197 * Quick return if possible 00198 * 00199 IF( N.EQ.0 ) 00200 $ RETURN 00201 * 00202 IF( N.EQ.1 ) THEN 00203 W( 1 ) = AB( 1, 1 ) 00204 IF( WANTZ ) 00205 $ Z( 1, 1 ) = ONE 00206 RETURN 00207 END IF 00208 * 00209 * Get machine constants. 00210 * 00211 SAFMIN = DLAMCH( 'Safe minimum' ) 00212 EPS = DLAMCH( 'Precision' ) 00213 SMLNUM = SAFMIN / EPS 00214 BIGNUM = ONE / SMLNUM 00215 RMIN = SQRT( SMLNUM ) 00216 RMAX = SQRT( BIGNUM ) 00217 * 00218 * Scale matrix to allowable range, if necessary. 00219 * 00220 ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK ) 00221 ISCALE = 0 00222 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00223 ISCALE = 1 00224 SIGMA = RMIN / ANRM 00225 ELSE IF( ANRM.GT.RMAX ) THEN 00226 ISCALE = 1 00227 SIGMA = RMAX / ANRM 00228 END IF 00229 IF( ISCALE.EQ.1 ) THEN 00230 IF( LOWER ) THEN 00231 CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 00232 ELSE 00233 CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 00234 END IF 00235 END IF 00236 * 00237 * Call DSBTRD to reduce symmetric band matrix to tridiagonal form. 00238 * 00239 INDE = 1 00240 INDWRK = INDE + N 00241 INDWK2 = INDWRK + N*N 00242 LLWRK2 = LWORK - INDWK2 + 1 00243 CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, WORK( INDE ), Z, LDZ, 00244 $ WORK( INDWRK ), IINFO ) 00245 * 00246 * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC. 00247 * 00248 IF( .NOT.WANTZ ) THEN 00249 CALL DSTERF( N, W, WORK( INDE ), INFO ) 00250 ELSE 00251 CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N, 00252 $ WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO ) 00253 CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N, 00254 $ ZERO, WORK( INDWK2 ), N ) 00255 CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ ) 00256 END IF 00257 * 00258 * If matrix was scaled, then rescale eigenvalues appropriately. 00259 * 00260 IF( ISCALE.EQ.1 ) 00261 $ CALL DSCAL( N, ONE / SIGMA, W, 1 ) 00262 * 00263 WORK( 1 ) = LWMIN 00264 IWORK( 1 ) = LIWMIN 00265 RETURN 00266 * 00267 * End of DSBEVD 00268 * 00269 END
1.7.3 
