00001 SUBROUTINE CSTEDC( 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 REAL D( * ), E( * ), RWORK( * ) 00016 COMPLEX WORK( * ), Z( LDZ, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CSTEDC 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 CHETRD or CHPTRD or CHBTRD 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 SLAED3 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) REAL 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) REAL 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 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 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) REAL array, dimension (MAX(1,LRWORK)) 00087 * On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 00088 * 00089 * LRWORK (input) INTEGER 00090 * The dimension of the array RWORK. 00091 * If COMPZ = 'N' or N <= 1, LRWORK must be at least 1. 00092 * If COMPZ = 'V' and N > 1, LRWORK must be at least 00093 * 1 + 3*N + 2*N*lg N + 3*N**2 , 00094 * where lg( N ) = smallest integer k such 00095 * that 2**k >= N. 00096 * If COMPZ = 'I' and N > 1, LRWORK must be at least 00097 * 1 + 4*N + 2*N**2 . 00098 * Note that for COMPZ = 'I' or 'V', then if N is less than or 00099 * equal to the minimum divide size, usually 25, then LRWORK 00100 * need only be max(1,2*(N-1)). 00101 * 00102 * If LRWORK = -1, then a workspace query is assumed; the 00103 * routine only calculates the optimal sizes of the WORK, RWORK 00104 * and IWORK arrays, returns these values as the first entries 00105 * of the WORK, RWORK and IWORK arrays, and no error message 00106 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00107 * 00108 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00109 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00110 * 00111 * LIWORK (input) INTEGER 00112 * The dimension of the array IWORK. 00113 * If COMPZ = 'N' or N <= 1, LIWORK must be at least 1. 00114 * If COMPZ = 'V' or N > 1, LIWORK must be at least 00115 * 6 + 6*N + 5*N*lg N. 00116 * If COMPZ = 'I' or N > 1, LIWORK must be at least 00117 * 3 + 5*N . 00118 * Note that for COMPZ = 'I' or 'V', then if N is less than or 00119 * equal to the minimum divide size, usually 25, then LIWORK 00120 * need only be 1. 00121 * 00122 * If LIWORK = -1, then a workspace query is assumed; the 00123 * routine only calculates the optimal sizes of the WORK, RWORK 00124 * and IWORK arrays, returns these values as the first entries 00125 * of the WORK, RWORK and IWORK arrays, and no error message 00126 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00127 * 00128 * INFO (output) INTEGER 00129 * = 0: successful exit. 00130 * < 0: if INFO = -i, the i-th argument had an illegal value. 00131 * > 0: The algorithm failed to compute an eigenvalue while 00132 * working on the submatrix lying in rows and columns 00133 * INFO/(N+1) through mod(INFO,N+1). 00134 * 00135 * Further Details 00136 * =============== 00137 * 00138 * Based on contributions by 00139 * Jeff Rutter, Computer Science Division, University of California 00140 * at Berkeley, USA 00141 * 00142 * ===================================================================== 00143 * 00144 * .. Parameters .. 00145 REAL ZERO, ONE, TWO 00146 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 00147 * .. 00148 * .. Local Scalars .. 00149 LOGICAL LQUERY 00150 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL, 00151 $ LRWMIN, LWMIN, M, SMLSIZ, START 00152 REAL EPS, ORGNRM, P, TINY 00153 * .. 00154 * .. External Functions .. 00155 LOGICAL LSAME 00156 INTEGER ILAENV 00157 REAL SLAMCH, SLANST 00158 EXTERNAL ILAENV, LSAME, SLAMCH, SLANST 00159 * .. 00160 * .. External Subroutines .. 00161 EXTERNAL XERBLA, CLACPY, CLACRM, CLAED0, CSTEQR, CSWAP, 00162 $ SLASCL, SLASET, SSTEDC, SSTEQR, SSTERF 00163 * .. 00164 * .. Intrinsic Functions .. 00165 INTRINSIC ABS, INT, LOG, MAX, MOD, REAL, SQRT 00166 * .. 00167 * .. Executable Statements .. 00168 * 00169 * Test the input parameters. 00170 * 00171 INFO = 0 00172 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00173 * 00174 IF( LSAME( COMPZ, 'N' ) ) THEN 00175 ICOMPZ = 0 00176 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00177 ICOMPZ = 1 00178 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00179 ICOMPZ = 2 00180 ELSE 00181 ICOMPZ = -1 00182 END IF 00183 IF( ICOMPZ.LT.0 ) THEN 00184 INFO = -1 00185 ELSE IF( N.LT.0 ) THEN 00186 INFO = -2 00187 ELSE IF( ( LDZ.LT.1 ) .OR. 00188 $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN 00189 INFO = -6 00190 END IF 00191 * 00192 IF( INFO.EQ.0 ) THEN 00193 * 00194 * Compute the workspace requirements 00195 * 00196 SMLSIZ = ILAENV( 9, 'CSTEDC', ' ', 0, 0, 0, 0 ) 00197 IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN 00198 LWMIN = 1 00199 LIWMIN = 1 00200 LRWMIN = 1 00201 ELSE IF( N.LE.SMLSIZ ) THEN 00202 LWMIN = 1 00203 LIWMIN = 1 00204 LRWMIN = 2*( N - 1 ) 00205 ELSE IF( ICOMPZ.EQ.1 ) THEN 00206 LGN = INT( LOG( REAL( N ) ) / LOG( TWO ) ) 00207 IF( 2**LGN.LT.N ) 00208 $ LGN = LGN + 1 00209 IF( 2**LGN.LT.N ) 00210 $ LGN = LGN + 1 00211 LWMIN = N*N 00212 LRWMIN = 1 + 3*N + 2*N*LGN + 3*N**2 00213 LIWMIN = 6 + 6*N + 5*N*LGN 00214 ELSE IF( ICOMPZ.EQ.2 ) THEN 00215 LWMIN = 1 00216 LRWMIN = 1 + 4*N + 2*N**2 00217 LIWMIN = 3 + 5*N 00218 END IF 00219 WORK( 1 ) = LWMIN 00220 RWORK( 1 ) = LRWMIN 00221 IWORK( 1 ) = LIWMIN 00222 * 00223 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00224 INFO = -8 00225 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00226 INFO = -10 00227 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00228 INFO = -12 00229 END IF 00230 END IF 00231 * 00232 IF( INFO.NE.0 ) THEN 00233 CALL XERBLA( 'CSTEDC', -INFO ) 00234 RETURN 00235 ELSE IF( LQUERY ) THEN 00236 RETURN 00237 END IF 00238 * 00239 * Quick return if possible 00240 * 00241 IF( N.EQ.0 ) 00242 $ RETURN 00243 IF( N.EQ.1 ) THEN 00244 IF( ICOMPZ.NE.0 ) 00245 $ Z( 1, 1 ) = ONE 00246 RETURN 00247 END IF 00248 * 00249 * If the following conditional clause is removed, then the routine 00250 * will use the Divide and Conquer routine to compute only the 00251 * eigenvalues, which requires (3N + 3N**2) real workspace and 00252 * (2 + 5N + 2N lg(N)) integer workspace. 00253 * Since on many architectures SSTERF is much faster than any other 00254 * algorithm for finding eigenvalues only, it is used here 00255 * as the default. If the conditional clause is removed, then 00256 * information on the size of workspace needs to be changed. 00257 * 00258 * If COMPZ = 'N', use SSTERF to compute the eigenvalues. 00259 * 00260 IF( ICOMPZ.EQ.0 ) THEN 00261 CALL SSTERF( N, D, E, INFO ) 00262 GO TO 70 00263 END IF 00264 * 00265 * If N is smaller than the minimum divide size (SMLSIZ+1), then 00266 * solve the problem with another solver. 00267 * 00268 IF( N.LE.SMLSIZ ) THEN 00269 * 00270 CALL CSTEQR( COMPZ, N, D, E, Z, LDZ, RWORK, INFO ) 00271 * 00272 ELSE 00273 * 00274 * If COMPZ = 'I', we simply call SSTEDC instead. 00275 * 00276 IF( ICOMPZ.EQ.2 ) THEN 00277 CALL SLASET( 'Full', N, N, ZERO, ONE, RWORK, N ) 00278 LL = N*N + 1 00279 CALL SSTEDC( 'I', N, D, E, RWORK, N, 00280 $ RWORK( LL ), LRWORK-LL+1, IWORK, LIWORK, INFO ) 00281 DO 20 J = 1, N 00282 DO 10 I = 1, N 00283 Z( I, J ) = RWORK( ( J-1 )*N+I ) 00284 10 CONTINUE 00285 20 CONTINUE 00286 GO TO 70 00287 END IF 00288 * 00289 * From now on, only option left to be handled is COMPZ = 'V', 00290 * i.e. ICOMPZ = 1. 00291 * 00292 * Scale. 00293 * 00294 ORGNRM = SLANST( 'M', N, D, E ) 00295 IF( ORGNRM.EQ.ZERO ) 00296 $ GO TO 70 00297 * 00298 EPS = SLAMCH( 'Epsilon' ) 00299 * 00300 START = 1 00301 * 00302 * while ( START <= N ) 00303 * 00304 30 CONTINUE 00305 IF( START.LE.N ) THEN 00306 * 00307 * Let FINISH be the position of the next subdiagonal entry 00308 * such that E( FINISH ) <= TINY or FINISH = N if no such 00309 * subdiagonal exists. The matrix identified by the elements 00310 * between START and FINISH constitutes an independent 00311 * sub-problem. 00312 * 00313 FINISH = START 00314 40 CONTINUE 00315 IF( FINISH.LT.N ) THEN 00316 TINY = EPS*SQRT( ABS( D( FINISH ) ) )* 00317 $ SQRT( ABS( D( FINISH+1 ) ) ) 00318 IF( ABS( E( FINISH ) ).GT.TINY ) THEN 00319 FINISH = FINISH + 1 00320 GO TO 40 00321 END IF 00322 END IF 00323 * 00324 * (Sub) Problem determined. Compute its size and solve it. 00325 * 00326 M = FINISH - START + 1 00327 IF( M.GT.SMLSIZ ) THEN 00328 * 00329 * Scale. 00330 * 00331 ORGNRM = SLANST( 'M', M, D( START ), E( START ) ) 00332 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M, 00333 $ INFO ) 00334 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ), 00335 $ M-1, INFO ) 00336 * 00337 CALL CLAED0( N, M, D( START ), E( START ), Z( 1, START ), 00338 $ LDZ, WORK, N, RWORK, IWORK, INFO ) 00339 IF( INFO.GT.0 ) THEN 00340 INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) + 00341 $ MOD( INFO, ( M+1 ) ) + START - 1 00342 GO TO 70 00343 END IF 00344 * 00345 * Scale back. 00346 * 00347 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M, 00348 $ INFO ) 00349 * 00350 ELSE 00351 CALL SSTEQR( 'I', M, D( START ), E( START ), RWORK, M, 00352 $ RWORK( M*M+1 ), INFO ) 00353 CALL CLACRM( N, M, Z( 1, START ), LDZ, RWORK, M, WORK, N, 00354 $ RWORK( M*M+1 ) ) 00355 CALL CLACPY( 'A', N, M, WORK, N, Z( 1, START ), LDZ ) 00356 IF( INFO.GT.0 ) THEN 00357 INFO = START*( N+1 ) + FINISH 00358 GO TO 70 00359 END IF 00360 END IF 00361 * 00362 START = FINISH + 1 00363 GO TO 30 00364 END IF 00365 * 00366 * endwhile 00367 * 00368 * If the problem split any number of times, then the eigenvalues 00369 * will not be properly ordered. Here we permute the eigenvalues 00370 * (and the associated eigenvectors) into ascending order. 00371 * 00372 IF( M.NE.N ) THEN 00373 * 00374 * Use Selection Sort to minimize swaps of eigenvectors 00375 * 00376 DO 60 II = 2, N 00377 I = II - 1 00378 K = I 00379 P = D( I ) 00380 DO 50 J = II, N 00381 IF( D( J ).LT.P ) THEN 00382 K = J 00383 P = D( J ) 00384 END IF 00385 50 CONTINUE 00386 IF( K.NE.I ) THEN 00387 D( K ) = D( I ) 00388 D( I ) = P 00389 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) 00390 END IF 00391 60 CONTINUE 00392 END IF 00393 END IF 00394 * 00395 70 CONTINUE 00396 WORK( 1 ) = LWMIN 00397 RWORK( 1 ) = LRWMIN 00398 IWORK( 1 ) = LIWMIN 00399 * 00400 RETURN 00401 * 00402 * End of CSTEDC 00403 * 00404 END
1.7.2