1*> \brief \b ZSTEDC 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZSTEDC + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstedc.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstedc.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstedc.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, 22* LRWORK, IWORK, LIWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER COMPZ 26* INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N 27* .. 28* .. Array Arguments .. 29* INTEGER IWORK( * ) 30* DOUBLE PRECISION D( * ), E( * ), RWORK( * ) 31* COMPLEX*16 WORK( * ), Z( LDZ, * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> ZSTEDC computes all eigenvalues and, optionally, eigenvectors of a 41*> symmetric tridiagonal matrix using the divide and conquer method. 42*> The eigenvectors of a full or band complex Hermitian matrix can also 43*> be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this 44*> matrix to tridiagonal form. 45*> 46*> This code makes very mild assumptions about floating point 47*> arithmetic. It will work on machines with a guard digit in 48*> add/subtract, or on those binary machines without guard digits 49*> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. 50*> It could conceivably fail on hexadecimal or decimal machines 51*> without guard digits, but we know of none. See DLAED3 for details. 52*> \endverbatim 53* 54* Arguments: 55* ========== 56* 57*> \param[in] COMPZ 58*> \verbatim 59*> COMPZ is CHARACTER*1 60*> = 'N': Compute eigenvalues only. 61*> = 'I': Compute eigenvectors of tridiagonal matrix also. 62*> = 'V': Compute eigenvectors of original Hermitian matrix 63*> also. On entry, Z contains the unitary matrix used 64*> to reduce the original matrix to tridiagonal form. 65*> \endverbatim 66*> 67*> \param[in] N 68*> \verbatim 69*> N is INTEGER 70*> The dimension of the symmetric tridiagonal matrix. N >= 0. 71*> \endverbatim 72*> 73*> \param[in,out] D 74*> \verbatim 75*> D is DOUBLE PRECISION array, dimension (N) 76*> On entry, the diagonal elements of the tridiagonal matrix. 77*> On exit, if INFO = 0, the eigenvalues in ascending order. 78*> \endverbatim 79*> 80*> \param[in,out] E 81*> \verbatim 82*> E is DOUBLE PRECISION array, dimension (N-1) 83*> On entry, the subdiagonal elements of the tridiagonal matrix. 84*> On exit, E has been destroyed. 85*> \endverbatim 86*> 87*> \param[in,out] Z 88*> \verbatim 89*> Z is COMPLEX*16 array, dimension (LDZ,N) 90*> On entry, if COMPZ = 'V', then Z contains the unitary 91*> matrix used in the reduction to tridiagonal form. 92*> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the 93*> orthonormal eigenvectors of the original Hermitian matrix, 94*> and if COMPZ = 'I', Z contains the orthonormal eigenvectors 95*> of the symmetric tridiagonal matrix. 96*> If COMPZ = 'N', then Z is not referenced. 97*> \endverbatim 98*> 99*> \param[in] LDZ 100*> \verbatim 101*> LDZ is INTEGER 102*> The leading dimension of the array Z. LDZ >= 1. 103*> If eigenvectors are desired, then LDZ >= max(1,N). 104*> \endverbatim 105*> 106*> \param[out] WORK 107*> \verbatim 108*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 109*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 110*> \endverbatim 111*> 112*> \param[in] LWORK 113*> \verbatim 114*> LWORK is INTEGER 115*> The dimension of the array WORK. 116*> If COMPZ = 'N' or 'I', or N <= 1, LWORK must be at least 1. 117*> If COMPZ = 'V' and N > 1, LWORK must be at least N*N. 118*> Note that for COMPZ = 'V', then if N is less than or 119*> equal to the minimum divide size, usually 25, then LWORK need 120*> only be 1. 121*> 122*> If LWORK = -1, then a workspace query is assumed; the routine 123*> only calculates the optimal sizes of the WORK, RWORK and 124*> IWORK arrays, returns these values as the first entries of 125*> the WORK, RWORK and IWORK arrays, and no error message 126*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 127*> \endverbatim 128*> 129*> \param[out] RWORK 130*> \verbatim 131*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) 132*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 133*> \endverbatim 134*> 135*> \param[in] LRWORK 136*> \verbatim 137*> LRWORK is INTEGER 138*> The dimension of the array RWORK. 139*> If COMPZ = 'N' or N <= 1, LRWORK must be at least 1. 140*> If COMPZ = 'V' and N > 1, LRWORK must be at least 141*> 1 + 3*N + 2*N*lg N + 4*N**2 , 142*> where lg( N ) = smallest integer k such 143*> that 2**k >= N. 144*> If COMPZ = 'I' and N > 1, LRWORK must be at least 145*> 1 + 4*N + 2*N**2 . 146*> Note that for COMPZ = 'I' or 'V', then if N is less than or 147*> equal to the minimum divide size, usually 25, then LRWORK 148*> need only be max(1,2*(N-1)). 149*> 150*> If LRWORK = -1, then a workspace query is assumed; the 151*> routine only calculates the optimal sizes of the WORK, RWORK 152*> and IWORK arrays, returns these values as the first entries 153*> of the WORK, RWORK and IWORK arrays, and no error message 154*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 155*> \endverbatim 156*> 157*> \param[out] IWORK 158*> \verbatim 159*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 160*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 161*> \endverbatim 162*> 163*> \param[in] LIWORK 164*> \verbatim 165*> LIWORK is INTEGER 166*> The dimension of the array IWORK. 167*> If COMPZ = 'N' or N <= 1, LIWORK must be at least 1. 168*> If COMPZ = 'V' or N > 1, LIWORK must be at least 169*> 6 + 6*N + 5*N*lg N. 170*> If COMPZ = 'I' or N > 1, LIWORK must be at least 171*> 3 + 5*N . 172*> Note that for COMPZ = 'I' or 'V', then if N is less than or 173*> equal to the minimum divide size, usually 25, then LIWORK 174*> need only be 1. 175*> 176*> If LIWORK = -1, then a workspace query is assumed; the 177*> routine only calculates the optimal sizes of the WORK, RWORK 178*> and IWORK arrays, returns these values as the first entries 179*> of the WORK, RWORK and IWORK arrays, and no error message 180*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 181*> \endverbatim 182*> 183*> \param[out] INFO 184*> \verbatim 185*> INFO is INTEGER 186*> = 0: successful exit. 187*> < 0: if INFO = -i, the i-th argument had an illegal value. 188*> > 0: The algorithm failed to compute an eigenvalue while 189*> working on the submatrix lying in rows and columns 190*> INFO/(N+1) through mod(INFO,N+1). 191*> \endverbatim 192* 193* Authors: 194* ======== 195* 196*> \author Univ. of Tennessee 197*> \author Univ. of California Berkeley 198*> \author Univ. of Colorado Denver 199*> \author NAG Ltd. 200* 201*> \ingroup complex16OTHERcomputational 202* 203*> \par Contributors: 204* ================== 205*> 206*> Jeff Rutter, Computer Science Division, University of California 207*> at Berkeley, USA 208* 209* ===================================================================== 210 SUBROUTINE ZSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, 211 $ LRWORK, IWORK, LIWORK, INFO ) 212* 213* -- LAPACK computational routine -- 214* -- LAPACK is a software package provided by Univ. of Tennessee, -- 215* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 216* 217* .. Scalar Arguments .. 218 CHARACTER COMPZ 219 INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N 220* .. 221* .. Array Arguments .. 222 INTEGER IWORK( * ) 223 DOUBLE PRECISION D( * ), E( * ), RWORK( * ) 224 COMPLEX*16 WORK( * ), Z( LDZ, * ) 225* .. 226* 227* ===================================================================== 228* 229* .. Parameters .. 230 DOUBLE PRECISION ZERO, ONE, TWO 231 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 232* .. 233* .. Local Scalars .. 234 LOGICAL LQUERY 235 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL, 236 $ LRWMIN, LWMIN, M, SMLSIZ, START 237 DOUBLE PRECISION EPS, ORGNRM, P, TINY 238* .. 239* .. External Functions .. 240 LOGICAL LSAME 241 INTEGER ILAENV 242 DOUBLE PRECISION DLAMCH, DLANST 243 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST 244* .. 245* .. External Subroutines .. 246 EXTERNAL DLASCL, DLASET, DSTEDC, DSTEQR, DSTERF, XERBLA, 247 $ ZLACPY, ZLACRM, ZLAED0, ZSTEQR, ZSWAP 248* .. 249* .. Intrinsic Functions .. 250 INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT 251* .. 252* .. Executable Statements .. 253* 254* Test the input parameters. 255* 256 INFO = 0 257 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 258* 259 IF( LSAME( COMPZ, 'N' ) ) THEN 260 ICOMPZ = 0 261 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 262 ICOMPZ = 1 263 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 264 ICOMPZ = 2 265 ELSE 266 ICOMPZ = -1 267 END IF 268 IF( ICOMPZ.LT.0 ) THEN 269 INFO = -1 270 ELSE IF( N.LT.0 ) THEN 271 INFO = -2 272 ELSE IF( ( LDZ.LT.1 ) .OR. 273 $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN 274 INFO = -6 275 END IF 276* 277 IF( INFO.EQ.0 ) THEN 278* 279* Compute the workspace requirements 280* 281 SMLSIZ = ILAENV( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 ) 282 IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN 283 LWMIN = 1 284 LIWMIN = 1 285 LRWMIN = 1 286 ELSE IF( N.LE.SMLSIZ ) THEN 287 LWMIN = 1 288 LIWMIN = 1 289 LRWMIN = 2*( N - 1 ) 290 ELSE IF( ICOMPZ.EQ.1 ) THEN 291 LGN = INT( LOG( DBLE( N ) ) / LOG( TWO ) ) 292 IF( 2**LGN.LT.N ) 293 $ LGN = LGN + 1 294 IF( 2**LGN.LT.N ) 295 $ LGN = LGN + 1 296 LWMIN = N*N 297 LRWMIN = 1 + 3*N + 2*N*LGN + 4*N**2 298 LIWMIN = 6 + 6*N + 5*N*LGN 299 ELSE IF( ICOMPZ.EQ.2 ) THEN 300 LWMIN = 1 301 LRWMIN = 1 + 4*N + 2*N**2 302 LIWMIN = 3 + 5*N 303 END IF 304 WORK( 1 ) = LWMIN 305 RWORK( 1 ) = LRWMIN 306 IWORK( 1 ) = LIWMIN 307* 308 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 309 INFO = -8 310 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 311 INFO = -10 312 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 313 INFO = -12 314 END IF 315 END IF 316* 317 IF( INFO.NE.0 ) THEN 318 CALL XERBLA( 'ZSTEDC', -INFO ) 319 RETURN 320 ELSE IF( LQUERY ) THEN 321 RETURN 322 END IF 323* 324* Quick return if possible 325* 326 IF( N.EQ.0 ) 327 $ RETURN 328 IF( N.EQ.1 ) THEN 329 IF( ICOMPZ.NE.0 ) 330 $ Z( 1, 1 ) = ONE 331 RETURN 332 END IF 333* 334* If the following conditional clause is removed, then the routine 335* will use the Divide and Conquer routine to compute only the 336* eigenvalues, which requires (3N + 3N**2) real workspace and 337* (2 + 5N + 2N lg(N)) integer workspace. 338* Since on many architectures DSTERF is much faster than any other 339* algorithm for finding eigenvalues only, it is used here 340* as the default. If the conditional clause is removed, then 341* information on the size of workspace needs to be changed. 342* 343* If COMPZ = 'N', use DSTERF to compute the eigenvalues. 344* 345 IF( ICOMPZ.EQ.0 ) THEN 346 CALL DSTERF( N, D, E, INFO ) 347 GO TO 70 348 END IF 349* 350* If N is smaller than the minimum divide size (SMLSIZ+1), then 351* solve the problem with another solver. 352* 353 IF( N.LE.SMLSIZ ) THEN 354* 355 CALL ZSTEQR( COMPZ, N, D, E, Z, LDZ, RWORK, INFO ) 356* 357 ELSE 358* 359* If COMPZ = 'I', we simply call DSTEDC instead. 360* 361 IF( ICOMPZ.EQ.2 ) THEN 362 CALL DLASET( 'Full', N, N, ZERO, ONE, RWORK, N ) 363 LL = N*N + 1 364 CALL DSTEDC( 'I', N, D, E, RWORK, N, 365 $ RWORK( LL ), LRWORK-LL+1, IWORK, LIWORK, INFO ) 366 DO 20 J = 1, N 367 DO 10 I = 1, N 368 Z( I, J ) = RWORK( ( J-1 )*N+I ) 369 10 CONTINUE 370 20 CONTINUE 371 GO TO 70 372 END IF 373* 374* From now on, only option left to be handled is COMPZ = 'V', 375* i.e. ICOMPZ = 1. 376* 377* Scale. 378* 379 ORGNRM = DLANST( 'M', N, D, E ) 380 IF( ORGNRM.EQ.ZERO ) 381 $ GO TO 70 382* 383 EPS = DLAMCH( 'Epsilon' ) 384* 385 START = 1 386* 387* while ( START <= N ) 388* 389 30 CONTINUE 390 IF( START.LE.N ) THEN 391* 392* Let FINISH be the position of the next subdiagonal entry 393* such that E( FINISH ) <= TINY or FINISH = N if no such 394* subdiagonal exists. The matrix identified by the elements 395* between START and FINISH constitutes an independent 396* sub-problem. 397* 398 FINISH = START 399 40 CONTINUE 400 IF( FINISH.LT.N ) THEN 401 TINY = EPS*SQRT( ABS( D( FINISH ) ) )* 402 $ SQRT( ABS( D( FINISH+1 ) ) ) 403 IF( ABS( E( FINISH ) ).GT.TINY ) THEN 404 FINISH = FINISH + 1 405 GO TO 40 406 END IF 407 END IF 408* 409* (Sub) Problem determined. Compute its size and solve it. 410* 411 M = FINISH - START + 1 412 IF( M.GT.SMLSIZ ) THEN 413* 414* Scale. 415* 416 ORGNRM = DLANST( 'M', M, D( START ), E( START ) ) 417 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M, 418 $ INFO ) 419 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ), 420 $ M-1, INFO ) 421* 422 CALL ZLAED0( N, M, D( START ), E( START ), Z( 1, START ), 423 $ LDZ, WORK, N, RWORK, IWORK, INFO ) 424 IF( INFO.GT.0 ) THEN 425 INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) + 426 $ MOD( INFO, ( M+1 ) ) + START - 1 427 GO TO 70 428 END IF 429* 430* Scale back. 431* 432 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M, 433 $ INFO ) 434* 435 ELSE 436 CALL DSTEQR( 'I', M, D( START ), E( START ), RWORK, M, 437 $ RWORK( M*M+1 ), INFO ) 438 CALL ZLACRM( N, M, Z( 1, START ), LDZ, RWORK, M, WORK, N, 439 $ RWORK( M*M+1 ) ) 440 CALL ZLACPY( 'A', N, M, WORK, N, Z( 1, START ), LDZ ) 441 IF( INFO.GT.0 ) THEN 442 INFO = START*( N+1 ) + FINISH 443 GO TO 70 444 END IF 445 END IF 446* 447 START = FINISH + 1 448 GO TO 30 449 END IF 450* 451* endwhile 452* 453* 454* Use Selection Sort to minimize swaps of eigenvectors 455* 456 DO 60 II = 2, N 457 I = II - 1 458 K = I 459 P = D( I ) 460 DO 50 J = II, N 461 IF( D( J ).LT.P ) THEN 462 K = J 463 P = D( J ) 464 END IF 465 50 CONTINUE 466 IF( K.NE.I ) THEN 467 D( K ) = D( I ) 468 D( I ) = P 469 CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) 470 END IF 471 60 CONTINUE 472 END IF 473* 474 70 CONTINUE 475 WORK( 1 ) = LWMIN 476 RWORK( 1 ) = LRWMIN 477 IWORK( 1 ) = LIWMIN 478* 479 RETURN 480* 481* End of ZSTEDC 482* 483 END 484