1 SUBROUTINE PDHSEQR( JOB, COMPZ, N, ILO, IHI, H, DESCH, WR, WI, Z, 2 $ DESCZ, WORK, LWORK, IWORK, LIWORK, INFO ) 3* 4* Contribution from the Department of Computing Science and HPC2N, 5* Umea University, Sweden 6* 7* -- ScaLAPACK driver routine (version 2.0.1) -- 8* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 9* Univ. of Colorado Denver and University of California, Berkeley. 10* January, 2012 11* 12 IMPLICIT NONE 13* 14* .. Scalar Arguments .. 15 INTEGER IHI, ILO, INFO, LWORK, LIWORK, N 16 CHARACTER COMPZ, JOB 17* .. 18* .. Array Arguments .. 19 INTEGER DESCH( * ) , DESCZ( * ), IWORK( * ) 20 DOUBLE PRECISION H( * ), WI( N ), WORK( * ), WR( N ), Z( * ) 21* .. 22* Purpose 23* ======= 24* 25* PDHSEQR computes the eigenvalues of an upper Hessenberg matrix H 26* and, optionally, the matrices T and Z from the Schur decomposition 27* H = Z*T*Z**T, where T is an upper quasi-triangular matrix (the 28* Schur form), and Z is the orthogonal matrix of Schur vectors. 29* 30* Optionally Z may be postmultiplied into an input orthogonal 31* matrix Q so that this routine can give the Schur factorization 32* of a matrix A which has been reduced to the Hessenberg form H 33* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. 34* 35* Notes 36* ===== 37* 38* Each global data object is described by an associated description 39* vector. This vector stores the information required to establish 40* the mapping between an object element and its corresponding process 41* and memory location. 42* 43* Let A be a generic term for any 2D block cyclicly distributed array. 44* Such a global array has an associated description vector DESCA. 45* In the following comments, the character _ should be read as 46* "of the global array". 47* 48* NOTATION STORED IN EXPLANATION 49* --------------- -------------- -------------------------------------- 50* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 51* DTYPE_A = 1. 52* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 53* the BLACS process grid A is distribu- 54* ted over. The context itself is glo- 55* bal, but the handle (the integer 56* value) may vary. 57* M_A (global) DESCA( M_ ) The number of rows in the global 58* array A. 59* N_A (global) DESCA( N_ ) The number of columns in the global 60* array A. 61* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 62* the rows of the array. 63* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 64* the columns of the array. 65* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 66* row of the array A is distributed. 67* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 68* first column of the array A is 69* distributed. 70* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 71* array. LLD_A >= MAX(1,LOCr(M_A)). 72* 73* Let K be the number of rows or columns of a distributed matrix, 74* and assume that its process grid has dimension p x q. 75* LOCr( K ) denotes the number of elements of K that a process 76* would receive if K were distributed over the p processes of its 77* process column. 78* Similarly, LOCc( K ) denotes the number of elements of K that a 79* process would receive if K were distributed over the q processes of 80* its process row. 81* The values of LOCr() and LOCc() may be determined via a call to the 82* ScaLAPACK tool function, NUMROC: 83* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 84* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 85* An upper bound for these quantities may be computed by: 86* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 87* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 88* 89* Arguments 90* ========= 91* 92* JOB (global input) CHARACTER*1 93* = 'E': compute eigenvalues only; 94* = 'S': compute eigenvalues and the Schur form T. 95* 96* COMPZ (global input) CHARACTER*1 97* = 'N': no Schur vectors are computed; 98* = 'I': Z is initialized to the unit matrix and the matrix Z 99* of Schur vectors of H is returned; 100* = 'V': Z must contain an orthogonal matrix Q on entry, and 101* the product Q*Z is returned. 102* 103* N (global input) INTEGER 104* The order of the Hessenberg matrix H (and Z if WANTZ). 105* N >= 0. 106* 107* ILO (global input) INTEGER 108* IHI (global input) INTEGER 109* It is assumed that H is already upper triangular in rows 110* and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally 111* set by a previous call to PDGEBAL, and then passed to PDGEHRD 112* when the matrix output by PDGEBAL is reduced to Hessenberg 113* form. Otherwise ILO and IHI should be set to 1 and N 114* respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. 115* If N = 0, then ILO = 1 and IHI = 0. 116* 117* H (global input/output) DOUBLE PRECISION array, dimension 118* (DESCH(LLD_),*) 119* On entry, the upper Hessenberg matrix H. 120* On exit, if JOB = 'S', H is upper quasi-triangular in 121* rows and columns ILO:IHI, with 1-by-1 and 2-by-2 blocks on 122* the main diagonal. The 2-by-2 diagonal blocks (corresponding 123* to complex conjugate pairs of eigenvalues) are returned in 124* standard form, with H(i,i) = H(i+1,i+1) and 125* H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the 126* contents of H are unspecified on exit. 127* 128* DESCH (global and local input) INTEGER array of dimension DLEN_. 129* The array descriptor for the distributed matrix H. 130* 131* WR (global output) DOUBLE PRECISION array, dimension (N) 132* WI (global output) DOUBLE PRECISION array, dimension (N) 133* The real and imaginary parts, respectively, of the computed 134* eigenvalues ILO to IHI are stored in the corresponding 135* elements of WR and WI. If two eigenvalues are computed as a 136* complex conjugate pair, they are stored in consecutive 137* elements of WR and WI, say the i-th and (i+1)th, with 138* WI(i) > 0 and WI(i+1) < 0. If JOB = 'S', the 139* eigenvalues are stored in the same order as on the diagonal 140* of the Schur form returned in H. 141* 142* Z (global input/output) DOUBLE PRECISION array. 143* If COMPZ = 'V', on entry Z must contain the current 144* matrix Z of accumulated transformations from, e.g., PDGEHRD, 145* and on exit Z has been updated; transformations are applied 146* only to the submatrix Z(ILO:IHI,ILO:IHI). 147* If COMPZ = 'N', Z is not referenced. 148* If COMPZ = 'I', on entry Z need not be set and on exit, 149* if INFO = 0, Z contains the orthogonal matrix Z of the Schur 150* vectors of H. 151* 152* DESCZ (global and local input) INTEGER array of dimension DLEN_. 153* The array descriptor for the distributed matrix Z. 154* 155* WORK (local workspace) DOUBLE PRECISION array, dimension(LWORK) 156* 157* LWORK (local input) INTEGER 158* The length of the workspace array WORK. 159* 160* IWORK (local workspace) INTEGER array, dimension (LIWORK) 161* 162* LIWORK (local input) INTEGER 163* The length of the workspace array IWORK. 164* 165* INFO (output) INTEGER 166* = 0: successful exit 167* .LT. 0: if INFO = -i, the i-th argument had an illegal 168* value (see also below for -7777 and -8888). 169* .GT. 0: if INFO = i, PDHSEQR failed to compute all of 170* the eigenvalues. Elements 1:ilo-1 and i+1:n of WR 171* and WI contain those eigenvalues which have been 172* successfully computed. (Failures are rare.) 173* 174* If INFO .GT. 0 and JOB = 'E', then on exit, the 175* remaining unconverged eigenvalues are the eigen- 176* values of the upper Hessenberg matrix rows and 177* columns ILO through INFO of the final, output 178* value of H. 179* 180* If INFO .GT. 0 and JOB = 'S', then on exit 181* 182* (*) (initial value of H)*U = U*(final value of H) 183* 184* where U is an orthogonal matrix. The final 185* value of H is upper Hessenberg and quasi-triangular 186* in rows and columns INFO+1 through IHI. 187* 188* If INFO .GT. 0 and COMPZ = 'V', then on exit 189* 190* (final value of Z) = (initial value of Z)*U 191* 192* where U is the orthogonal matrix in (*) (regard- 193* less of the value of JOB.) 194* 195* If INFO .GT. 0 and COMPZ = 'I', then on exit 196* (final value of Z) = U 197* where U is the orthogonal matrix in (*) (regard- 198* less of the value of JOB.) 199* 200* If INFO .GT. 0 and COMPZ = 'N', then Z is not 201* accessed. 202* 203* = -7777: PDLAQR0 failed to converge and PDLAQR1 was called 204* instead. This could happen. Mostly due to a bug. 205* Please, send a bug report to the authors. 206* = -8888: PDLAQR1 failed to converge and PDLAQR0 was called 207* instead. This should not happen. 208* 209* ================================================================ 210* Based on contributions by 211* Robert Granat, Department of Computing Science and HPC2N, 212* Umea University, Sweden. 213* ================================================================ 214* 215* Restrictions: The block size in H and Z must be square and larger 216* than or equal to six (6) due to restrictions in PDLAQR1, PDLAQR5 217* and DLAQR6. Moreover, H and Z need to be distributed identically 218* with the same context. 219* 220* ================================================================ 221* References: 222* K. Braman, R. Byers, and R. Mathias, 223* The Multi-Shift QR Algorithm Part I: Maintaining Well Focused 224* Shifts, and Level 3 Performance. 225* SIAM J. Matrix Anal. Appl., 23(4):929--947, 2002. 226* 227* K. Braman, R. Byers, and R. Mathias, 228* The Multi-Shift QR Algorithm Part II: Aggressive Early 229* Deflation. 230* SIAM J. Matrix Anal. Appl., 23(4):948--973, 2002. 231* 232* R. Granat, B. Kagstrom, and D. Kressner, 233* A Novel Parallel QR Algorithm for Hybrid Distributed Momory HPC 234* Systems. 235* SIAM J. Sci. Comput., 32(4):2345--2378, 2010. 236* 237* ================================================================ 238* .. Parameters .. 239 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 240 $ LLD_, MB_, M_, NB_, N_, RSRC_ 241 LOGICAL CRSOVER 242 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 243 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 244 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9, 245 $ CRSOVER = .TRUE. ) 246 INTEGER NTINY 247 PARAMETER ( NTINY = 11 ) 248 INTEGER NL 249 PARAMETER ( NL = 49 ) 250 DOUBLE PRECISION ZERO, ONE 251 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 252* .. 253* .. Local Scalars .. 254 INTEGER I, KBOT, NMIN, LLDH, LLDZ, ICTXT, NPROW, NPCOL, 255 $ MYROW, MYCOL, HROWS, HCOLS, IPW, NH, NB, 256 $ II, JJ, HRSRC, HCSRC, NPROCS, ILOC1, JLOC1, 257 $ HRSRC1, HCSRC1, K, ILOC2, JLOC2, ILOC3, JLOC3, 258 $ ILOC4, JLOC4, HRSRC2, HCSRC2, HRSRC3, HCSRC3, 259 $ HRSRC4, HCSRC4, LIWKOPT 260 LOGICAL INITZ, LQUERY, WANTT, WANTZ, PAIR, BORDER 261 DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4, DUM1, DUM2, DUM3, 262 $ DUM4, ELEM1, ELEM2, ELEM3, ELEM4, 263 $ CS, SN, ELEM5, TMP, LWKOPT 264* .. 265* .. Local Arrays .. 266 INTEGER DESCH2( DLEN_ ) 267* .. 268* .. External Functions .. 269 INTEGER PILAENVX, NUMROC, ICEIL 270 LOGICAL LSAME 271 EXTERNAL PILAENVX, LSAME, NUMROC, ICEIL 272* .. 273* .. External Subroutines .. 274 EXTERNAL PDLACPY, PDLAQR1, PDLAQR0, PDLASET, PXERBLA 275* .. 276* .. Intrinsic Functions .. 277 INTRINSIC DBLE, MAX, MIN 278* .. 279* .. Executable Statements .. 280* 281* Decode and check the input parameters. 282* 283 INFO = 0 284 ICTXT = DESCH( CTXT_ ) 285 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 286 NPROCS = NPROW*NPCOL 287 IF( NPROW.EQ.-1 ) INFO = -(600+CTXT_) 288 IF( INFO.EQ.0 ) THEN 289 WANTT = LSAME( JOB, 'S' ) 290 INITZ = LSAME( COMPZ, 'I' ) 291 WANTZ = INITZ .OR. LSAME( COMPZ, 'V' ) 292 LLDH = DESCH( LLD_ ) 293 LLDZ = DESCZ( LLD_ ) 294 NB = DESCH( MB_ ) 295 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 296* 297 IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN 298 INFO = -1 299 ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN 300 INFO = -2 301 ELSE IF( N.LT.0 ) THEN 302 INFO = -3 303 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN 304 INFO = -4 305 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN 306 INFO = -5 307 ELSEIF( DESCZ( CTXT_ ).NE.DESCH( CTXT_ ) ) THEN 308 INFO = -( 1000+CTXT_ ) 309 ELSEIF( DESCH( MB_ ).NE.DESCH( NB_ ) ) THEN 310 INFO = -( 700+NB_ ) 311 ELSEIF( DESCZ( MB_ ).NE.DESCZ( NB_ ) ) THEN 312 INFO = -( 1000+NB_ ) 313 ELSEIF( DESCH( MB_ ).NE.DESCZ( MB_ ) ) THEN 314 INFO = -( 1000+MB_ ) 315 ELSEIF( DESCH( MB_ ).LT.6 ) THEN 316 INFO = -( 700+NB_ ) 317 ELSEIF( DESCZ( MB_ ).LT.6 ) THEN 318 INFO = -( 1000+MB_ ) 319 ELSE 320 CALL CHK1MAT( N, 3, N, 3, 1, 1, DESCH, 7, INFO ) 321 IF( INFO.EQ.0 ) 322 $ CALL CHK1MAT( N, 3, N, 3, 1, 1, DESCZ, 11, INFO ) 323 IF( INFO.EQ.0 ) 324 $ CALL PCHK2MAT( N, 3, N, 3, 1, 1, DESCH, 7, N, 3, N, 3, 325 $ 1, 1, DESCZ, 11, 0, IWORK, IWORK, INFO ) 326 END IF 327 END IF 328* 329* Compute required workspace. 330* 331 CALL PDLAQR1( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI, 332 $ ILO, IHI, Z, DESCZ, WORK, -1, IWORK, -1, INFO ) 333 LWKOPT = WORK(1) 334 LIWKOPT = IWORK(1) 335 CALL PDLAQR0( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI, 336 $ ILO, IHI, Z, DESCZ, WORK, -1, IWORK, -1, INFO, 0 ) 337 IF( N.LT.NL ) THEN 338 HROWS = NUMROC( NL, NB, MYROW, DESCH(RSRC_), NPROW ) 339 HCOLS = NUMROC( NL, NB, MYCOL, DESCH(CSRC_), NPCOL ) 340 WORK(1) = WORK(1) + DBLE(2*HROWS*HCOLS) 341 END IF 342 LWKOPT = MAX( LWKOPT, WORK(1) ) 343 LIWKOPT = MAX( LIWKOPT, IWORK(1) ) 344 WORK(1) = LWKOPT 345 IWORK(1) = LIWKOPT 346* 347 IF( .NOT.LQUERY .AND. LWORK.LT.INT(LWKOPT) ) THEN 348 INFO = -13 349 ELSEIF( .NOT.LQUERY .AND. LIWORK.LT.LIWKOPT ) THEN 350 INFO = -15 351 END IF 352* 353 IF( INFO.NE.0 ) THEN 354* 355* Quick return in case of invalid argument. 356* 357 CALL PXERBLA( ICTXT, 'PDHSEQR', -INFO ) 358 RETURN 359* 360 ELSE IF( N.EQ.0 ) THEN 361* 362* Quick return in case N = 0; nothing to do. 363* 364 RETURN 365* 366 ELSE IF( LQUERY ) THEN 367* 368* Quick return in case of a workspace query. 369* 370 RETURN 371* 372 ELSE 373* 374* Copy eigenvalues isolated by PDGEBAL. 375* 376 DO 10 I = 1, ILO - 1 377 CALL INFOG2L( I, I, DESCH, NPROW, NPCOL, MYROW, MYCOL, II, 378 $ JJ, HRSRC, HCSRC ) 379 IF( MYROW.EQ.HRSRC .AND. MYCOL.EQ.HCSRC ) THEN 380 WR( I ) = H( (JJ-1)*LLDH + II ) 381 ELSE 382 WR( I ) = ZERO 383 END IF 384 WI( I ) = ZERO 385 10 CONTINUE 386 IF( ILO.GT.1 ) 387 $ CALL DGSUM2D( ICTXT, 'All', '1-Tree', ILO-1, 1, WR, N, -1, 388 $ -1 ) 389 DO 20 I = IHI + 1, N 390 CALL INFOG2L( I, I, DESCH, NPROW, NPCOL, MYROW, MYCOL, II, 391 $ JJ, HRSRC, HCSRC ) 392 IF( MYROW.EQ.HRSRC .AND. MYCOL.EQ.HCSRC ) THEN 393 WR( I ) = H( (JJ-1)*LLDH + II ) 394 ELSE 395 WR( I ) = ZERO 396 END IF 397 WI( I ) = ZERO 398 20 CONTINUE 399 IF( IHI.LT.N ) 400 $ CALL DGSUM2D( ICTXT, 'All', '1-Tree', N-IHI, 1, WR(IHI+1), 401 $ N, -1, -1 ) 402* 403* Initialize Z, if requested. 404* 405 IF( INITZ ) 406 $ CALL PDLASET( 'A', N, N, ZERO, ONE, Z, 1, 1, DESCZ ) 407* 408* Quick return if possible. 409* 410 NPROCS = NPROW*NPCOL 411 IF( ILO.EQ.IHI ) THEN 412 CALL INFOG2L( ILO, ILO, DESCH, NPROW, NPCOL, MYROW, 413 $ MYCOL, II, JJ, HRSRC, HCSRC ) 414 IF( MYROW.EQ.HRSRC .AND. MYCOL.EQ.HCSRC ) THEN 415 WR( ILO ) = H( (JJ-1)*LLDH + II ) 416 IF( NPROCS.GT.1 ) 417 $ CALL DGEBS2D( ICTXT, 'All', '1-Tree', 1, 1, WR(ILO), 418 $ 1 ) 419 ELSE 420 CALL DGEBR2D( ICTXT, 'All', '1-Tree', 1, 1, WR(ILO), 421 $ 1, HRSRC, HCSRC ) 422 END IF 423 WI( ILO ) = ZERO 424 RETURN 425 END IF 426* 427* PDLAQR1/PDLAQR0 crossover point. 428* 429 NH = IHI-ILO+1 430 NMIN = PILAENVX( ICTXT, 12, 'PDHSEQR', 431 $ JOB( : 1 ) // COMPZ( : 1 ), N, ILO, IHI, LWORK ) 432 NMIN = MAX( NTINY, NMIN ) 433* 434* PDLAQR0 for big matrices; PDLAQR1 for small ones. 435* 436 IF( (.NOT. CRSOVER .AND. NH.GT.NTINY) .OR. NH.GT.NMIN .OR. 437 $ DESCH(RSRC_).NE.0 .OR. DESCH(CSRC_).NE.0 ) THEN 438 CALL PDLAQR0( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI, 439 $ ILO, IHI, Z, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO, 440 $ 0 ) 441 IF( INFO.GT.0 .AND. ( DESCH(RSRC_).NE.0 .OR. 442 $ DESCH(CSRC_).NE.0 ) ) THEN 443* 444* A rare PDLAQR0 failure! PDLAQR1 sometimes succeeds 445* when PDLAQR0 fails. 446* 447 KBOT = INFO 448 CALL PDLAQR1( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, 449 $ WI, ILO, IHI, Z, DESCZ, WORK, LWORK, IWORK, 450 $ LIWORK, INFO ) 451 INFO = -7777 452 END IF 453 ELSE 454* 455* Small matrix. 456* 457 CALL PDLAQR1( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI, 458 $ ILO, IHI, Z, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO ) 459* 460 IF( INFO.GT.0 ) THEN 461* 462* A rare PDLAQR1 failure! PDLAQR0 sometimes succeeds 463* when PDLAQR1 fails. 464* 465 KBOT = INFO 466* 467 IF( N.GE.NL ) THEN 468* 469* Larger matrices have enough subdiagonal scratch 470* space to call PDLAQR0 directly. 471* 472 CALL PDLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, DESCH, 473 $ WR, WI, ILO, IHI, Z, DESCZ, WORK, LWORK, 474 $ IWORK, LIWORK, INFO, 0 ) 475 ELSE 476* 477* Tiny matrices don't have enough subdiagonal 478* scratch space to benefit from PDLAQR0. Hence, 479* tiny matrices must be copied into a larger 480* array before calling PDLAQR0. 481* 482 HROWS = NUMROC( NL, NB, MYROW, DESCH(RSRC_), NPROW ) 483 HCOLS = NUMROC( NL, NB, MYCOL, DESCH(CSRC_), NPCOL ) 484 CALL DESCINIT( DESCH2, NL, NL, NB, NB, DESCH(RSRC_), 485 $ DESCH(CSRC_), ICTXT, MAX(1, HROWS), INFO ) 486 CALL PDLACPY( 'All', N, N, H, 1, 1, DESCH, WORK, 1, 487 $ 1, DESCH2 ) 488 CALL PDELSET( WORK, N+1, N, DESCH2, ZERO ) 489 CALL PDLASET( 'All', NL, NL-N, ZERO, ZERO, WORK, 1, 490 $ N+1, DESCH2 ) 491 IPW = 1 + DESCH2(LLD_)*HCOLS 492 CALL PDLAQR0( WANTT, WANTZ, NL, ILO, KBOT, WORK, 493 $ DESCH2, WR, WI, ILO, IHI, Z, DESCZ, 494 $ WORK(IPW), LWORK-IPW+1, IWORK, 495 $ LIWORK, INFO, 0 ) 496 IF( WANTT .OR. INFO.NE.0 ) 497 $ CALL PDLACPY( 'All', N, N, WORK, 1, 1, DESCH2, 498 $ H, 1, 1, DESCH ) 499 END IF 500 INFO = -8888 501 END IF 502 END IF 503* 504* Clear out the trash, if necessary. 505* 506 IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 ) 507 $ CALL PDLASET( 'L', N-2, N-2, ZERO, ZERO, H, 3, 1, DESCH ) 508* 509* Force any 2-by-2 blocks to be complex conjugate pairs of 510* eigenvalues by removing false such blocks. 511* 512 DO 30 I = ILO, IHI-1 513 CALL PDELGET( 'All', ' ', TMP3, H, I+1, I, DESCH ) 514 IF( TMP3.NE.0.0D+00 ) THEN 515 CALL PDELGET( 'All', ' ', TMP1, H, I, I, DESCH ) 516 CALL PDELGET( 'All', ' ', TMP2, H, I, I+1, DESCH ) 517 CALL PDELGET( 'All', ' ', TMP4, H, I+1, I+1, DESCH ) 518 CALL DLANV2( TMP1, TMP2, TMP3, TMP4, DUM1, DUM2, DUM3, 519 $ DUM4, CS, SN ) 520 IF( TMP3.EQ.0.0D+00 ) THEN 521 IF( WANTT ) THEN 522 IF( I+2.LE.N ) 523 $ CALL PDROT( N-I-1, H, I, I+2, DESCH, 524 $ DESCH(M_), H, I+1, I+2, DESCH, DESCH(M_), 525 $ CS, SN, WORK, LWORK, INFO ) 526 CALL PDROT( I-1, H, 1, I, DESCH, 1, H, 1, I+1, 527 $ DESCH, 1, CS, SN, WORK, LWORK, INFO ) 528 END IF 529 IF( WANTZ ) THEN 530 CALL PDROT( N, Z, 1, I, DESCZ, 1, Z, 1, I+1, DESCZ, 531 $ 1, CS, SN, WORK, LWORK, INFO ) 532 END IF 533 CALL PDELSET( H, I, I, DESCH, TMP1 ) 534 CALL PDELSET( H, I, I+1, DESCH, TMP2 ) 535 CALL PDELSET( H, I+1, I, DESCH, TMP3 ) 536 CALL PDELSET( H, I+1, I+1, DESCH, TMP4 ) 537 END IF 538 END IF 539 30 CONTINUE 540* 541* Read out eigenvalues: first let all the processes compute the 542* eigenvalue inside their diagonal blocks in parallel, except for 543* the eigenvalue located next to a block border. After that, 544* compute all eigenvalues located next to the block borders. 545* Finally, do a global summation over WR and WI so that all 546* processors receive the result. 547* 548 DO 40 K = ILO, IHI 549 WR( K ) = ZERO 550 WI( K ) = ZERO 551 40 CONTINUE 552 NB = DESCH( MB_ ) 553* 554* Loop 50: extract eigenvalues from the blocks which are not laid 555* out across a border of the processor mesh, except for those 1x1 556* blocks on the border. 557* 558 PAIR = .FALSE. 559 DO 50 K = ILO, IHI 560 IF( .NOT. PAIR ) THEN 561 BORDER = MOD( K, NB ).EQ.0 .OR. ( K.NE.1 .AND. 562 $ MOD( K, NB ).EQ.1 ) 563 IF( .NOT. BORDER ) THEN 564 CALL INFOG2L( K, K, DESCH, NPROW, NPCOL, MYROW, 565 $ MYCOL, ILOC1, JLOC1, HRSRC1, HCSRC1 ) 566 IF( MYROW.EQ.HRSRC1 .AND. MYCOL.EQ.HCSRC1 ) THEN 567 ELEM1 = H((JLOC1-1)*LLDH+ILOC1) 568 IF( K.LT.N ) THEN 569 ELEM3 = H((JLOC1-1)*LLDH+ILOC1+1) 570 ELSE 571 ELEM3 = ZERO 572 END IF 573 IF( ELEM3.NE.ZERO ) THEN 574 ELEM2 = H((JLOC1)*LLDH+ILOC1) 575 ELEM4 = H((JLOC1)*LLDH+ILOC1+1) 576 CALL DLANV2( ELEM1, ELEM2, ELEM3, ELEM4, 577 $ WR( K ), WI( K ), WR( K+1 ), WI( K+1 ), 578 $ SN, CS ) 579 PAIR = .TRUE. 580 ELSE 581 IF( K.GT.1 ) THEN 582 TMP = H((JLOC1-2)*LLDH+ILOC1) 583 IF( TMP.NE.ZERO ) THEN 584 ELEM1 = H((JLOC1-2)*LLDH+ILOC1-1) 585 ELEM2 = H((JLOC1-1)*LLDH+ILOC1-1) 586 ELEM3 = H((JLOC1-2)*LLDH+ILOC1) 587 ELEM4 = H((JLOC1-1)*LLDH+ILOC1) 588 CALL DLANV2( ELEM1, ELEM2, ELEM3, 589 $ ELEM4, WR( K-1 ), WI( K-1 ), 590 $ WR( K ), WI( K ), SN, CS ) 591 ELSE 592 WR( K ) = ELEM1 593 END IF 594 ELSE 595 WR( K ) = ELEM1 596 END IF 597 END IF 598 END IF 599 END IF 600 ELSE 601 PAIR = .FALSE. 602 END IF 603 50 CONTINUE 604* 605* Loop 60: extract eigenvalues from the blocks which are laid 606* out across a border of the processor mesh. The processors are 607* numbered as below: 608* 609* 1 | 2 610* --+-- 611* 3 | 4 612* 613 DO 60 K = ICEIL(ILO,NB)*NB, IHI-1, NB 614 CALL INFOG2L( K, K, DESCH, NPROW, NPCOL, MYROW, MYCOL, 615 $ ILOC1, JLOC1, HRSRC1, HCSRC1 ) 616 CALL INFOG2L( K, K+1, DESCH, NPROW, NPCOL, MYROW, MYCOL, 617 $ ILOC2, JLOC2, HRSRC2, HCSRC2 ) 618 CALL INFOG2L( K+1, K, DESCH, NPROW, NPCOL, MYROW, MYCOL, 619 $ ILOC3, JLOC3, HRSRC3, HCSRC3 ) 620 CALL INFOG2L( K+1, K+1, DESCH, NPROW, NPCOL, MYROW, MYCOL, 621 $ ILOC4, JLOC4, HRSRC4, HCSRC4 ) 622 IF( MYROW.EQ.HRSRC2 .AND. MYCOL.EQ.HCSRC2 ) THEN 623 ELEM2 = H((JLOC2-1)*LLDH+ILOC2) 624 IF( HRSRC1.NE.HRSRC2 .OR. HCSRC1.NE.HCSRC2 ) 625 $ CALL DGESD2D( ICTXT, 1, 1, ELEM2, 1, HRSRC1, HCSRC1) 626 END IF 627 IF( MYROW.EQ.HRSRC3 .AND. MYCOL.EQ.HCSRC3 ) THEN 628 ELEM3 = H((JLOC3-1)*LLDH+ILOC3) 629 IF( HRSRC1.NE.HRSRC3 .OR. HCSRC1.NE.HCSRC3 ) 630 $ CALL DGESD2D( ICTXT, 1, 1, ELEM3, 1, HRSRC1, HCSRC1) 631 END IF 632 IF( MYROW.EQ.HRSRC4 .AND. MYCOL.EQ.HCSRC4 ) THEN 633 WORK(1) = H((JLOC4-1)*LLDH+ILOC4) 634 IF( K+1.LT.N ) THEN 635 WORK(2) = H((JLOC4-1)*LLDH+ILOC4+1) 636 ELSE 637 WORK(2) = ZERO 638 END IF 639 IF( HRSRC1.NE.HRSRC4 .OR. HCSRC1.NE.HCSRC4 ) 640 $ CALL DGESD2D( ICTXT, 2, 1, WORK, 2, HRSRC1, HCSRC1 ) 641 END IF 642 IF( MYROW.EQ.HRSRC1 .AND. MYCOL.EQ.HCSRC1 ) THEN 643 ELEM1 = H((JLOC1-1)*LLDH+ILOC1) 644 IF( HRSRC1.NE.HRSRC2 .OR. HCSRC1.NE.HCSRC2 ) 645 $ CALL DGERV2D( ICTXT, 1, 1, ELEM2, 1, HRSRC2, HCSRC2) 646 IF( HRSRC1.NE.HRSRC3 .OR. HCSRC1.NE.HCSRC3 ) 647 $ CALL DGERV2D( ICTXT, 1, 1, ELEM3, 1, HRSRC3, HCSRC3) 648 IF( HRSRC1.NE.HRSRC4 .OR. HCSRC1.NE.HCSRC4 ) 649 $ CALL DGERV2D( ICTXT, 2, 1, WORK, 2, HRSRC4, HCSRC4 ) 650 ELEM4 = WORK(1) 651 ELEM5 = WORK(2) 652 IF( ELEM5.EQ.ZERO ) THEN 653 IF( WR( K ).EQ.ZERO .AND. WI( K ).EQ.ZERO ) THEN 654 CALL DLANV2( ELEM1, ELEM2, ELEM3, ELEM4, WR( K ), 655 $ WI( K ), WR( K+1 ), WI( K+1 ), SN, CS ) 656 ELSEIF( WR( K+1 ).EQ.ZERO .AND. WI( K+1 ).EQ.ZERO ) 657 $ THEN 658 WR( K+1 ) = ELEM4 659 END IF 660 ELSEIF( WR( K ).EQ.ZERO .AND. WI( K ).EQ.ZERO ) 661 $ THEN 662 WR( K ) = ELEM1 663 END IF 664 END IF 665 60 CONTINUE 666* 667 IF( NPROCS.GT.1 ) THEN 668 CALL DGSUM2D( ICTXT, 'All', ' ', IHI-ILO+1, 1, WR(ILO), N, 669 $ -1, -1 ) 670 CALL DGSUM2D( ICTXT, 'All', ' ', IHI-ILO+1, 1, WI(ILO), N, 671 $ -1, -1 ) 672 END IF 673* 674 END IF 675* 676 WORK(1) = LWKOPT 677 IWORK(1) = LIWKOPT 678 RETURN 679* 680* End of PDHSEQR 681* 682 END 683