1 SUBROUTINE PCHEEV( JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, 2 $ DESCZ, WORK, LWORK, RWORK, LRWORK, INFO ) 3* 4* -- ScaLAPACK routine (version 1.7) -- 5* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 6* and University of California, Berkeley. 7* August 14, 2001 8* 9* .. Scalar Arguments .. 10 CHARACTER JOBZ, UPLO 11 INTEGER IA, INFO, IZ, JA, JZ, LRWORK, LWORK, N 12* .. 13* .. Array Arguments .. 14 INTEGER DESCA( * ), DESCZ( * ) 15 REAL RWORK( * ), W( * ) 16 COMPLEX A( * ), WORK( * ), Z( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* PCHEEV computes selected eigenvalues and, optionally, eigenvectors 23* of a real Hermitian matrix A by calling the recommended sequence 24* of ScaLAPACK routines. 25* 26* In its present form, PCHEEV assumes a homogeneous system and makes 27* only spot checks of the consistency of the eigenvalues across the 28* different processes. Because of this, it is possible that a 29* heterogeneous system may return incorrect results without any error 30* messages. 31* 32* Notes 33* ===== 34* A description vector is associated with each 2D block-cyclicly dis- 35* tributed matrix. This vector stores the information required to 36* establish the mapping between a matrix entry and its corresponding 37* process and memory location. 38* 39* In the following comments, the character _ should be read as 40* "of the distributed matrix". Let A be a generic term for any 2D 41* block cyclicly distributed matrix. Its description vector is DESCA: 42* 43* NOTATION STORED IN EXPLANATION 44* --------------- -------------- -------------------------------------- 45* DTYPE_A(global) DESCA( DTYPE_) The descriptor type. 46* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 47* the BLACS process grid A is distribu- 48* ted over. The context itself is glo- 49* bal, but the handle (the integer 50* value) may vary. 51* M_A (global) DESCA( M_ ) The number of rows in the distributed 52* matrix A. 53* N_A (global) DESCA( N_ ) The number of columns in the distri- 54* buted matrix A. 55* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 56* the rows of A. 57* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 58* the columns of A. 59* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 60* row of the matrix A is distributed. 61* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 62* first column of A is distributed. 63* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 64* array storing the local blocks of the 65* distributed matrix A. 66* LLD_A >= MAX(1,LOCr(M_A)). 67* 68* Let K be the number of rows or columns of a distributed matrix, 69* and assume that its process grid has dimension p x q. 70* LOCr( K ) denotes the number of elements of K that a process 71* would receive if K were distributed over the p processes of its 72* process column. 73* Similarly, LOCc( K ) denotes the number of elements of K that a 74* process would receive if K were distributed over the q processes of 75* its process row. 76* The values of LOCr() and LOCc() may be determined via a call to the 77* ScaLAPACK tool function, NUMROC: 78* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 79* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 80* 81* Arguments 82* ========= 83* 84* NP = the number of rows local to a given process. 85* NQ = the number of columns local to a given process. 86* 87* JOBZ (global input) CHARACTER*1 88* Specifies whether or not to compute the eigenvectors: 89* = 'N': Compute eigenvalues only. 90* = 'V': Compute eigenvalues and eigenvectors. 91* 92* UPLO (global input) CHARACTER*1 93* Specifies whether the upper or lower triangular part of the 94* Hermitian matrix A is stored: 95* = 'U': Upper triangular 96* = 'L': Lower triangular 97* 98* N (global input) INTEGER 99* The number of rows and columns of the matrix A. N >= 0. 100* 101* A (local input/workspace) block cyclic COMPLEX array, 102* global dimension (N, N), local dimension ( LLD_A, 103* LOCc(JA+N-1) ) 104* 105* On entry, the Hermitian matrix A. If UPLO = 'U', only the 106* upper triangular part of A is used to define the elements of 107* the Hermitian matrix. If UPLO = 'L', only the lower 108* triangular part of A is used to define the elements of the 109* Hermitian matrix. 110* 111* On exit, the lower triangle (if UPLO='L') or the upper 112* triangle (if UPLO='U') of A, including the diagonal, is 113* destroyed. 114* 115* IA (global input) INTEGER 116* A's global row index, which points to the beginning of the 117* submatrix which is to be operated on. 118* 119* JA (global input) INTEGER 120* A's global column index, which points to the beginning of 121* the submatrix which is to be operated on. 122* 123* DESCA (global and local input) INTEGER array of dimension DLEN_. 124* The array descriptor for the distributed matrix A. 125* If DESCA( CTXT_ ) is incorrect, PCHEEV cannot guarantee 126* correct error reporting. 127* 128* W (global output) REAL array, dimension (N) 129* If INFO=0, the eigenvalues in ascending order. 130* 131* Z (local output) COMPLEX array, 132* global dimension (N, N), 133* local dimension (LLD_Z, LOCc(JZ+N-1)) 134* If JOBZ = 'V', then on normal exit the first M columns of Z 135* contain the orthonormal eigenvectors of the matrix 136* corresponding to the selected eigenvalues. 137* If JOBZ = 'N', then Z is not referenced. 138* 139* IZ (global input) INTEGER 140* Z's global row index, which points to the beginning of the 141* submatrix which is to be operated on. 142* 143* JZ (global input) INTEGER 144* Z's global column index, which points to the beginning of 145* the submatrix which is to be operated on. 146* 147* DESCZ (global and local input) INTEGER array of dimension DLEN_. 148* The array descriptor for the distributed matrix Z. 149* DESCZ( CTXT_ ) must equal DESCA( CTXT_ ) 150* 151* WORK (local workspace/output) COMPLEX array, 152* dimension (LWORK) 153* On output, WORK(1) returns the workspace needed to guarantee 154* completion. If the input parameters are incorrect, WORK(1) 155* may also be incorrect. 156* 157* If JOBZ='N' WORK(1) = minimal workspace for eigenvalues only. 158* If JOBZ='V' WORK(1) = minimal workspace required to 159* generate all the eigenvectors. 160* 161* 162* LWORK (local input) INTEGER 163* See below for definitions of variables used to define LWORK. 164* If no eigenvectors are requested (JOBZ = 'N') then 165* LWORK >= MAX( NB*( NP0+1 ), 3 ) +3*N 166* If eigenvectors are requested (JOBZ = 'V' ) then 167* the amount of workspace required: 168* LWORK >= (NP0 + NQ0 + NB)*NB + 3*N + N^2 169* 170* Variable definitions: 171* NB = DESCA( MB_ ) = DESCA( NB_ ) = 172* DESCZ( MB_ ) = DESCZ( NB_ ) 173* NP0 = NUMROC( NN, NB, 0, 0, NPROW ) 174* NQ0 = NUMROC( MAX( N, NB, 2 ), NB, 0, 0, NPCOL ) 175* 176* If LWORK = -1, the LWORK is global input and a workspace 177* query is assumed; the routine only calculates the minimum 178* size for the WORK array. The required workspace is returned 179* as the first element of WORK and no error message is issued 180* by PXERBLA. 181* 182* RWORK (local workspace/output) COMPLEX array, 183* dimension (LRWORK) 184* On output RWORK(1) returns the 185* REAL workspace needed to 186* guarantee completion. If the input parameters are incorrect, 187* RWORK(1) may also be incorrect. 188* 189* LRWORK (local input) INTEGER 190* Size of RWORK array. 191* If eigenvectors are desired (JOBZ = 'V') then 192* LRWORK >= 2*N + 2*N-2 193* If eigenvectors are not desired (JOBZ = 'N') then 194* LRWORK >= 2*N 195* 196* If LRWORK = -1, the LRWORK is global input and a workspace 197* query is assumed; the routine only calculates the minimum 198* size for the RWORK array. The required workspace is returned 199* as the first element of RWORK and no error message is issued 200* by PXERBLA. 201* 202* INFO (global output) INTEGER 203* = 0: successful exit 204* < 0: If the i-th argument is an array and the j-entry had 205* an illegal value, then INFO = -(i*100+j), if the i-th 206* argument is a scalar and had an illegal value, then 207* INFO = -i. 208* > 0: If INFO = 1 through N, the i(th) eigenvalue did not 209* converge in CSTEQR2 after a total of 30*N iterations. 210* If INFO = N+1, then PCHEEV has detected heterogeneity 211* by finding that eigenvalues were not identical across 212* the process grid. In this case, the accuracy of 213* the results from PCHEEV cannot be guaranteed. 214* 215* Alignment requirements 216* ====================== 217* 218* The distributed submatrices A(IA:*, JA:*) and C(IC:IC+M-1,JC:JC+N-1) 219* must verify some alignment properties, namely the following 220* expressions should be true: 221* 222* ( MB_A.EQ.NB_A.EQ.MB_Z .AND. IROFFA.EQ.IROFFZ .AND. IROFFA.EQ.0 .AND. 223* IAROW.EQ.IZROW ) 224* where 225* IROFFA = MOD( IA-1, MB_A ) and ICOFFA = MOD( JA-1, NB_A ). 226* 227* ===================================================================== 228* 229* Version 1.4 limitations: 230* DESCA(MB_) = DESCA(NB_) 231* DESCA(M_) = DESCZ(M_) 232* DESCA(N_) = DESCZ(N_) 233* DESCA(MB_) = DESCZ(MB_) 234* DESCA(NB_) = DESCZ(NB_) 235* DESCA(RSRC_) = DESCZ(RSRC_) 236* 237* .. Parameters .. 238 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, 239 $ MB_, NB_, RSRC_, CSRC_, LLD_ 240 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 241 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 242 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 243 REAL ZERO, ONE 244 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 245 COMPLEX CZERO, CONE 246 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 247 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 248 INTEGER ITHVAL 249 PARAMETER ( ITHVAL = 10 ) 250* .. 251* .. Local Scalars .. 252 LOGICAL LOWER, WANTZ 253 INTEGER CONTEXTC, CSRC_A, I, IACOL, IAROW, ICOFFA, 254 $ IINFO, INDD, INDE, INDRD, INDRE, INDRWORK, 255 $ INDTAU, INDWORK, INDWORK2, IROFFA, IROFFZ, 256 $ ISCALE, IZROW, J, K, LDC, LLRWORK, LLWORK, 257 $ LRMIN, LRWMIN, LWMIN, MB_A, MB_Z, MYCOL, 258 $ MYPCOLC, MYPROWC, MYROW, NB, NB_A, NB_Z, NP0, 259 $ NPCOL, NPCOLC, NPROCS, NPROW, NPROWC, NQ0, NRC, 260 $ RSIZECSTEQR2, RSRC_A, RSRC_Z, SIZECSTEQR2, 261 $ SIZEPCHETRD, SIZEPCUNMTR 262 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 263 $ SMLNUM 264* .. 265* .. Local Arrays .. 266 INTEGER DESCQR( 10 ), IDUM1( 3 ), IDUM2( 3 ) 267* .. 268* .. External Functions .. 269 LOGICAL LSAME 270 INTEGER INDXG2P, NUMROC, SL_GRIDRESHAPE 271 REAL PCLANHE, PSLAMCH 272 EXTERNAL LSAME, INDXG2P, NUMROC, SL_GRIDRESHAPE, 273 $ PCLANHE, PSLAMCH 274* .. 275* .. External Subroutines .. 276 EXTERNAL BLACS_GRIDEXIT, BLACS_GRIDINFO, CHK1MAT, 277 $ CSTEQR2, DESCINIT, PCELGET, PCGEMR2D, PCHETRD, 278 $ PCHK1MAT, PCHK2MAT, PCLASCL, PCLASET, PCUNMTR, 279 $ PXERBLA, SCOPY, SGAMN2D, SGAMX2D, SSCAL 280* .. 281* .. Intrinsic Functions .. 282 INTRINSIC ABS, CMPLX, ICHAR, INT, MAX, MIN, MOD, REAL, 283 $ SQRT 284* .. 285* .. Executable Statements .. 286* This is just to keep ftnchek and toolpack/1 happy 287 IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_* 288 $ RSRC_.LT.0 )RETURN 289* 290* Quick return 291* 292 IF( N.EQ.0 ) 293 $ RETURN 294* 295* Test the input arguments. 296* 297 CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL ) 298 INFO = 0 299* 300* Initialize pointer to some safe value 301* 302 INDTAU = 1 303 INDD = 1 304 INDE = 1 305 INDWORK = 1 306 INDWORK2 = 1 307* 308 INDRE = 1 309 INDRD = 1 310 INDRWORK = 1 311* 312 WANTZ = LSAME( JOBZ, 'V' ) 313 IF( NPROW.EQ.-1 ) THEN 314 INFO = -( 700+CTXT_ ) 315 ELSE IF( WANTZ ) THEN 316 IF( DESCA( CTXT_ ).NE.DESCZ( CTXT_ ) ) THEN 317 INFO = -( 1200+CTXT_ ) 318 END IF 319 END IF 320 IF( INFO.EQ.0 ) THEN 321 CALL CHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, INFO ) 322 IF( WANTZ ) 323 $ CALL CHK1MAT( N, 3, N, 3, IZ, JZ, DESCZ, 12, INFO ) 324* 325 IF( INFO.EQ.0 ) THEN 326* 327* Get machine constants. 328* 329 SAFMIN = PSLAMCH( DESCA( CTXT_ ), 'Safe minimum' ) 330 EPS = PSLAMCH( DESCA( CTXT_ ), 'Precision' ) 331 SMLNUM = SAFMIN / EPS 332 BIGNUM = ONE / SMLNUM 333 RMIN = SQRT( SMLNUM ) 334 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 335* 336 NPROCS = NPROW*NPCOL 337 NB_A = DESCA( NB_ ) 338 MB_A = DESCA( MB_ ) 339 NB = NB_A 340 LOWER = LSAME( UPLO, 'L' ) 341* 342 RSRC_A = DESCA( RSRC_ ) 343 CSRC_A = DESCA( CSRC_ ) 344 IROFFA = MOD( IA-1, MB_A ) 345 ICOFFA = MOD( JA-1, NB_A ) 346 IAROW = INDXG2P( 1, NB_A, MYROW, RSRC_A, NPROW ) 347 IACOL = INDXG2P( 1, MB_A, MYCOL, CSRC_A, NPCOL ) 348 NP0 = NUMROC( N+IROFFA, NB, MYROW, IAROW, NPROW ) 349 NQ0 = NUMROC( N+ICOFFA, NB, MYCOL, IACOL, NPCOL ) 350 IF( WANTZ ) THEN 351 NB_Z = DESCZ( NB_ ) 352 MB_Z = DESCZ( MB_ ) 353 RSRC_Z = DESCZ( RSRC_ ) 354 IROFFZ = MOD( IZ-1, MB_A ) 355 IZROW = INDXG2P( 1, NB_A, MYROW, RSRC_Z, NPROW ) 356 ELSE 357 IROFFZ = 0 358 IZROW = 0 359 END IF 360* 361* COMPLEX work space for PCHETRD 362* 363 CALL PCHETRD( UPLO, N, A, IA, JA, DESCA, RWORK( INDD ), 364 $ RWORK( INDE ), WORK( INDTAU ), 365 $ WORK( INDWORK ), -1, IINFO ) 366 SIZEPCHETRD = INT( ABS( WORK( 1 ) ) ) 367* 368* COMPLEX work space for PCUNMTR 369* 370 IF( WANTZ ) THEN 371 CALL PCUNMTR( 'L', UPLO, 'N', N, N, A, IA, JA, DESCA, 372 $ WORK( INDTAU ), Z, IZ, JZ, DESCZ, 373 $ WORK( INDWORK ), -1, IINFO ) 374 SIZEPCUNMTR = INT( ABS( WORK( 1 ) ) ) 375 ELSE 376 SIZEPCUNMTR = 0 377 END IF 378* 379* REAL work space for CSTEQR2 380* 381 IF( WANTZ ) THEN 382 RSIZECSTEQR2 = MIN( 1, 2*N-2 ) 383 ELSE 384 RSIZECSTEQR2 = 0 385 END IF 386* 387* Initialize the context of the single column distributed 388* matrix required by CSTEQR2. This specific distribution 389* allows each process to do 1/pth of the work updating matrix 390* Q during CSTEQR2 and achieve some parallelization to an 391* otherwise serial subroutine. 392* 393 LDC = 0 394 IF( WANTZ ) THEN 395 CONTEXTC = SL_GRIDRESHAPE( DESCA( CTXT_ ), 0, 1, 1, 396 $ NPROCS, 1 ) 397 CALL BLACS_GRIDINFO( CONTEXTC, NPROWC, NPCOLC, MYPROWC, 398 $ MYPCOLC ) 399 NRC = NUMROC( N, NB_A, MYPROWC, 0, NPROCS ) 400 LDC = MAX( 1, NRC ) 401 CALL DESCINIT( DESCQR, N, N, NB, NB, 0, 0, CONTEXTC, LDC, 402 $ INFO ) 403 END IF 404* 405* COMPLEX work space for CSTEQR2 406* 407 IF( WANTZ ) THEN 408 SIZECSTEQR2 = N*LDC 409 ELSE 410 SIZECSTEQR2 = 0 411 END IF 412* 413* Set up pointers into the WORK array 414* 415 INDTAU = 1 416 INDD = INDTAU + N 417 INDE = INDD + N 418 INDWORK = INDE + N 419 INDWORK2 = INDWORK + N*LDC 420 LLWORK = LWORK - INDWORK + 1 421* 422* Set up pointers into the RWORK array 423* 424 INDRE = 1 425 INDRD = INDRE + N 426 INDRWORK = INDRD + N 427 LLRWORK = LRWORK - INDRWORK + 1 428* 429* Compute the total amount of space needed 430* 431 LRWMIN = 2*N + RSIZECSTEQR2 432 LWMIN = 3*N + MAX( SIZEPCHETRD, SIZEPCUNMTR, SIZECSTEQR2 ) 433* 434 END IF 435 IF( INFO.EQ.0 ) THEN 436 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 437 INFO = -1 438 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 439 INFO = -2 440 ELSE IF( LWORK.LT.LWMIN .AND. LWORK.NE.-1 ) THEN 441 INFO = -14 442 ELSE IF( LRWORK.LT.LRWMIN .AND. LRWORK.NE.-1 ) THEN 443 INFO = -16 444 ELSE IF( IROFFA.NE.0 ) THEN 445 INFO = -5 446 ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN 447 INFO = -( 700+NB_ ) 448 END IF 449 IF( WANTZ ) THEN 450 IF( IROFFA.NE.IROFFZ ) THEN 451 INFO = -10 452 ELSE IF( IAROW.NE.IZROW ) THEN 453 INFO = -10 454 ELSE IF( DESCA( M_ ).NE.DESCZ( M_ ) ) THEN 455 INFO = -( 1200+M_ ) 456 ELSE IF( DESCA( N_ ).NE.DESCZ( N_ ) ) THEN 457 INFO = -( 1200+N_ ) 458 ELSE IF( DESCA( MB_ ).NE.DESCZ( MB_ ) ) THEN 459 INFO = -( 1200+MB_ ) 460 ELSE IF( DESCA( NB_ ).NE.DESCZ( NB_ ) ) THEN 461 INFO = -( 1200+NB_ ) 462 ELSE IF( DESCA( RSRC_ ).NE.DESCZ( RSRC_ ) ) THEN 463 INFO = -( 1200+RSRC_ ) 464 ELSE IF( DESCA( CTXT_ ).NE.DESCZ( CTXT_ ) ) THEN 465 INFO = -( 1200+CTXT_ ) 466 END IF 467 END IF 468 END IF 469 IF( WANTZ ) THEN 470 IDUM1( 1 ) = ICHAR( 'V' ) 471 ELSE 472 IDUM1( 1 ) = ICHAR( 'N' ) 473 END IF 474 IDUM2( 1 ) = 1 475 IF( LOWER ) THEN 476 IDUM1( 2 ) = ICHAR( 'L' ) 477 ELSE 478 IDUM1( 2 ) = ICHAR( 'U' ) 479 END IF 480 IDUM2( 2 ) = 2 481 IF( LWORK.EQ.-1 ) THEN 482 IDUM1( 3 ) = -1 483 ELSE 484 IDUM1( 3 ) = 1 485 END IF 486 IDUM2( 3 ) = 3 487 IF( WANTZ ) THEN 488 CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 7, N, 3, N, 3, IZ, 489 $ JZ, DESCZ, 12, 3, IDUM1, IDUM2, INFO ) 490 ELSE 491 CALL PCHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, 3, IDUM1, 492 $ IDUM2, INFO ) 493 END IF 494 WORK( 1 ) = CMPLX( LWMIN ) 495 RWORK( 1 ) = REAL( LRWMIN ) 496 END IF 497* 498 IF( INFO.NE.0 ) THEN 499 CALL PXERBLA( DESCA( CTXT_ ), 'PCHEEV', -INFO ) 500 IF( WANTZ ) 501 $ CALL BLACS_GRIDEXIT( CONTEXTC ) 502 RETURN 503 ELSE IF( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 ) THEN 504 IF( WANTZ ) 505 $ CALL BLACS_GRIDEXIT( CONTEXTC ) 506 RETURN 507 END IF 508* 509* Scale matrix to allowable range, if necessary. 510* 511 ISCALE = 0 512* 513 ANRM = PCLANHE( 'M', UPLO, N, A, IA, JA, DESCA, 514 $ RWORK( INDRWORK ) ) 515* 516* 517 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 518 ISCALE = 1 519 SIGMA = RMIN / ANRM 520 ELSE IF( ANRM.GT.RMAX ) THEN 521 ISCALE = 1 522 SIGMA = RMAX / ANRM 523 END IF 524* 525 IF( ISCALE.EQ.1 ) THEN 526 CALL PCLASCL( UPLO, ONE, SIGMA, N, N, A, IA, JA, DESCA, IINFO ) 527 END IF 528* 529* Reduce Hermitian matrix to tridiagonal form. 530* 531 CALL PCHETRD( UPLO, N, A, IA, JA, DESCA, RWORK( INDRD ), 532 $ RWORK( INDRE ), WORK( INDTAU ), WORK( INDWORK ), 533 $ LLWORK, IINFO ) 534* 535* Copy the values of D, E to all processes. 536* 537 DO 10 I = 1, N 538 CALL PCELGET( 'A', ' ', WORK( INDD+I-1 ), A, I+IA-1, I+JA-1, 539 $ DESCA ) 540 RWORK( INDRD+I-1 ) = REAL( WORK( INDD+I-1 ) ) 541 10 CONTINUE 542 IF( LSAME( UPLO, 'U' ) ) THEN 543 DO 20 I = 1, N - 1 544 CALL PCELGET( 'A', ' ', WORK( INDE+I-1 ), A, I+IA-1, I+JA, 545 $ DESCA ) 546 RWORK( INDRE+I-1 ) = REAL( WORK( INDE+I-1 ) ) 547 20 CONTINUE 548 ELSE 549 DO 30 I = 1, N - 1 550 CALL PCELGET( 'A', ' ', WORK( INDE+I-1 ), A, I+IA, I+JA-1, 551 $ DESCA ) 552 RWORK( INDRE+I-1 ) = REAL( WORK( INDE+I-1 ) ) 553 30 CONTINUE 554 END IF 555* 556 IF( WANTZ ) THEN 557* 558 CALL PCLASET( 'Full', N, N, CZERO, CONE, WORK( INDWORK ), 1, 1, 559 $ DESCQR ) 560* 561* CSTEQR2 is a modified version of LAPACK's CSTEQR. The 562* modifications allow each process to perform partial updates 563* to matrix Q. 564* 565 CALL CSTEQR2( 'I', N, RWORK( INDRD ), RWORK( INDRE ), 566 $ WORK( INDWORK ), LDC, NRC, RWORK( INDRWORK ), 567 $ INFO ) 568* 569 CALL PCGEMR2D( N, N, WORK( INDWORK ), 1, 1, DESCQR, Z, IA, JA, 570 $ DESCZ, CONTEXTC ) 571* 572 CALL PCUNMTR( 'L', UPLO, 'N', N, N, A, IA, JA, DESCA, 573 $ WORK( INDTAU ), Z, IZ, JZ, DESCZ, 574 $ WORK( INDWORK ), LLWORK, IINFO ) 575* 576 ELSE 577* 578 CALL CSTEQR2( 'N', N, RWORK( INDRD ), RWORK( INDRE ), 579 $ WORK( INDWORK ), 1, 1, RWORK( INDRWORK ), INFO ) 580 END IF 581* 582* Copy eigenvalues from workspace to output array 583* 584 CALL SCOPY( N, RWORK( INDD ), 1, W, 1 ) 585* 586* If matrix was scaled, then rescale eigenvalues appropriately. 587* 588 IF( ISCALE.EQ.1 ) THEN 589 CALL SSCAL( N, ONE / SIGMA, W, 1 ) 590 END IF 591* 592 WORK( 1 ) = REAL( LWMIN ) 593* 594* Free up resources 595* 596 IF( WANTZ ) THEN 597 CALL BLACS_GRIDEXIT( CONTEXTC ) 598 END IF 599* 600* Compare every ith eigenvalue, or all if there are only a few, 601* across the process grid to check for heterogeneity. 602* 603 IF( N.LE.ITHVAL ) THEN 604 J = N 605 K = 1 606 ELSE 607 J = N / ITHVAL 608 K = ITHVAL 609 END IF 610* 611 LRMIN = INT( RWORK( 1 ) ) 612 INDTAU = 0 613 INDE = INDTAU + J 614 DO 40 I = 1, J 615 RWORK( I+INDTAU ) = W( ( I-1 )*K+1 ) 616 RWORK( I+INDE ) = W( ( I-1 )*K+1 ) 617 40 CONTINUE 618* 619 CALL SGAMN2D( DESCA( CTXT_ ), 'All', ' ', J, 1, RWORK( 1+INDTAU ), 620 $ J, 1, 1, -1, -1, 0 ) 621 CALL SGAMX2D( DESCA( CTXT_ ), 'All', ' ', J, 1, RWORK( 1+INDE ), 622 $ J, 1, 1, -1, -1, 0 ) 623* 624 DO 50 I = 1, J 625 IF( INFO.EQ.0 .AND. ( RWORK( I+INDTAU )-RWORK( I+INDE ).NE. 626 $ ZERO ) ) THEN 627 INFO = N + 1 628 END IF 629 50 CONTINUE 630 RWORK( 1 ) = LRMIN 631* 632 RETURN 633* 634* End of PCHEEV 635* 636 END 637