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