 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, 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 COMPZ 00011 INTEGER INFO, LDZ, LIWORK, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DSTEDC computes all eigenvalues and, optionally, eigenvectors of a 00022 * symmetric tridiagonal matrix using the divide and conquer method. 00023 * The eigenvectors of a full or band real symmetric matrix can also be 00024 * found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this 00025 * matrix to tridiagonal form. 00026 * 00027 * This code makes very mild assumptions about floating point 00028 * arithmetic. It will work on machines with a guard digit in 00029 * add/subtract, or on those binary machines without guard digits 00030 * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. 00031 * It could conceivably fail on hexadecimal or decimal machines 00032 * without guard digits, but we know of none. See DLAED3 for details. 00033 * 00034 * Arguments 00035 * ========= 00036 * 00037 * COMPZ (input) CHARACTER*1 00038 * = 'N': Compute eigenvalues only. 00039 * = 'I': Compute eigenvectors of tridiagonal matrix also. 00040 * = 'V': Compute eigenvectors of original dense symmetric 00041 * matrix also. On entry, Z contains the orthogonal 00042 * matrix used to reduce the original matrix to 00043 * tridiagonal form. 00044 * 00045 * N (input) INTEGER 00046 * The dimension of the symmetric tridiagonal matrix. N >= 0. 00047 * 00048 * D (input/output) DOUBLE PRECISION array, dimension (N) 00049 * On entry, the diagonal elements of the tridiagonal matrix. 00050 * On exit, if INFO = 0, the eigenvalues in ascending order. 00051 * 00052 * E (input/output) DOUBLE PRECISION array, dimension (N-1) 00053 * On entry, the subdiagonal elements of the tridiagonal matrix. 00054 * On exit, E has been destroyed. 00055 * 00056 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) 00057 * On entry, if COMPZ = 'V', then Z contains the orthogonal 00058 * matrix used in the reduction to tridiagonal form. 00059 * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the 00060 * orthonormal eigenvectors of the original symmetric matrix, 00061 * and if COMPZ = 'I', Z contains the orthonormal eigenvectors 00062 * of the symmetric tridiagonal matrix. 00063 * If COMPZ = 'N', then Z is not referenced. 00064 * 00065 * LDZ (input) INTEGER 00066 * The leading dimension of the array Z. LDZ >= 1. 00067 * If eigenvectors are desired, then LDZ >= max(1,N). 00068 * 00069 * WORK (workspace/output) DOUBLE PRECISION array, 00070 * dimension (LWORK) 00071 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00072 * 00073 * LWORK (input) INTEGER 00074 * The dimension of the array WORK. 00075 * If COMPZ = 'N' or N <= 1 then LWORK must be at least 1. 00076 * If COMPZ = 'V' and N > 1 then LWORK must be at least 00077 * ( 1 + 3*N + 2*N*lg N + 3*N**2 ), 00078 * where lg( N ) = smallest integer k such 00079 * that 2**k >= N. 00080 * If COMPZ = 'I' and N > 1 then LWORK must be at least 00081 * ( 1 + 4*N + N**2 ). 00082 * Note that for COMPZ = 'I' or 'V', then if N is less than or 00083 * equal to the minimum divide size, usually 25, then LWORK need 00084 * only be max(1,2*(N-1)). 00085 * 00086 * If LWORK = -1, then a workspace query is assumed; the routine 00087 * only calculates the optimal size of the WORK array, returns 00088 * this value as the first entry of the WORK array, and no error 00089 * message related to LWORK is issued by XERBLA. 00090 * 00091 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00092 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00093 * 00094 * LIWORK (input) INTEGER 00095 * The dimension of the array IWORK. 00096 * If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1. 00097 * If COMPZ = 'V' and N > 1 then LIWORK must be at least 00098 * ( 6 + 6*N + 5*N*lg N ). 00099 * If COMPZ = 'I' and N > 1 then LIWORK must be at least 00100 * ( 3 + 5*N ). 00101 * Note that for COMPZ = 'I' or 'V', then if N is less than or 00102 * equal to the minimum divide size, usually 25, then LIWORK 00103 * need only be 1. 00104 * 00105 * If LIWORK = -1, then a workspace query is assumed; the 00106 * routine only calculates the optimal size of the IWORK array, 00107 * returns this value as the first entry of the IWORK array, and 00108 * no error message related to LIWORK is issued by XERBLA. 00109 * 00110 * INFO (output) INTEGER 00111 * = 0: successful exit. 00112 * < 0: if INFO = -i, the i-th argument had an illegal value. 00113 * > 0: The algorithm failed to compute an eigenvalue while 00114 * working on the submatrix lying in rows and columns 00115 * INFO/(N+1) through mod(INFO,N+1). 00116 * 00117 * Further Details 00118 * =============== 00119 * 00120 * Based on contributions by 00121 * Jeff Rutter, Computer Science Division, University of California 00122 * at Berkeley, USA 00123 * Modified by Francoise Tisseur, University of Tennessee. 00124 * 00125 * ===================================================================== 00126 * 00127 * .. Parameters .. 00128 DOUBLE PRECISION ZERO, ONE, TWO 00129 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00130 * .. 00131 * .. Local Scalars .. 00132 LOGICAL LQUERY 00133 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, 00134 $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW 00135 DOUBLE PRECISION EPS, ORGNRM, P, TINY 00136 * .. 00137 * .. External Functions .. 00138 LOGICAL LSAME 00139 INTEGER ILAENV 00140 DOUBLE PRECISION DLAMCH, DLANST 00141 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST 00142 * .. 00143 * .. External Subroutines .. 00144 EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT, 00145 $ DSTEQR, DSTERF, DSWAP, XERBLA 00146 * .. 00147 * .. Intrinsic Functions .. 00148 INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT 00149 * .. 00150 * .. Executable Statements .. 00151 * 00152 * Test the input parameters. 00153 * 00154 INFO = 0 00155 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00156 * 00157 IF( LSAME( COMPZ, 'N' ) ) THEN 00158 ICOMPZ = 0 00159 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00160 ICOMPZ = 1 00161 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00162 ICOMPZ = 2 00163 ELSE 00164 ICOMPZ = -1 00165 END IF 00166 IF( ICOMPZ.LT.0 ) THEN 00167 INFO = -1 00168 ELSE IF( N.LT.0 ) THEN 00169 INFO = -2 00170 ELSE IF( ( LDZ.LT.1 ) .OR. 00171 $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN 00172 INFO = -6 00173 END IF 00174 * 00175 IF( INFO.EQ.0 ) THEN 00176 * 00177 * Compute the workspace requirements 00178 * 00179 SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 ) 00180 IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN 00181 LIWMIN = 1 00182 LWMIN = 1 00183 ELSE IF( N.LE.SMLSIZ ) THEN 00184 LIWMIN = 1 00185 LWMIN = 2*( N - 1 ) 00186 ELSE 00187 LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) ) 00188 IF( 2**LGN.LT.N ) 00189 $ LGN = LGN + 1 00190 IF( 2**LGN.LT.N ) 00191 $ LGN = LGN + 1 00192 IF( ICOMPZ.EQ.1 ) THEN 00193 LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2 00194 LIWMIN = 6 + 6*N + 5*N*LGN 00195 ELSE IF( ICOMPZ.EQ.2 ) THEN 00196 LWMIN = 1 + 4*N + N**2 00197 LIWMIN = 3 + 5*N 00198 END IF 00199 END IF 00200 WORK( 1 ) = LWMIN 00201 IWORK( 1 ) = LIWMIN 00202 * 00203 IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN 00204 INFO = -8 00205 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN 00206 INFO = -10 00207 END IF 00208 END IF 00209 * 00210 IF( INFO.NE.0 ) THEN 00211 CALL XERBLA( 'DSTEDC', -INFO ) 00212 RETURN 00213 ELSE IF (LQUERY) THEN 00214 RETURN 00215 END IF 00216 * 00217 * Quick return if possible 00218 * 00219 IF( N.EQ.0 ) 00220 $ RETURN 00221 IF( N.EQ.1 ) THEN 00222 IF( ICOMPZ.NE.0 ) 00223 $ Z( 1, 1 ) = ONE 00224 RETURN 00225 END IF 00226 * 00227 * If the following conditional clause is removed, then the routine 00228 * will use the Divide and Conquer routine to compute only the 00229 * eigenvalues, which requires (3N + 3N**2) real workspace and 00230 * (2 + 5N + 2N lg(N)) integer workspace. 00231 * Since on many architectures DSTERF is much faster than any other 00232 * algorithm for finding eigenvalues only, it is used here 00233 * as the default. If the conditional clause is removed, then 00234 * information on the size of workspace needs to be changed. 00235 * 00236 * If COMPZ = 'N', use DSTERF to compute the eigenvalues. 00237 * 00238 IF( ICOMPZ.EQ.0 ) THEN 00239 CALL DSTERF( N, D, E, INFO ) 00240 GO TO 50 00241 END IF 00242 * 00243 * If N is smaller than the minimum divide size (SMLSIZ+1), then 00244 * solve the problem with another solver. 00245 * 00246 IF( N.LE.SMLSIZ ) THEN 00247 * 00248 CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO ) 00249 * 00250 ELSE 00251 * 00252 * If COMPZ = 'V', the Z matrix must be stored elsewhere for later 00253 * use. 00254 * 00255 IF( ICOMPZ.EQ.1 ) THEN 00256 STOREZ = 1 + N*N 00257 ELSE 00258 STOREZ = 1 00259 END IF 00260 * 00261 IF( ICOMPZ.EQ.2 ) THEN 00262 CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) 00263 END IF 00264 * 00265 * Scale. 00266 * 00267 ORGNRM = DLANST( 'M', N, D, E ) 00268 IF( ORGNRM.EQ.ZERO ) 00269 $ GO TO 50 00270 * 00271 EPS = DLAMCH( 'Epsilon' ) 00272 * 00273 START = 1 00274 * 00275 * while ( START <= N ) 00276 * 00277 10 CONTINUE 00278 IF( START.LE.N ) THEN 00279 * 00280 * Let FINISH be the position of the next subdiagonal entry 00281 * such that E( FINISH ) <= TINY or FINISH = N if no such 00282 * subdiagonal exists. The matrix identified by the elements 00283 * between START and FINISH constitutes an independent 00284 * sub-problem. 00285 * 00286 FINISH = START 00287 20 CONTINUE 00288 IF( FINISH.LT.N ) THEN 00289 TINY = EPS*SQRT( ABS( D( FINISH ) ) )* 00290 $ SQRT( ABS( D( FINISH+1 ) ) ) 00291 IF( ABS( E( FINISH ) ).GT.TINY ) THEN 00292 FINISH = FINISH + 1 00293 GO TO 20 00294 END IF 00295 END IF 00296 * 00297 * (Sub) Problem determined. Compute its size and solve it. 00298 * 00299 M = FINISH - START + 1 00300 IF( M.EQ.1 ) THEN 00301 START = FINISH + 1 00302 GO TO 10 00303 END IF 00304 IF( M.GT.SMLSIZ ) THEN 00305 * 00306 * Scale. 00307 * 00308 ORGNRM = DLANST( 'M', M, D( START ), E( START ) ) 00309 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M, 00310 $ INFO ) 00311 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ), 00312 $ M-1, INFO ) 00313 * 00314 IF( ICOMPZ.EQ.1 ) THEN 00315 STRTRW = 1 00316 ELSE 00317 STRTRW = START 00318 END IF 00319 CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ), 00320 $ Z( STRTRW, START ), LDZ, WORK( 1 ), N, 00321 $ WORK( STOREZ ), IWORK, INFO ) 00322 IF( INFO.NE.0 ) THEN 00323 INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) + 00324 $ MOD( INFO, ( M+1 ) ) + START - 1 00325 GO TO 50 00326 END IF 00327 * 00328 * Scale back. 00329 * 00330 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M, 00331 $ INFO ) 00332 * 00333 ELSE 00334 IF( ICOMPZ.EQ.1 ) THEN 00335 * 00336 * Since QR won't update a Z matrix which is larger than 00337 * the length of D, we must solve the sub-problem in a 00338 * workspace and then multiply back into Z. 00339 * 00340 CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M, 00341 $ WORK( M*M+1 ), INFO ) 00342 CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ, 00343 $ WORK( STOREZ ), N ) 00344 CALL DGEMM( 'N', 'N', N, M, M, ONE, 00345 $ WORK( STOREZ ), N, WORK, M, ZERO, 00346 $ Z( 1, START ), LDZ ) 00347 ELSE IF( ICOMPZ.EQ.2 ) THEN 00348 CALL DSTEQR( 'I', M, D( START ), E( START ), 00349 $ Z( START, START ), LDZ, WORK, INFO ) 00350 ELSE 00351 CALL DSTERF( M, D( START ), E( START ), INFO ) 00352 END IF 00353 IF( INFO.NE.0 ) THEN 00354 INFO = START*( N+1 ) + FINISH 00355 GO TO 50 00356 END IF 00357 END IF 00358 * 00359 START = FINISH + 1 00360 GO TO 10 00361 END IF 00362 * 00363 * endwhile 00364 * 00365 * If the problem split any number of times, then the eigenvalues 00366 * will not be properly ordered. Here we permute the eigenvalues 00367 * (and the associated eigenvectors) into ascending order. 00368 * 00369 IF( M.NE.N ) THEN 00370 IF( ICOMPZ.EQ.0 ) THEN 00371 * 00372 * Use Quick Sort 00373 * 00374 CALL DLASRT( 'I', N, D, INFO ) 00375 * 00376 ELSE 00377 * 00378 * Use Selection Sort to minimize swaps of eigenvectors 00379 * 00380 DO 40 II = 2, N 00381 I = II - 1 00382 K = I 00383 P = D( I ) 00384 DO 30 J = II, N 00385 IF( D( J ).LT.P ) THEN 00386 K = J 00387 P = D( J ) 00388 END IF 00389 30 CONTINUE 00390 IF( K.NE.I ) THEN 00391 D( K ) = D( I ) 00392 D( I ) = P 00393 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) 00394 END IF 00395 40 CONTINUE 00396 END IF 00397 END IF 00398 END IF 00399 * 00400 50 CONTINUE 00401 WORK( 1 ) = LWMIN 00402 IWORK( 1 ) = LIWMIN 00403 * 00404 RETURN 00405 * 00406 * End of DSTEDC 00407 * 00408 END
1.7.3 
