1 SUBROUTINE PDSYEVR( JOBZ, RANGE, UPLO, N, A, IA, JA, 2 $ DESCA, VL, VU, IL, IU, M, NZ, W, Z, IZ, 3 $ JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, 4 $ INFO ) 5 6 IMPLICIT NONE 7* 8* -- ScaLAPACK routine (version 2.0.2) -- 9* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver 10* May 1 2012 11* 12* .. Scalar Arguments .. 13 CHARACTER JOBZ, RANGE, UPLO 14 INTEGER IA, IL, INFO, IU, IZ, JA, JZ, LIWORK, LWORK, M, 15 $ N, NZ 16 DOUBLE PRECISION VL, VU 17* .. 18* .. Array Arguments .. 19 INTEGER DESCA( * ), DESCZ( * ), IWORK( * ) 20 DOUBLE PRECISION A( * ), W( * ), WORK( * ), Z( * ) 21* .. 22* 23* Purpose 24* ======= 25* 26* PDSYEVR computes selected eigenvalues and, optionally, eigenvectors 27* of a real symmetric matrix A distributed in 2D blockcyclic format 28* by calling the recommended sequence of ScaLAPACK routines. 29* 30* First, the matrix A is reduced to real symmetric tridiagonal form. 31* Then, the eigenproblem is solved using the parallel MRRR algorithm. 32* Last, if eigenvectors have been computed, a backtransformation is done. 33* 34* Upon successful completion, each processor stores a copy of all computed 35* eigenvalues in W. The eigenvector matrix Z is stored in 36* 2D blockcyclic format distributed over all processors. 37* 38* Note that subsets of eigenvalues/vectors can be selected by 39* specifying a range of values or a range of indices for the desired 40* eigenvalues. 41* 42* For constructive feedback and comments, please contact cvoemel@lbl.gov 43* C. Voemel 44* 45* Arguments 46* ========= 47* 48* JOBZ (global input) CHARACTER*1 49* Specifies whether or not to compute the eigenvectors: 50* = 'N': Compute eigenvalues only. 51* = 'V': Compute eigenvalues and eigenvectors. 52* 53* RANGE (global input) CHARACTER*1 54* = 'A': all eigenvalues will be found. 55* = 'V': all eigenvalues in the interval [VL,VU] will be found. 56* = 'I': the IL-th through IU-th eigenvalues will be found. 57* 58* UPLO (global input) CHARACTER*1 59* Specifies whether the upper or lower triangular part of the 60* symmetric matrix A is stored: 61* = 'U': Upper triangular 62* = 'L': Lower triangular 63* 64* N (global input) INTEGER 65* The number of rows and columns of the matrix A. N >= 0 66* 67* A (local input/workspace) 2D block cyclic DOUBLE PRECISION array, 68* global dimension (N, N), 69* local dimension ( LLD_A, LOCc(JA+N-1) ), 70* (see Notes below for more detailed explanation of 2d arrays) 71* 72* On entry, the symmetric matrix A. If UPLO = 'U', only the 73* upper triangular part of A is used to define the elements of 74* the symmetric matrix. If UPLO = 'L', only the lower 75* triangular part of A is used to define the elements of the 76* symmetric matrix. 77* 78* On exit, the lower triangle (if UPLO='L') or the upper 79* triangle (if UPLO='U') of A, including the diagonal, is 80* destroyed. 81* 82* IA (global input) INTEGER 83* A's global row index, which points to the beginning of the 84* submatrix which is to be operated on. 85* It should be set to 1 when operating on a full matrix. 86* 87* JA (global input) INTEGER 88* A's global column index, which points to the beginning of 89* the submatrix which is to be operated on. 90* It should be set to 1 when operating on a full matrix. 91* 92* DESCA (global and local input) INTEGER array of dimension DLEN=9. 93* The array descriptor for the distributed matrix A. 94* The descriptor stores details about the 2D block-cyclic 95* storage, see the notes below. 96* If DESCA is incorrect, PDSYEVR cannot guarantee 97* correct error reporting. 98* Also note the array alignment requirements specified below. 99* 100* VL (global input) DOUBLE PRECISION 101* If RANGE='V', the lower bound of the interval to be searched 102* for eigenvalues. Not referenced if RANGE = 'A' or 'I'. 103* 104* VU (global input) DOUBLE PRECISION 105* If RANGE='V', the upper bound of the interval to be searched 106* for eigenvalues. Not referenced if RANGE = 'A' or 'I'. 107* 108* IL (global input) INTEGER 109* If RANGE='I', the index (from smallest to largest) of the 110* smallest eigenvalue to be returned. IL >= 1. 111* Not referenced if RANGE = 'A'. 112* 113* IU (global input) INTEGER 114* If RANGE='I', the index (from smallest to largest) of the 115* largest eigenvalue to be returned. min(IL,N) <= IU <= N. 116* Not referenced if RANGE = 'A'. 117* 118* M (global output) INTEGER 119* Total number of eigenvalues found. 0 <= M <= N. 120* 121* NZ (global output) INTEGER 122* Total number of eigenvectors computed. 0 <= NZ <= M. 123* The number of columns of Z that are filled. 124* If JOBZ .NE. 'V', NZ is not referenced. 125* If JOBZ .EQ. 'V', NZ = M 126* 127* W (global output) DOUBLE PRECISION array, dimension (N) 128* Upon successful exit, the first M entries contain the selected 129* eigenvalues in ascending order. 130* 131* Z (local output) DOUBLE PRECISION array, 132* global dimension (N, N), 133* local dimension ( LLD_Z, LOCc(JZ+N-1) ) 134* (see Notes below for more detailed explanation of 2d arrays) 135* If JOBZ = 'V', then on normal exit the first M columns of Z 136* contain the orthonormal eigenvectors of the matrix 137* corresponding to the selected eigenvalues. 138* If JOBZ = 'N', then Z is not referenced. 139* 140* IZ (global input) INTEGER 141* Z's global row index, which points to the beginning of the 142* submatrix which is to be operated on. 143* It should be set to 1 when operating on a full matrix. 144* 145* JZ (global input) INTEGER 146* Z's global column index, which points to the beginning of 147* the submatrix which is to be operated on. 148* It should be set to 1 when operating on a full matrix. 149* 150* DESCZ (global and local input) INTEGER array of dimension DLEN_. 151* The array descriptor for the distributed matrix Z. 152* The context DESCZ( CTXT_ ) must equal DESCA( CTXT_ ). 153* Also note the array alignment requirements specified below. 154* 155* WORK (local workspace/output) DOUBLE PRECISION array, 156* dimension (LWORK) 157* On return, WORK(1) contains the optimal amount of 158* workspace required for efficient execution. 159* if JOBZ='N' WORK(1) = optimal amount of workspace 160* required to compute the eigenvalues. 161* if JOBZ='V' WORK(1) = optimal amount of workspace 162* required to compute eigenvalues and eigenvectors. 163* 164* LWORK (local input) INTEGER 165* Size of WORK, must be at least 3. 166* See below for definitions of variables used to define LWORK. 167* If no eigenvectors are requested (JOBZ = 'N') then 168* LWORK >= 2 + 5*N + MAX( 12 * NN, NB * ( NP0 + 1 ) ) 169* If eigenvectors are requested (JOBZ = 'V' ) then 170* the amount of workspace required is: 171* LWORK >= 2 + 5*N + MAX( 18*NN, NP0 * MQ0 + 2 * NB * NB ) + 172* (2 + ICEIL( NEIG, NPROW*NPCOL))*NN 173* 174* Variable definitions: 175* NEIG = number of eigenvectors requested 176* NB = DESCA( MB_ ) = DESCA( NB_ ) = 177* DESCZ( MB_ ) = DESCZ( NB_ ) 178* NN = MAX( N, NB, 2 ) 179* DESCA( RSRC_ ) = DESCA( NB_ ) = DESCZ( RSRC_ ) = 180* DESCZ( CSRC_ ) = 0 181* NP0 = NUMROC( NN, NB, 0, 0, NPROW ) 182* MQ0 = NUMROC( MAX( NEIG, NB, 2 ), NB, 0, 0, NPCOL ) 183* ICEIL( X, Y ) is a ScaLAPACK function returning 184* ceiling(X/Y) 185* 186* If LWORK = -1, then LWORK is global input and a workspace 187* query is assumed; the routine only calculates the size 188* required for optimal performance for all work arrays. Each of 189* these values is returned in the first entry of the 190* corresponding work arrays, and no error message is issued by 191* PXERBLA. 192* Note that in a workspace query, for performance the optimal 193* workspace LWOPT is returned rather than the minimum necessary 194* WORKSPACE LWMIN. For very small matrices, LWOPT >> LWMIN. 195* 196* IWORK (local workspace) INTEGER array 197* On return, IWORK(1) contains the amount of integer workspace 198* required. 199* 200* LIWORK (local input) INTEGER 201* size of IWORK 202* 203* Let NNP = MAX( N, NPROW*NPCOL + 1, 4 ). Then: 204* LIWORK >= 12*NNP + 2*N when the eigenvectors are desired 205* LIWORK >= 10*NNP + 2*N when only the eigenvalues have to be computed 206* 207* If LIWORK = -1, then LIWORK is global input and a workspace 208* query is assumed; the routine only calculates the minimum 209* and optimal size for all work arrays. Each of these 210* values is returned in the first entry of the corresponding 211* work array, and no error message is issued by PXERBLA. 212* 213* INFO (global output) INTEGER 214* = 0: successful exit 215* < 0: If the i-th argument is an array and the j-entry had 216* an illegal value, then INFO = -(i*100+j), if the i-th 217* argument is a scalar and had an illegal value, then 218* INFO = -i. 219* 220* Notes 221* ===== 222* 223* Each global data object is described by an associated description 224* vector. This vector stores the information required to establish 225* the mapping between an object element and its corresponding process 226* and memory location. 227* 228* Let A be a generic term for any 2D block cyclicly distributed array. 229* Such a global array has an associated description vector DESCA, 230* or DESCZ for the descriptor of Z, etc. 231* The length of a ScaLAPACK descriptor is nine. 232* In the following comments, the character _ should be read as 233* "of the global array". 234* 235* NOTATION STORED IN EXPLANATION 236* --------------- -------------- -------------------------------------- 237* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 238* DTYPE_A = 1. 239* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 240* the BLACS process grid A is distribu- 241* ted over. The context itself is glo- 242* bal, but the handle (the integer 243* value) may vary. 244* M_A (global) DESCA( M_ ) The number of rows in the global 245* array A. 246* N_A (global) DESCA( N_ ) The number of columns in the global 247* array A. 248* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 249* the rows of the array. 250* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 251* the columns of the array. 252* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 253* row of the array A is distributed. 254* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 255* first column of the array A is 256* distributed. 257* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 258* array. LLD_A >= MAX(1,LOCr(M_A)). 259* 260* Let K be the number of rows or columns of a distributed matrix, 261* and assume that its process grid has dimension p x q. 262* LOCr( K ) denotes the number of elements of K that a process 263* would receive if K were distributed over the p processes of its 264* process column. 265* Similarly, LOCc( K ) denotes the number of elements of K that a 266* process would receive if K were distributed over the q processes of 267* its process row. 268* The values of LOCr() and LOCc() may be determined via a call to the 269* ScaLAPACK tool function, NUMROC: 270* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 271* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 272* An upper bound for these quantities may be computed by: 273* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 274* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 275* 276* PDSYEVR assumes IEEE 754 standard compliant arithmetic. 277* 278* Alignment requirements 279* ====================== 280* 281* The distributed submatrices A(IA:*, JA:*) and Z(IZ:IZ+M-1,JZ:JZ+N-1) 282* must satisfy the following alignment properties: 283* 284* 1.Identical (quadratic) dimension: 285* DESCA(M_) = DESCZ(M_) = DESCA(N_) = DESCZ(N_) 286* 2.Quadratic conformal blocking: 287* DESCA(MB_) = DESCA(NB_) = DESCZ(MB_) = DESCZ(NB_) 288* DESCA(RSRC_) = DESCZ(RSRC_) 289* 3.MOD( IA-1, MB_A ) = MOD( IZ-1, MB_Z ) = 0 290* 4.IAROW = IZROW 291* 292* 293* .. Parameters .. 294 INTEGER CTXT_, M_, N_, 295 $ MB_, NB_, RSRC_, CSRC_ 296 PARAMETER ( CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 297 $ RSRC_ = 7, CSRC_ = 8 ) 298 DOUBLE PRECISION ZERO 299 PARAMETER ( ZERO = 0.0D0 ) 300* .. 301* .. Local Scalars .. 302 LOGICAL ALLEIG, COLBRT, DOBCST, FINISH, FIRST, INDEIG, 303 $ LOWER, LQUERY, VALEIG, VSTART, WANTZ 304 INTEGER ANB, DOL, DOU, DSTCOL, DSTROW, EIGCNT, FRSTCL, 305 $ I, IAROW, ICTXT, IIL, IINDERR, IINDWLC, IINFO, 306 $ IIU, IM, INDD, INDD2, INDE, INDE2, INDERR, 307 $ INDILU, INDRW, INDTAU, INDWLC, INDWORK, IPIL, 308 $ IPIU, IPROC, IZROW, LASTCL, LENGTHI, LENGTHI2, 309 $ LIWMIN, LLWORK, LWMIN, LWOPT, MAXCLS, MQ00, 310 $ MYCOL, MYIL, MYIU, MYPROC, MYROW, MZ, NB, 311 $ NDEPTH, NEEDIL, NEEDIU, NNP, NP00, NPCOL, 312 $ NPROCS, NPROW, NPS, NSPLIT, NSYTRD_LWOPT, 313 $ OFFSET, PARITY, RLENGTHI, RLENGTHI2, RSTARTI, 314 $ SIZE1, SIZE2, SQNPC, SRCCOL, SRCROW, STARTI, 315 $ ZOFFSET 316 317 DOUBLE PRECISION PIVMIN, SAFMIN, SCALE, VLL, VUU, WL, 318 $ WU 319* 320* .. Local Arrays .. 321 INTEGER IDUM1( 4 ), IDUM2( 4 ) 322* .. 323* .. External Functions .. 324 LOGICAL LSAME 325 INTEGER ICEIL, INDXG2P, NUMROC, PJLAENV 326 DOUBLE PRECISION PDLAMCH 327 EXTERNAL ICEIL, INDXG2P, LSAME, NUMROC, PDLAMCH, 328 $ PJLAENV 329* .. 330* .. External Subroutines .. 331 EXTERNAL BLACS_GRIDINFO, CHK1MAT, DCOPY, DGEBR2D, 332 $ DGEBS2D, DGERV2D, DGESD2D, DLARRC, DLASRT2, 333 $ DSTEGR2A, DSTEGR2B, DSTEGR2, IGEBR2D, 334 $ IGEBS2D, IGERV2D, IGESD2D, IGSUM2D, PCHK1MAT, 335 $ PCHK2MAT, PDELGET, PDLAEVSWP, PDLARED1D, 336 $ PDORMTR, PDSYNTRD, PXERBLA 337* .. 338* .. Intrinsic Functions .. 339 INTRINSIC ABS, DBLE, ICHAR, INT, MAX, MIN, MOD, SQRT 340* .. 341* .. Executable Statements .. 342* 343 344 345 INFO = 0 346*********************************************************************** 347* 348* Decode character arguments to find out what the code should do 349* 350*********************************************************************** 351 WANTZ = LSAME( JOBZ, 'V' ) 352 LOWER = LSAME( UPLO, 'L' ) 353 ALLEIG = LSAME( RANGE, 'A' ) 354 VALEIG = LSAME( RANGE, 'V' ) 355 INDEIG = LSAME( RANGE, 'I' ) 356 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 357 358*********************************************************************** 359* 360* GET MACHINE PARAMETERS 361* 362*********************************************************************** 363 ICTXT = DESCA( CTXT_ ) 364 SAFMIN = PDLAMCH( ICTXT, 'Safe minimum' ) 365 366*********************************************************************** 367* 368* Set up pointers into the WORK array 369* 370*********************************************************************** 371 INDTAU = 1 372 INDD = INDTAU + N 373 INDE = INDD + N + 1 374 INDD2 = INDE + N + 1 375 INDE2 = INDD2 + N 376 INDWORK = INDE2 + N 377 LLWORK = LWORK - INDWORK + 1 378 379*********************************************************************** 380* 381* BLACS PROCESSOR GRID SETUP 382* 383*********************************************************************** 384 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 385 386 387 NPROCS = NPROW * NPCOL 388 MYPROC = MYROW * NPCOL + MYCOL 389 IF( NPROW.EQ.-1 ) THEN 390 INFO = -( 800+CTXT_ ) 391 ELSE IF( WANTZ ) THEN 392 IF( ICTXT.NE.DESCZ( CTXT_ ) ) THEN 393 INFO = -( 2100+CTXT_ ) 394 END IF 395 END IF 396 397*********************************************************************** 398* 399* COMPUTE REAL WORKSPACE 400* 401*********************************************************************** 402 IF ( ALLEIG ) THEN 403 MZ = N 404 ELSE IF ( INDEIG ) THEN 405 MZ = IU - IL + 1 406 ELSE 407* Take upper bound for VALEIG case 408 MZ = N 409 END IF 410* 411 NB = DESCA( NB_ ) 412 IF ( WANTZ ) THEN 413 NP00 = NUMROC( N, NB, 0, 0, NPROW ) 414 MQ00 = NUMROC( MZ, NB, 0, 0, NPCOL ) 415 INDRW = INDWORK + MAX(18*N, NP00*MQ00 + 2*NB*NB) 416 LWMIN = INDRW - 1 + (ICEIL(MZ, NPROCS) + 2)*N 417 ELSE 418 INDRW = INDWORK + 12*N 419 LWMIN = INDRW - 1 420 END IF 421* The code that validates the input requires 3 workspace entries 422 LWMIN = MAX(3, LWMIN) 423 LWOPT = LWMIN 424 ANB = PJLAENV( ICTXT, 3, 'PDSYTTRD', 'L', 0, 0, 0, 0 ) 425 SQNPC = INT( SQRT( DBLE( NPROCS ) ) ) 426 NPS = MAX( NUMROC( N, 1, 0, 0, SQNPC ), 2*ANB ) 427 NSYTRD_LWOPT = 2*( ANB+1 )*( 4*NPS+2 ) + ( NPS+4 )*NPS 428 LWOPT = MAX( LWOPT, 5*N+NSYTRD_LWOPT ) 429* 430 SIZE1 = INDRW - INDWORK 431 432*********************************************************************** 433* 434* COMPUTE INTEGER WORKSPACE 435* 436*********************************************************************** 437 NNP = MAX( N, NPROCS+1, 4 ) 438 IF ( WANTZ ) THEN 439 LIWMIN = 12*NNP + 2*N 440 ELSE 441 LIWMIN = 10*NNP + 2*N 442 END IF 443 444*********************************************************************** 445* 446* Set up pointers into the IWORK array 447* 448*********************************************************************** 449* Pointer to eigenpair distribution over processors 450 INDILU = LIWMIN - 2*NPROCS + 1 451 SIZE2 = INDILU - 2*N 452 453 454*********************************************************************** 455* 456* Test the input arguments. 457* 458*********************************************************************** 459 IF( INFO.EQ.0 ) THEN 460 CALL CHK1MAT( N, 4, N, 4, IA, JA, DESCA, 8, INFO ) 461 IF( WANTZ ) 462 $ CALL CHK1MAT( N, 4, N, 4, IZ, JZ, DESCZ, 21, INFO ) 463* 464 IF( INFO.EQ.0 ) THEN 465 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 466 INFO = -1 467 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 468 INFO = -2 469 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 470 INFO = -3 471 ELSE IF( MOD( IA-1, DESCA( MB_ ) ).NE.0 ) THEN 472 INFO = -6 473 ELSE IF( VALEIG .AND. N.GT.0 .AND. VU.LE.VL ) THEN 474 INFO = -10 475 ELSE IF( INDEIG .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) 476 $ THEN 477 INFO = -11 478 ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) 479 $ THEN 480 INFO = -12 481 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 482 INFO = -21 483 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 484 INFO = -23 485 ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN 486 INFO = -( 800+NB_ ) 487 END IF 488 IF( WANTZ ) THEN 489 IAROW = INDXG2P( 1, DESCA( NB_ ), MYROW, 490 $ DESCA( RSRC_ ), NPROW ) 491 IZROW = INDXG2P( 1, DESCA( NB_ ), MYROW, 492 $ DESCZ( RSRC_ ), NPROW ) 493 IF( IAROW.NE.IZROW ) THEN 494 INFO = -19 495 ELSE IF( MOD( IA-1, DESCA( MB_ ) ).NE. 496 $ MOD( IZ-1, DESCZ( MB_ ) ) ) THEN 497 INFO = -19 498 ELSE IF( DESCA( M_ ).NE.DESCZ( M_ ) ) THEN 499 INFO = -( 2100+M_ ) 500 ELSE IF( DESCA( N_ ).NE.DESCZ( N_ ) ) THEN 501 INFO = -( 2100+N_ ) 502 ELSE IF( DESCA( MB_ ).NE.DESCZ( MB_ ) ) THEN 503 INFO = -( 2100+MB_ ) 504 ELSE IF( DESCA( NB_ ).NE.DESCZ( NB_ ) ) THEN 505 INFO = -( 2100+NB_ ) 506 ELSE IF( DESCA( RSRC_ ).NE.DESCZ( RSRC_ ) ) THEN 507 INFO = -( 2100+RSRC_ ) 508 ELSE IF( DESCA( CSRC_ ).NE.DESCZ( CSRC_ ) ) THEN 509 INFO = -( 2100+CSRC_ ) 510 ELSE IF( ICTXT.NE.DESCZ( CTXT_ ) ) THEN 511 INFO = -( 2100+CTXT_ ) 512 END IF 513 END IF 514 END IF 515 IDUM2( 1 ) = 1 516 IF( LOWER ) THEN 517 IDUM1( 2 ) = ICHAR( 'L' ) 518 ELSE 519 IDUM1( 2 ) = ICHAR( 'U' ) 520 END IF 521 IDUM2( 2 ) = 2 522 IF( ALLEIG ) THEN 523 IDUM1( 3 ) = ICHAR( 'A' ) 524 ELSE IF( INDEIG ) THEN 525 IDUM1( 3 ) = ICHAR( 'I' ) 526 ELSE 527 IDUM1( 3 ) = ICHAR( 'V' ) 528 END IF 529 IDUM2( 3 ) = 3 530 IF( LQUERY ) THEN 531 IDUM1( 4 ) = -1 532 ELSE 533 IDUM1( 4 ) = 1 534 END IF 535 IDUM2( 4 ) = 4 536 IF( WANTZ ) THEN 537 IDUM1( 1 ) = ICHAR( 'V' ) 538 CALL PCHK2MAT( N, 4, N, 4, IA, JA, DESCA, 8, N, 4, N, 4, IZ, 539 $ JZ, DESCZ, 21, 4, IDUM1, IDUM2, INFO ) 540 ELSE 541 IDUM1( 1 ) = ICHAR( 'N' ) 542 CALL PCHK1MAT( N, 4, N, 4, IA, JA, DESCA, 8, 4, IDUM1, 543 $ IDUM2, INFO ) 544 END IF 545 WORK( 1 ) = DBLE( LWOPT ) 546 IWORK( 1 ) = LIWMIN 547 END IF 548* 549 IF( INFO.NE.0 ) THEN 550 CALL PXERBLA( ICTXT, 'PDSYEVR', -INFO ) 551 RETURN 552 ELSE IF( LQUERY ) THEN 553 RETURN 554 END IF 555 556*********************************************************************** 557* 558* Quick return if possible 559* 560*********************************************************************** 561 IF( N.EQ.0 ) THEN 562 IF( WANTZ ) THEN 563 NZ = 0 564 END IF 565 M = 0 566 WORK( 1 ) = DBLE( LWOPT ) 567 IWORK( 1 ) = LIWMIN 568 RETURN 569 END IF 570 571 IF( VALEIG ) THEN 572 VLL = VL 573 VUU = VU 574 ELSE 575 VLL = ZERO 576 VUU = ZERO 577 END IF 578* 579* No scaling done here, leave this to MRRR kernel. 580* Scale tridiagonal rather than full matrix. 581* 582*********************************************************************** 583* 584* REDUCE SYMMETRIC MATRIX TO TRIDIAGONAL FORM. 585* 586*********************************************************************** 587 588 589 CALL PDSYNTRD( UPLO, N, A, IA, JA, DESCA, WORK( INDD ), 590 $ WORK( INDE ), WORK( INDTAU ), WORK( INDWORK ), 591 $ LLWORK, IINFO ) 592 593 594 IF (IINFO .NE. 0) THEN 595 CALL PXERBLA( ICTXT, 'PDSYNTRD', -IINFO ) 596 RETURN 597 END IF 598 599*********************************************************************** 600* 601* DISTRIBUTE TRIDIAGONAL TO ALL PROCESSORS 602* 603*********************************************************************** 604 OFFSET = 0 605 IF( IA.EQ.1 .AND. JA.EQ.1 .AND. 606 $ DESCA( RSRC_ ).EQ.0 .AND. DESCA( CSRC_ ).EQ.0 ) 607 $ THEN 608 CALL PDLARED1D( N, IA, JA, DESCA, WORK( INDD ), WORK( INDD2 ), 609 $ WORK( INDWORK ), LLWORK ) 610* 611 CALL PDLARED1D( N, IA, JA, DESCA, WORK( INDE ), WORK( INDE2 ), 612 $ WORK( INDWORK ), LLWORK ) 613 IF( .NOT.LOWER ) 614 $ OFFSET = 1 615 ELSE 616 DO 10 I = 1, N 617 CALL PDELGET( 'A', ' ', WORK( INDD2+I-1 ), A, I+IA-1, 618 $ I+JA-1, DESCA ) 619 10 CONTINUE 620 IF( LSAME( UPLO, 'U' ) ) THEN 621 DO 20 I = 1, N - 1 622 CALL PDELGET( 'A', ' ', WORK( INDE2+I-1 ), A, I+IA-1, 623 $ I+JA, DESCA ) 624 20 CONTINUE 625 ELSE 626 DO 30 I = 1, N - 1 627 CALL PDELGET( 'A', ' ', WORK( INDE2+I-1 ), A, I+IA, 628 $ I+JA-1, DESCA ) 629 30 CONTINUE 630 END IF 631 END IF 632 633 634 635 636*********************************************************************** 637* 638* SET IIL, IIU 639* 640*********************************************************************** 641 IF ( ALLEIG ) THEN 642 IIL = 1 643 IIU = N 644 ELSE IF ( INDEIG ) THEN 645 IIL = IL 646 IIU = IU 647 ELSE IF ( VALEIG ) THEN 648 CALL DLARRC('T', N, VLL, VUU, WORK( INDD2 ), 649 $ WORK( INDE2 + OFFSET ), SAFMIN, EIGCNT, IIL, IIU, INFO) 650* Refine upper bound N that was taken 651 MZ = EIGCNT 652 IIL = IIL + 1 653 ENDIF 654 655 IF(MZ.EQ.0) THEN 656 M = 0 657 IF( WANTZ ) THEN 658 NZ = 0 659 END IF 660 WORK( 1 ) = DBLE( LWOPT ) 661 IWORK( 1 ) = LIWMIN 662 RETURN 663 END IF 664 665 MYIL = 0 666 MYIU = 0 667 M = 0 668 IM = 0 669 670*********************************************************************** 671* 672* COMPUTE WORK ASSIGNMENTS 673* 674*********************************************************************** 675* 676* Each processor computes the work assignments for all processors 677* 678 CALL PMPIM2( IIL, IIU, NPROCS, 679 $ IWORK(INDILU), IWORK(INDILU+NPROCS) ) 680* 681* Find local work assignment 682* 683 MYIL = IWORK(INDILU+MYPROC) 684 MYIU = IWORK(INDILU+NPROCS+MYPROC) 685 686 687 ZOFFSET = MAX(0, MYIL - IIL - 1) 688 FIRST = ( MYIL .EQ. IIL ) 689 690 691*********************************************************************** 692* 693* CALLS TO MRRR KERNEL 694* 695*********************************************************************** 696 IF(.NOT.WANTZ) THEN 697* 698* Compute eigenvalues only. 699* 700 IINFO = 0 701 IF ( MYIL.GT.0 ) THEN 702 DOL = 1 703 DOU = MYIU - MYIL + 1 704 CALL DSTEGR2( JOBZ, 'I', N, WORK( INDD2 ), 705 $ WORK( INDE2+OFFSET ), VLL, VUU, MYIL, MYIU, 706 $ IM, W( 1 ), WORK( INDRW ), N, 707 $ MYIU - MYIL + 1, 708 $ IWORK( 1 ), WORK( INDWORK ), SIZE1, 709 $ IWORK( 2*N+1 ), SIZE2, 710 $ DOL, DOU, ZOFFSET, IINFO ) 711* DSTEGR2 zeroes out the entire W array, so we can't just give 712* it the part of W we need. So here we copy the W entries into 713* their correct location 714 DO 49 I = 1, IM 715 W( MYIL-IIL+I ) = W( I ) 716 49 CONTINUE 717* W( MYIL ) is at W( MYIL - IIL + 1 ) 718* W( X ) is at W(X - IIL + 1 ) 719 END IF 720 IF (IINFO .NE. 0) THEN 721 CALL PXERBLA( ICTXT, 'DSTEGR2', -IINFO ) 722 RETURN 723 END IF 724 ELSEIF ( WANTZ .AND. NPROCS.EQ.1 ) THEN 725* 726* Compute eigenvalues and -vectors, but only on one processor 727* 728 IINFO = 0 729 IF ( MYIL.GT.0 ) THEN 730 DOL = MYIL - IIL + 1 731 DOU = MYIU - IIL + 1 732 CALL DSTEGR2( JOBZ, 'I', N, WORK( INDD2 ), 733 $ WORK( INDE2+OFFSET ), VLL, VUU, IIL, IIU, 734 $ IM, W( 1 ), WORK( INDRW ), N, 735 $ N, 736 $ IWORK( 1 ), WORK( INDWORK ), SIZE1, 737 $ IWORK( 2*N+1 ), SIZE2, DOL, DOU, 738 $ ZOFFSET, IINFO ) 739 ENDIF 740 IF (IINFO .NE. 0) THEN 741 CALL PXERBLA( ICTXT, 'DSTEGR2', -IINFO ) 742 RETURN 743 END IF 744 ELSEIF ( WANTZ ) THEN 745* 746* Compute representations in parallel. 747* Share eigenvalue computation for root between all processors 748* Then compute the eigenvectors. 749* 750 IINFO = 0 751* Part 1. compute root representations and root eigenvalues 752 IF ( MYIL.GT.0 ) THEN 753 DOL = MYIL - IIL + 1 754 DOU = MYIU - IIL + 1 755 CALL DSTEGR2A( JOBZ, 'I', N, WORK( INDD2 ), 756 $ WORK( INDE2+OFFSET ), VLL, VUU, IIL, IIU, 757 $ IM, W( 1 ), WORK( INDRW ), N, 758 $ N, WORK( INDWORK ), SIZE1, 759 $ IWORK( 2*N+1 ), SIZE2, DOL, 760 $ DOU, NEEDIL, NEEDIU, 761 $ INDERR, NSPLIT, PIVMIN, SCALE, WL, WU, 762 $ IINFO ) 763 ENDIF 764 IF (IINFO .NE. 0) THEN 765 CALL PXERBLA( ICTXT, 'DSTEGR2A', -IINFO ) 766 RETURN 767 END IF 768* 769* The second part of parallel MRRR, the representation tree 770* construction begins. Upon successful completion, the 771* eigenvectors have been computed. This is indicated by 772* the flag FINISH. 773* 774 VSTART = .TRUE. 775 FINISH = (MYIL.LE.0) 776C Part 2. Share eigenvalues and uncertainties between all processors 777 IINDERR = INDWORK + INDERR - 1 778 779* 780* 781* There are currently two ways to communicate eigenvalue information 782* using the BLACS. 783* 1.) BROADCAST 784* 2.) POINT2POINT between collaborators (those processors working 785* jointly on a cluster. 786* For efficiency, BROADCAST has been disabled. 787* At a later stage, other more efficient communication algorithms 788* might be implemented, e. g. group or tree-based communication. 789* 790 DOBCST = .FALSE. 791 IF(DOBCST) THEN 792* First gather everything on the first processor. 793* Then use BROADCAST-based communication 794 DO 45 I = 2, NPROCS 795 IF (MYPROC .EQ. (I - 1)) THEN 796 DSTROW = 0 797 DSTCOL = 0 798 STARTI = DOL 799 IWORK(1) = STARTI 800 IF(MYIL.GT.0) THEN 801 LENGTHI = MYIU - MYIL + 1 802 ELSE 803 LENGTHI = 0 804 ENDIF 805 IWORK(2) = LENGTHI 806 CALL IGESD2D( ICTXT, 2, 1, IWORK, 2, 807 $ DSTROW, DSTCOL ) 808 IF (( STARTI.GE.1 ) .AND. ( LENGTHI.GE.1 )) THEN 809 LENGTHI2 = 2*LENGTHI 810* Copy eigenvalues into communication buffer 811 CALL DCOPY(LENGTHI,W( STARTI ),1, 812 $ WORK( INDD ), 1) 813* Copy uncertainties into communication buffer 814 CALL DCOPY(LENGTHI,WORK( IINDERR+STARTI-1 ),1, 815 $ WORK( INDD+LENGTHI ), 1) 816* send buffer 817 CALL DGESD2D( ICTXT, LENGTHI2, 818 $ 1, WORK( INDD ), LENGTHI2, 819 $ DSTROW, DSTCOL ) 820 END IF 821 ELSE IF (MYPROC .EQ. 0) THEN 822 SRCROW = (I-1) / NPCOL 823 SRCCOL = MOD(I-1, NPCOL) 824 CALL IGERV2D( ICTXT, 2, 1, IWORK, 2, 825 $ SRCROW, SRCCOL ) 826 STARTI = IWORK(1) 827 LENGTHI = IWORK(2) 828 IF (( STARTI.GE.1 ) .AND. ( LENGTHI.GE.1 )) THEN 829 LENGTHI2 = 2*LENGTHI 830* receive buffer 831 CALL DGERV2D( ICTXT, LENGTHI2, 1, 832 $ WORK(INDD), LENGTHI2, SRCROW, SRCCOL ) 833* copy eigenvalues from communication buffer 834 CALL DCOPY( LENGTHI, WORK(INDD), 1, 835 $ W( STARTI ), 1) 836* copy uncertainties (errors) from communication buffer 837 CALL DCOPY(LENGTHI,WORK(INDD+LENGTHI),1, 838 $ WORK( IINDERR+STARTI-1 ), 1) 839 END IF 840 END IF 841 45 CONTINUE 842 LENGTHI = IIU - IIL + 1 843 LENGTHI2 = LENGTHI * 2 844 IF (MYPROC .EQ. 0) THEN 845* Broadcast eigenvalues and errors to all processors 846 CALL DCOPY(LENGTHI,W ,1, WORK( INDD ), 1) 847 CALL DCOPY(LENGTHI,WORK( IINDERR ),1, 848 $ WORK( INDD+LENGTHI ), 1) 849 CALL DGEBS2D( ICTXT, 'A', ' ', LENGTHI2, 1, 850 $ WORK(INDD), LENGTHI2 ) 851 ELSE 852 SRCROW = 0 853 SRCCOL = 0 854 CALL DGEBR2D( ICTXT, 'A', ' ', LENGTHI2, 1, 855 $ WORK(INDD), LENGTHI2, SRCROW, SRCCOL ) 856 CALL DCOPY( LENGTHI, WORK(INDD), 1, W, 1) 857 CALL DCOPY(LENGTHI,WORK(INDD+LENGTHI),1, 858 $ WORK( IINDERR ), 1) 859 END IF 860 ELSE 861* 862* Enable point2point communication between collaborators 863* 864* Find collaborators of MYPROC 865 IF( (NPROCS.GT.1).AND.(MYIL.GT.0) ) THEN 866 CALL PMPCOL( MYPROC, NPROCS, IIL, NEEDIL, NEEDIU, 867 $ IWORK(INDILU), IWORK(INDILU+NPROCS), 868 $ COLBRT, FRSTCL, LASTCL ) 869 ELSE 870 COLBRT = .FALSE. 871 ENDIF 872 873 IF(COLBRT) THEN 874* If the processor collaborates with others, 875* communicate information. 876 DO 47 IPROC = FRSTCL, LASTCL 877 IF (MYPROC .EQ. IPROC) THEN 878 STARTI = DOL 879 IWORK(1) = STARTI 880 LENGTHI = MYIU - MYIL + 1 881 IWORK(2) = LENGTHI 882 883 IF ((STARTI.GE.1) .AND. (LENGTHI.GE.1)) THEN 884* Copy eigenvalues into communication buffer 885 CALL DCOPY(LENGTHI,W( STARTI ),1, 886 $ WORK(INDD), 1) 887* Copy uncertainties into communication buffer 888 CALL DCOPY(LENGTHI, 889 $ WORK( IINDERR+STARTI-1 ),1, 890 $ WORK(INDD+LENGTHI), 1) 891 ENDIF 892 893 DO 46 I = FRSTCL, LASTCL 894 IF(I.EQ.MYPROC) GOTO 46 895 DSTROW = I/ NPCOL 896 DSTCOL = MOD(I, NPCOL) 897 CALL IGESD2D( ICTXT, 2, 1, IWORK, 2, 898 $ DSTROW, DSTCOL ) 899 IF ((STARTI.GE.1) .AND. (LENGTHI.GE.1)) THEN 900 LENGTHI2 = 2*LENGTHI 901* send buffer 902 CALL DGESD2D( ICTXT, LENGTHI2, 903 $ 1, WORK(INDD), LENGTHI2, 904 $ DSTROW, DSTCOL ) 905 END IF 906 46 CONTINUE 907 ELSE 908 SRCROW = IPROC / NPCOL 909 SRCCOL = MOD(IPROC, NPCOL) 910 CALL IGERV2D( ICTXT, 2, 1, IWORK, 2, 911 $ SRCROW, SRCCOL ) 912 RSTARTI = IWORK(1) 913 RLENGTHI = IWORK(2) 914 IF ((RSTARTI.GE.1 ) .AND. (RLENGTHI.GE.1 )) THEN 915 RLENGTHI2 = 2*RLENGTHI 916 CALL DGERV2D( ICTXT, RLENGTHI2, 1, 917 $ WORK(INDE), RLENGTHI2, 918 $ SRCROW, SRCCOL ) 919* copy eigenvalues from communication buffer 920 CALL DCOPY( RLENGTHI, WORK(INDE), 1, 921 $ W( RSTARTI ), 1) 922* copy uncertainties (errors) from communication buffer 923 CALL DCOPY(RLENGTHI,WORK(INDE+RLENGTHI),1, 924 $ WORK( IINDERR+RSTARTI-1 ), 1) 925 END IF 926 END IF 927 47 CONTINUE 928 ENDIF 929 ENDIF 930 931* 932* Part 3. Compute representation tree and eigenvectors. 933* What follows is a loop in which the tree 934* is constructed in parallel from top to bottom, 935* on level at a time, until all eigenvectors 936* have been computed. 937* 938 100 CONTINUE 939 IF ( MYIL.GT.0 ) THEN 940 CALL DSTEGR2B( JOBZ, N, WORK( INDD2 ), 941 $ WORK( INDE2+OFFSET ), 942 $ IM, W( 1 ), WORK( INDRW ), N, N, 943 $ IWORK( 1 ), WORK( INDWORK ), SIZE1, 944 $ IWORK( 2*N+1 ), SIZE2, DOL, 945 $ DOU, NEEDIL, NEEDIU, INDWLC, 946 $ PIVMIN, SCALE, WL, WU, 947 $ VSTART, FINISH, 948 $ MAXCLS, NDEPTH, PARITY, ZOFFSET, IINFO ) 949 IINDWLC = INDWORK + INDWLC - 1 950 IF(.NOT.FINISH) THEN 951 IF((NEEDIL.LT.DOL).OR.(NEEDIU.GT.DOU)) THEN 952 CALL PMPCOL( MYPROC, NPROCS, IIL, NEEDIL, NEEDIU, 953 $ IWORK(INDILU), IWORK(INDILU+NPROCS), 954 $ COLBRT, FRSTCL, LASTCL ) 955 ELSE 956 COLBRT = .FALSE. 957 FRSTCL = MYPROC 958 LASTCL = MYPROC 959 ENDIF 960* 961* Check if this processor collaborates, i.e. 962* communication is needed. 963* 964 IF(COLBRT) THEN 965 DO 147 IPROC = FRSTCL, LASTCL 966 IF (MYPROC .EQ. IPROC) THEN 967 STARTI = DOL 968 IWORK(1) = STARTI 969 IF(MYIL.GT.0) THEN 970 LENGTHI = MYIU - MYIL + 1 971 ELSE 972 LENGTHI = 0 973 ENDIF 974 IWORK(2) = LENGTHI 975 IF ((STARTI.GE.1).AND.(LENGTHI.GE.1)) THEN 976* Copy eigenvalues into communication buffer 977 CALL DCOPY(LENGTHI, 978 $ WORK( IINDWLC+STARTI-1 ),1, 979 $ WORK(INDD), 1) 980* Copy uncertainties into communication buffer 981 CALL DCOPY(LENGTHI, 982 $ WORK( IINDERR+STARTI-1 ),1, 983 $ WORK(INDD+LENGTHI), 1) 984 ENDIF 985 986 DO 146 I = FRSTCL, LASTCL 987 IF(I.EQ.MYPROC) GOTO 146 988 DSTROW = I/ NPCOL 989 DSTCOL = MOD(I, NPCOL) 990 CALL IGESD2D( ICTXT, 2, 1, IWORK, 2, 991 $ DSTROW, DSTCOL ) 992 IF ((STARTI.GE.1).AND.(LENGTHI.GE.1)) THEN 993 LENGTHI2 = 2*LENGTHI 994* send buffer 995 CALL DGESD2D( ICTXT, LENGTHI2, 996 $ 1, WORK(INDD), LENGTHI2, 997 $ DSTROW, DSTCOL ) 998 END IF 999 146 CONTINUE 1000 ELSE 1001 SRCROW = IPROC / NPCOL 1002 SRCCOL = MOD(IPROC, NPCOL) 1003 CALL IGERV2D( ICTXT, 2, 1, IWORK, 2, 1004 $ SRCROW, SRCCOL ) 1005 RSTARTI = IWORK(1) 1006 RLENGTHI = IWORK(2) 1007 IF ((RSTARTI.GE.1).AND.(RLENGTHI.GE.1)) THEN 1008 RLENGTHI2 = 2*RLENGTHI 1009 CALL DGERV2D( ICTXT,RLENGTHI2, 1, 1010 $ WORK(INDE),RLENGTHI2, 1011 $ SRCROW, SRCCOL ) 1012* copy eigenvalues from communication buffer 1013 CALL DCOPY(RLENGTHI, WORK(INDE), 1, 1014 $ WORK( IINDWLC+RSTARTI-1 ), 1) 1015* copy uncertainties (errors) from communication buffer 1016 CALL DCOPY(RLENGTHI,WORK(INDE+RLENGTHI),1, 1017 $ WORK( IINDERR+RSTARTI-1 ), 1) 1018 END IF 1019 END IF 1020 147 CONTINUE 1021 ENDIF 1022 GOTO 100 1023 ENDIF 1024 ENDIF 1025 IF (IINFO .NE. 0) THEN 1026 CALL PXERBLA( ICTXT, 'DSTEGR2B', -IINFO ) 1027 RETURN 1028 END IF 1029* 1030 ENDIF 1031 1032* 1033*********************************************************************** 1034* 1035* MAIN PART ENDS HERE 1036* 1037*********************************************************************** 1038* 1039*********************************************************************** 1040* 1041* ALLGATHER: EACH PROCESSOR SENDS ITS EIGENVALUES TO THE FIRST ONE, 1042* THEN THE FIRST PROCESSOR BROADCASTS ALL EIGENVALUES 1043* 1044*********************************************************************** 1045* 1046 DO 50 I = 2, NPROCS 1047 IF (MYPROC .EQ. (I - 1)) THEN 1048 DSTROW = 0 1049 DSTCOL = 0 1050 STARTI = MYIL - IIL + 1 1051 IWORK(1) = STARTI 1052 IF(MYIL.GT.0) THEN 1053 LENGTHI = MYIU - MYIL + 1 1054 ELSE 1055 LENGTHI = 0 1056 ENDIF 1057 IWORK(2) = LENGTHI 1058 CALL IGESD2D( ICTXT, 2, 1, IWORK, 2, 1059 $ DSTROW, DSTCOL ) 1060 IF ((STARTI.GE.1).AND.(LENGTHI.GE.1)) THEN 1061 CALL DGESD2D( ICTXT, LENGTHI, 1062 $ 1, W( STARTI ), LENGTHI, 1063 $ DSTROW, DSTCOL ) 1064 ENDIF 1065 ELSE IF (MYPROC .EQ. 0) THEN 1066 SRCROW = (I-1) / NPCOL 1067 SRCCOL = MOD(I-1, NPCOL) 1068 CALL IGERV2D( ICTXT, 2, 1, IWORK, 2, 1069 $ SRCROW, SRCCOL ) 1070 STARTI = IWORK(1) 1071 LENGTHI = IWORK(2) 1072 IF ((STARTI.GE.1).AND.(LENGTHI.GE.1)) THEN 1073 CALL DGERV2D( ICTXT, LENGTHI, 1, 1074 $ W( STARTI ), LENGTHI, SRCROW, SRCCOL ) 1075 ENDIF 1076 ENDIF 1077 50 CONTINUE 1078 1079* Accumulate M from all processors 1080 M = IM 1081 CALL IGSUM2D( ICTXT, 'A', ' ', 1, 1, M, 1, -1, -1 ) 1082 1083* Broadcast eigenvalues to all processors 1084 IF (MYPROC .EQ. 0) THEN 1085* Send eigenvalues 1086 CALL DGEBS2D( ICTXT, 'A', ' ', M, 1, W, M ) 1087 ELSE 1088 SRCROW = 0 1089 SRCCOL = 0 1090 CALL DGEBR2D( ICTXT, 'A', ' ', M, 1, 1091 $ W, M, SRCROW, SRCCOL ) 1092 END IF 1093* 1094* Sort the eigenvalues and keep permutation in IWORK to 1095* sort the eigenvectors accordingly 1096* 1097 DO 160 I = 1, M 1098 IWORK( NPROCS+1+I ) = I 1099 160 CONTINUE 1100 CALL DLASRT2( 'I', M, W, IWORK( NPROCS+2 ), IINFO ) 1101 IF (IINFO.NE.0) THEN 1102 CALL PXERBLA( ICTXT, 'DLASRT2', -IINFO ) 1103 RETURN 1104 END IF 1105 1106*********************************************************************** 1107* 1108* TRANSFORM Z FROM 1D WORKSPACE INTO 2D BLOCKCYCLIC STORAGE 1109* 1110*********************************************************************** 1111 IF ( WANTZ ) THEN 1112 DO 170 I = 1, M 1113 IWORK( M+NPROCS+1+IWORK( NPROCS+1+I ) ) = I 1114 170 CONTINUE 1115* Store NVS in IWORK(1:NPROCS+1) for PDLAEVSWP 1116 IWORK( 1 ) = 0 1117 DO 180 I = 1, NPROCS 1118* Find IL and IU for processor i-1 1119* Has already been computed by PMPIM2 and stored 1120 IPIL = IWORK(INDILU+I-1) 1121 IPIU = IWORK(INDILU+NPROCS+I-1) 1122 IF (IPIL .EQ. 0) THEN 1123 IWORK( I + 1 ) = IWORK( I ) 1124 ELSE 1125 IWORK( I + 1 ) = IWORK( I ) + IPIU - IPIL + 1 1126 ENDIF 1127 180 CONTINUE 1128 1129 IF ( FIRST ) THEN 1130 CALL PDLAEVSWP(N, WORK( INDRW ), N, Z, IZ, JZ, 1131 $ DESCZ, IWORK( 1 ), IWORK( NPROCS+M+2 ), WORK( INDWORK ), 1132 $ INDRW - INDWORK ) 1133 ELSE 1134 CALL PDLAEVSWP(N, WORK( INDRW + N ), N, Z, IZ, JZ, 1135 $ DESCZ, IWORK( 1 ), IWORK( NPROCS+M+2 ), WORK( INDWORK ), 1136 $ INDRW - INDWORK ) 1137 END IF 1138* 1139 NZ = M 1140* 1141 1142*********************************************************************** 1143* 1144* Compute eigenvectors of A from eigenvectors of T 1145* 1146*********************************************************************** 1147 IF( NZ.GT.0 ) THEN 1148 CALL PDORMTR( 'L', UPLO, 'N', N, NZ, A, IA, JA, DESCA, 1149 $ WORK( INDTAU ), Z, IZ, JZ, DESCZ, 1150 $ WORK( INDWORK ), SIZE1, IINFO ) 1151 END IF 1152 IF (IINFO.NE.0) THEN 1153 CALL PXERBLA( ICTXT, 'PDORMTR', -IINFO ) 1154 RETURN 1155 END IF 1156* 1157 1158 END IF 1159* 1160 WORK( 1 ) = DBLE( LWOPT ) 1161 IWORK( 1 ) = LIWMIN 1162 1163 RETURN 1164* 1165* End of PDSYEVR 1166* 1167 END 1168