![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CSTEDC 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CSTEDC + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cstedc.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cstedc.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cstedc.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, 00022 * LRWORK, IWORK, LIWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER COMPZ 00026 * INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * REAL D( * ), E( * ), RWORK( * ) 00031 * COMPLEX WORK( * ), Z( LDZ, * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> CSTEDC computes all eigenvalues and, optionally, eigenvectors of a 00041 *> symmetric tridiagonal matrix using the divide and conquer method. 00042 *> The eigenvectors of a full or band complex Hermitian matrix can also 00043 *> be found if CHETRD or CHPTRD or CHBTRD has been used to reduce this 00044 *> matrix to tridiagonal form. 00045 *> 00046 *> This code makes very mild assumptions about floating point 00047 *> arithmetic. It will work on machines with a guard digit in 00048 *> add/subtract, or on those binary machines without guard digits 00049 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. 00050 *> It could conceivably fail on hexadecimal or decimal machines 00051 *> without guard digits, but we know of none. See SLAED3 for details. 00052 *> \endverbatim 00053 * 00054 * Arguments: 00055 * ========== 00056 * 00057 *> \param[in] COMPZ 00058 *> \verbatim 00059 *> COMPZ is CHARACTER*1 00060 *> = 'N': Compute eigenvalues only. 00061 *> = 'I': Compute eigenvectors of tridiagonal matrix also. 00062 *> = 'V': Compute eigenvectors of original Hermitian matrix 00063 *> also. On entry, Z contains the unitary matrix used 00064 *> to reduce the original matrix to tridiagonal form. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] N 00068 *> \verbatim 00069 *> N is INTEGER 00070 *> The dimension of the symmetric tridiagonal matrix. N >= 0. 00071 *> \endverbatim 00072 *> 00073 *> \param[in,out] D 00074 *> \verbatim 00075 *> D is REAL array, dimension (N) 00076 *> On entry, the diagonal elements of the tridiagonal matrix. 00077 *> On exit, if INFO = 0, the eigenvalues in ascending order. 00078 *> \endverbatim 00079 *> 00080 *> \param[in,out] E 00081 *> \verbatim 00082 *> E is REAL array, dimension (N-1) 00083 *> On entry, the subdiagonal elements of the tridiagonal matrix. 00084 *> On exit, E has been destroyed. 00085 *> \endverbatim 00086 *> 00087 *> \param[in,out] Z 00088 *> \verbatim 00089 *> Z is COMPLEX array, dimension (LDZ,N) 00090 *> On entry, if COMPZ = 'V', then Z contains the unitary 00091 *> matrix used in the reduction to tridiagonal form. 00092 *> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the 00093 *> orthonormal eigenvectors of the original Hermitian matrix, 00094 *> and if COMPZ = 'I', Z contains the orthonormal eigenvectors 00095 *> of the symmetric tridiagonal matrix. 00096 *> If COMPZ = 'N', then Z is not referenced. 00097 *> \endverbatim 00098 *> 00099 *> \param[in] LDZ 00100 *> \verbatim 00101 *> LDZ is INTEGER 00102 *> The leading dimension of the array Z. LDZ >= 1. 00103 *> If eigenvectors are desired, then LDZ >= max(1,N). 00104 *> \endverbatim 00105 *> 00106 *> \param[out] WORK 00107 *> \verbatim 00108 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00109 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00110 *> \endverbatim 00111 *> 00112 *> \param[in] LWORK 00113 *> \verbatim 00114 *> LWORK is INTEGER 00115 *> The dimension of the array WORK. 00116 *> If COMPZ = 'N' or 'I', or N <= 1, LWORK must be at least 1. 00117 *> If COMPZ = 'V' and N > 1, LWORK must be at least N*N. 00118 *> Note that for COMPZ = 'V', then if N is less than or 00119 *> equal to the minimum divide size, usually 25, then LWORK need 00120 *> only be 1. 00121 *> 00122 *> If LWORK = -1, then a workspace query is assumed; the routine 00123 *> only calculates the optimal sizes of the WORK, RWORK and 00124 *> IWORK arrays, returns these values as the first entries of 00125 *> the WORK, RWORK and IWORK arrays, and no error message 00126 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00127 *> \endverbatim 00128 *> 00129 *> \param[out] RWORK 00130 *> \verbatim 00131 *> RWORK is REAL array, dimension (MAX(1,LRWORK)) 00132 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 00133 *> \endverbatim 00134 *> 00135 *> \param[in] LRWORK 00136 *> \verbatim 00137 *> LRWORK is INTEGER 00138 *> The dimension of the array RWORK. 00139 *> If COMPZ = 'N' or N <= 1, LRWORK must be at least 1. 00140 *> If COMPZ = 'V' and N > 1, LRWORK must be at least 00141 *> 1 + 3*N + 2*N*lg N + 4*N**2 , 00142 *> where lg( N ) = smallest integer k such 00143 *> that 2**k >= N. 00144 *> If COMPZ = 'I' and N > 1, LRWORK must be at least 00145 *> 1 + 4*N + 2*N**2 . 00146 *> Note that for COMPZ = 'I' or 'V', then if N is less than or 00147 *> equal to the minimum divide size, usually 25, then LRWORK 00148 *> need only be max(1,2*(N-1)). 00149 *> 00150 *> If LRWORK = -1, then a workspace query is assumed; the 00151 *> routine only calculates the optimal sizes of the WORK, RWORK 00152 *> and IWORK arrays, returns these values as the first entries 00153 *> of the WORK, RWORK and IWORK arrays, and no error message 00154 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00155 *> \endverbatim 00156 *> 00157 *> \param[out] IWORK 00158 *> \verbatim 00159 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 00160 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00161 *> \endverbatim 00162 *> 00163 *> \param[in] LIWORK 00164 *> \verbatim 00165 *> LIWORK is INTEGER 00166 *> The dimension of the array IWORK. 00167 *> If COMPZ = 'N' or N <= 1, LIWORK must be at least 1. 00168 *> If COMPZ = 'V' or N > 1, LIWORK must be at least 00169 *> 6 + 6*N + 5*N*lg N. 00170 *> If COMPZ = 'I' or N > 1, LIWORK must be at least 00171 *> 3 + 5*N . 00172 *> Note that for COMPZ = 'I' or 'V', then if N is less than or 00173 *> equal to the minimum divide size, usually 25, then LIWORK 00174 *> need only be 1. 00175 *> 00176 *> If LIWORK = -1, then a workspace query is assumed; the 00177 *> routine only calculates the optimal sizes of the WORK, RWORK 00178 *> and IWORK arrays, returns these values as the first entries 00179 *> of the WORK, RWORK and IWORK arrays, and no error message 00180 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00181 *> \endverbatim 00182 *> 00183 *> \param[out] INFO 00184 *> \verbatim 00185 *> INFO is INTEGER 00186 *> = 0: successful exit. 00187 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00188 *> > 0: The algorithm failed to compute an eigenvalue while 00189 *> working on the submatrix lying in rows and columns 00190 *> INFO/(N+1) through mod(INFO,N+1). 00191 *> \endverbatim 00192 * 00193 * Authors: 00194 * ======== 00195 * 00196 *> \author Univ. of Tennessee 00197 *> \author Univ. of California Berkeley 00198 *> \author Univ. of Colorado Denver 00199 *> \author NAG Ltd. 00200 * 00201 *> \date November 2011 00202 * 00203 *> \ingroup complexOTHERcomputational 00204 * 00205 *> \par Contributors: 00206 * ================== 00207 *> 00208 *> Jeff Rutter, Computer Science Division, University of California 00209 *> at Berkeley, USA 00210 * 00211 * ===================================================================== 00212 SUBROUTINE CSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, 00213 $ LRWORK, IWORK, LIWORK, INFO ) 00214 * 00215 * -- LAPACK computational routine (version 3.4.0) -- 00216 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00217 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00218 * November 2011 00219 * 00220 * .. Scalar Arguments .. 00221 CHARACTER COMPZ 00222 INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N 00223 * .. 00224 * .. Array Arguments .. 00225 INTEGER IWORK( * ) 00226 REAL D( * ), E( * ), RWORK( * ) 00227 COMPLEX WORK( * ), Z( LDZ, * ) 00228 * .. 00229 * 00230 * ===================================================================== 00231 * 00232 * .. Parameters .. 00233 REAL ZERO, ONE, TWO 00234 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 00235 * .. 00236 * .. Local Scalars .. 00237 LOGICAL LQUERY 00238 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL, 00239 $ LRWMIN, LWMIN, M, SMLSIZ, START 00240 REAL EPS, ORGNRM, P, TINY 00241 * .. 00242 * .. External Functions .. 00243 LOGICAL LSAME 00244 INTEGER ILAENV 00245 REAL SLAMCH, SLANST 00246 EXTERNAL ILAENV, LSAME, SLAMCH, SLANST 00247 * .. 00248 * .. External Subroutines .. 00249 EXTERNAL XERBLA, CLACPY, CLACRM, CLAED0, CSTEQR, CSWAP, 00250 $ SLASCL, SLASET, SSTEDC, SSTEQR, SSTERF 00251 * .. 00252 * .. Intrinsic Functions .. 00253 INTRINSIC ABS, INT, LOG, MAX, MOD, REAL, SQRT 00254 * .. 00255 * .. Executable Statements .. 00256 * 00257 * Test the input parameters. 00258 * 00259 INFO = 0 00260 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00261 * 00262 IF( LSAME( COMPZ, 'N' ) ) THEN 00263 ICOMPZ = 0 00264 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00265 ICOMPZ = 1 00266 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00267 ICOMPZ = 2 00268 ELSE 00269 ICOMPZ = -1 00270 END IF 00271 IF( ICOMPZ.LT.0 ) THEN 00272 INFO = -1 00273 ELSE IF( N.LT.0 ) THEN 00274 INFO = -2 00275 ELSE IF( ( LDZ.LT.1 ) .OR. 00276 $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN 00277 INFO = -6 00278 END IF 00279 * 00280 IF( INFO.EQ.0 ) THEN 00281 * 00282 * Compute the workspace requirements 00283 * 00284 SMLSIZ = ILAENV( 9, 'CSTEDC', ' ', 0, 0, 0, 0 ) 00285 IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN 00286 LWMIN = 1 00287 LIWMIN = 1 00288 LRWMIN = 1 00289 ELSE IF( N.LE.SMLSIZ ) THEN 00290 LWMIN = 1 00291 LIWMIN = 1 00292 LRWMIN = 2*( N - 1 ) 00293 ELSE IF( ICOMPZ.EQ.1 ) THEN 00294 LGN = INT( LOG( REAL( N ) ) / LOG( TWO ) ) 00295 IF( 2**LGN.LT.N ) 00296 $ LGN = LGN + 1 00297 IF( 2**LGN.LT.N ) 00298 $ LGN = LGN + 1 00299 LWMIN = N*N 00300 LRWMIN = 1 + 3*N + 2*N*LGN + 4*N**2 00301 LIWMIN = 6 + 6*N + 5*N*LGN 00302 ELSE IF( ICOMPZ.EQ.2 ) THEN 00303 LWMIN = 1 00304 LRWMIN = 1 + 4*N + 2*N**2 00305 LIWMIN = 3 + 5*N 00306 END IF 00307 WORK( 1 ) = LWMIN 00308 RWORK( 1 ) = LRWMIN 00309 IWORK( 1 ) = LIWMIN 00310 * 00311 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00312 INFO = -8 00313 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00314 INFO = -10 00315 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00316 INFO = -12 00317 END IF 00318 END IF 00319 * 00320 IF( INFO.NE.0 ) THEN 00321 CALL XERBLA( 'CSTEDC', -INFO ) 00322 RETURN 00323 ELSE IF( LQUERY ) THEN 00324 RETURN 00325 END IF 00326 * 00327 * Quick return if possible 00328 * 00329 IF( N.EQ.0 ) 00330 $ RETURN 00331 IF( N.EQ.1 ) THEN 00332 IF( ICOMPZ.NE.0 ) 00333 $ Z( 1, 1 ) = ONE 00334 RETURN 00335 END IF 00336 * 00337 * If the following conditional clause is removed, then the routine 00338 * will use the Divide and Conquer routine to compute only the 00339 * eigenvalues, which requires (3N + 3N**2) real workspace and 00340 * (2 + 5N + 2N lg(N)) integer workspace. 00341 * Since on many architectures SSTERF is much faster than any other 00342 * algorithm for finding eigenvalues only, it is used here 00343 * as the default. If the conditional clause is removed, then 00344 * information on the size of workspace needs to be changed. 00345 * 00346 * If COMPZ = 'N', use SSTERF to compute the eigenvalues. 00347 * 00348 IF( ICOMPZ.EQ.0 ) THEN 00349 CALL SSTERF( N, D, E, INFO ) 00350 GO TO 70 00351 END IF 00352 * 00353 * If N is smaller than the minimum divide size (SMLSIZ+1), then 00354 * solve the problem with another solver. 00355 * 00356 IF( N.LE.SMLSIZ ) THEN 00357 * 00358 CALL CSTEQR( COMPZ, N, D, E, Z, LDZ, RWORK, INFO ) 00359 * 00360 ELSE 00361 * 00362 * If COMPZ = 'I', we simply call SSTEDC instead. 00363 * 00364 IF( ICOMPZ.EQ.2 ) THEN 00365 CALL SLASET( 'Full', N, N, ZERO, ONE, RWORK, N ) 00366 LL = N*N + 1 00367 CALL SSTEDC( 'I', N, D, E, RWORK, N, 00368 $ RWORK( LL ), LRWORK-LL+1, IWORK, LIWORK, INFO ) 00369 DO 20 J = 1, N 00370 DO 10 I = 1, N 00371 Z( I, J ) = RWORK( ( J-1 )*N+I ) 00372 10 CONTINUE 00373 20 CONTINUE 00374 GO TO 70 00375 END IF 00376 * 00377 * From now on, only option left to be handled is COMPZ = 'V', 00378 * i.e. ICOMPZ = 1. 00379 * 00380 * Scale. 00381 * 00382 ORGNRM = SLANST( 'M', N, D, E ) 00383 IF( ORGNRM.EQ.ZERO ) 00384 $ GO TO 70 00385 * 00386 EPS = SLAMCH( 'Epsilon' ) 00387 * 00388 START = 1 00389 * 00390 * while ( START <= N ) 00391 * 00392 30 CONTINUE 00393 IF( START.LE.N ) THEN 00394 * 00395 * Let FINISH be the position of the next subdiagonal entry 00396 * such that E( FINISH ) <= TINY or FINISH = N if no such 00397 * subdiagonal exists. The matrix identified by the elements 00398 * between START and FINISH constitutes an independent 00399 * sub-problem. 00400 * 00401 FINISH = START 00402 40 CONTINUE 00403 IF( FINISH.LT.N ) THEN 00404 TINY = EPS*SQRT( ABS( D( FINISH ) ) )* 00405 $ SQRT( ABS( D( FINISH+1 ) ) ) 00406 IF( ABS( E( FINISH ) ).GT.TINY ) THEN 00407 FINISH = FINISH + 1 00408 GO TO 40 00409 END IF 00410 END IF 00411 * 00412 * (Sub) Problem determined. Compute its size and solve it. 00413 * 00414 M = FINISH - START + 1 00415 IF( M.GT.SMLSIZ ) THEN 00416 * 00417 * Scale. 00418 * 00419 ORGNRM = SLANST( 'M', M, D( START ), E( START ) ) 00420 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M, 00421 $ INFO ) 00422 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ), 00423 $ M-1, INFO ) 00424 * 00425 CALL CLAED0( N, M, D( START ), E( START ), Z( 1, START ), 00426 $ LDZ, WORK, N, RWORK, IWORK, INFO ) 00427 IF( INFO.GT.0 ) THEN 00428 INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) + 00429 $ MOD( INFO, ( M+1 ) ) + START - 1 00430 GO TO 70 00431 END IF 00432 * 00433 * Scale back. 00434 * 00435 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M, 00436 $ INFO ) 00437 * 00438 ELSE 00439 CALL SSTEQR( 'I', M, D( START ), E( START ), RWORK, M, 00440 $ RWORK( M*M+1 ), INFO ) 00441 CALL CLACRM( N, M, Z( 1, START ), LDZ, RWORK, M, WORK, N, 00442 $ RWORK( M*M+1 ) ) 00443 CALL CLACPY( 'A', N, M, WORK, N, Z( 1, START ), LDZ ) 00444 IF( INFO.GT.0 ) THEN 00445 INFO = START*( N+1 ) + FINISH 00446 GO TO 70 00447 END IF 00448 END IF 00449 * 00450 START = FINISH + 1 00451 GO TO 30 00452 END IF 00453 * 00454 * endwhile 00455 * 00456 * If the problem split any number of times, then the eigenvalues 00457 * will not be properly ordered. Here we permute the eigenvalues 00458 * (and the associated eigenvectors) into ascending order. 00459 * 00460 IF( M.NE.N ) THEN 00461 * 00462 * Use Selection Sort to minimize swaps of eigenvectors 00463 * 00464 DO 60 II = 2, N 00465 I = II - 1 00466 K = I 00467 P = D( I ) 00468 DO 50 J = II, N 00469 IF( D( J ).LT.P ) THEN 00470 K = J 00471 P = D( J ) 00472 END IF 00473 50 CONTINUE 00474 IF( K.NE.I ) THEN 00475 D( K ) = D( I ) 00476 D( I ) = P 00477 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) 00478 END IF 00479 60 CONTINUE 00480 END IF 00481 END IF 00482 * 00483 70 CONTINUE 00484 WORK( 1 ) = LWMIN 00485 RWORK( 1 ) = LRWMIN 00486 IWORK( 1 ) = LIWMIN 00487 * 00488 RETURN 00489 * 00490 * End of CSTEDC 00491 * 00492 END