|
SRC\pdlaconsb.f |
|
| #lines: 574 size: 20 Kb creation: 18/01/2006 23:36:04 last modification: 08/05/2008 18:37:53 attribute: ARCH Find Reload | |
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54: 55: 56: 57: 58: 59: 60: 61: 62: 63: 64: 65: 66: 67: 68: 69: 70: 71: 72: 73: 74: 75: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: 116: 117: 118: 119: 120: 121: 122: 123: 124: 125: 126: 127: 128: 129: 130: 131: 132: 133: 134: 135: 136: 137: 138: 139: 140: 141: 142: 143: 144: 145: 146: 147: 148: 149: 150: 151: 152: 153: 154: 155: 156: 157: 158: 159: 160: 161: 162: 163: 164: 165: 166: 167: 168: 169: 170: 171: 172: 173: 174: 175: 176: 177: 178: 179: 180: 181: 182: 183: 184: 185: 186: 187: 188: 189: 190: 191: 192: 193: 194: 195: 196: 197: 198: 199: 200: 201: 202: 203: 204: 205: 206: 207: 208: 209: 210: 211: 212: 213: 214: 215: 216: 217: 218: 219: 220: 221: 222: 223: 224: 225: 226: 227: 228: 229: 230: 231: 232: 233: 234: 235: 236: 237: 238: 239: 240: 241: 242: 243: 244: 245: 246: 247: 248: 249: 250: 251: 252: 253: 254: 255: 256: 257: 258: 259: 260: 261: 262: 263: 264: 265: 266: 267: 268: 269: 270: 271: 272: 273: 274: 275: 276: 277: 278: 279: 280: 281: 282: 283: 284: 285: 286: 287: 288: 289: 290: 291: 292: 293: 294: 295: 296: 297: 298: 299: 300: 301: 302: 303: 304: 305: 306: 307: 308: 309: 310: 311: 312: 313: 314: 315: 316: 317: 318: 319: 320: 321: 322: 323: 324: 325: 326: 327: 328: 329: 330: 331: 332: 333: 334: 335: 336: 337: 338: 339: 340: 341: 342: 343: 344: 345: 346: 347: 348: 349: 350: 351: 352: 353: 354: 355: 356: 357: 358: 359: 360: 361: 362: 363: 364: 365: 366: 367: 368: 369: 370: 371: 372: 373: 374: 375: 376: 377: 378: 379: 380: 381: 382: 383: 384: 385: 386: 387: 388: 389: 390: 391: 392: 393: 394: 395: 396: 397: 398: 399: 400: 401: 402: 403: 404: 405: 406: 407: 408: 409: 410: 411: 412: 413: 414: 415: 416: 417: 418: 419: 420: 421: 422: 423: 424: 425: 426: 427: 428: 429: 430: 431: 432: 433: 434: 435: 436: 437: 438: 439: 440: 441: 442: 443: 444: 445: 446: 447: 448: 449: 450: 451: 452: 453: 454: 455: 456: 457: 458: 459: 460: 461: 462: 463: 464: 465: 466: 467: 468: 469: 470: 471: 472: 473: 474: 475: 476: 477: 478: 479: 480: 481: 482: 483: 484: 485: 486: 487: 488: 489: 490: 491: 492: 493: 494: 495: 496: 497: 498: 499: 500: 501: 502: 503: 504: 505: 506: 507: 508: 509: 510: 511: 512: 513: 514: 515: 516: 517: 518: 519: 520: 521: 522: 523: 524: 525: 526: 527: 528: 529: 530: 531: 532: 533: 534: 535: 536: 537: 538: 539: 540: 541: 542: 543: 544: 545: 546: 547: 548: 549: 550: 551: 552: 553: 554: 555: 556: 557: 558: 559: 560: 561: 562: 563: 564: 565: 566: 567: 568: 569: 570: 571: 572: 573: 574: |
SUBROUTINE PDLACONSB( A, DESCA, I, L, M, H44, H33, H43H34, BUF,
$ LWORK )
*
* -- ScaLAPACK routine (version 1.7) --
* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
* and University of California, Berkeley.
* May 1, 1997
*
* .. Scalar Arguments ..
INTEGER I, L, LWORK, M
DOUBLE PRECISION H33, H43H34, H44
* ..
* .. Array Arguments ..
INTEGER DESCA( * )
DOUBLE PRECISION A( * ), BUF( * )
* ..
*
* Purpose
* =======
*
* PDLACONSB looks for two consecutive small subdiagonal elements by
* seeing the effect of starting a double shift QR iteration
* given by H44, H33, & H43H34 and see if this would make a
* subdiagonal negligible.
*
* Notes
* =====
*
* Each global data object is described by an associated description
* vector. This vector stores the information required to establish
* the mapping between an object element and its corresponding process
* and memory location.
*
* Let A be a generic term for any 2D block cyclicly distributed array.
* Such a global array has an associated description vector DESCA.
* In the following comments, the character _ should be read as
* "of the global array".
*
* NOTATION STORED IN EXPLANATION
* --------------- -------------- --------------------------------------
* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
* DTYPE_A = 1.
* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
* the BLACS process grid A is distribu-
* ted over. The context itself is glo-
* bal, but the handle (the integer
* value) may vary.
* M_A (global) DESCA( M_ ) The number of rows in the global
* array A.
* N_A (global) DESCA( N_ ) The number of columns in the global
* array A.
* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
* the rows of the array.
* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
* the columns of the array.
* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
* row of the array A is distributed.
* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
* first column of the array A is
* distributed.
* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
* array. LLD_A >= MAX(1,LOCr(M_A)).
*
* Let K be the number of rows or columns of a distributed matrix,
* and assume that its process grid has dimension p x q.
* LOCr( K ) denotes the number of elements of K that a process
* would receive if K were distributed over the p processes of its
* process column.
* Similarly, LOCc( K ) denotes the number of elements of K that a
* process would receive if K were distributed over the q processes of
* its process row.
* The values of LOCr() and LOCc() may be determined via a call to the
* ScaLAPACK tool function, NUMROC:
* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
* An upper bound for these quantities may be computed by:
* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
* Arguments
* =========
*
* A (global input) DOUBLE PRECISION array, dimension
* (DESCA(LLD_),*)
* On entry, the Hessenberg matrix whose tridiagonal part is
* being scanned.
* Unchanged on exit.
*
* DESCA (global and local input) INTEGER array of dimension DLEN_.
* The array descriptor for the distributed matrix A.
*
* I (global input) INTEGER
* The global location of the bottom of the unreduced
* submatrix of A.
* Unchanged on exit.
*
* L (global input) INTEGER
* The global location of the top of the unreduced submatrix
* of A.
* Unchanged on exit.
*
* M (global output) INTEGER
* On exit, this yields the starting location of the QR double
* shift. This will satisfy: L <= M <= I-2.
*
* H44
* H33
* H43H34 (global input) DOUBLE PRECISION
* These three values are for the double shift QR iteration.
*
* BUF (local output) DOUBLE PRECISION array of size LWORK.
*
* LWORK (global input) INTEGER
* On exit, LWORK is the size of the work buffer.
* This must be at least 7*Ceil( Ceil( (I-L)/HBL ) /
* LCM(NPROW,NPCOL) )
* Here LCM is least common multiple, and NPROWxNPCOL is the
* logical grid size.
*
* Logic:
* ======
*
* Two consecutive small subdiagonal elements will stall
* convergence of a double shift if their product is small
* relatively even if each is not very small. Thus it is
* necessary to scan the "tridiagonal portion of the matrix." In
* the LAPACK algorithm DLAHQR, a loop of M goes from I-2 down to
* L and examines
* H(m,m),H(m+1,m+1),H(m+1,m),H(m,m+1),H(m-1,m-1),H(m,m-1), and
* H(m+2,m-1). Since these elements may be on separate
* processors, the first major loop (10) goes over the tridiagonal
* and has each node store whatever values of the 7 it has that
* the node owning H(m,m) does not. This will occur on a border
* and can happen in no more than 3 locations per block assuming
* square blocks. There are 5 buffers that each node stores these
* values: a buffer to send diagonally down and right, a buffer
* to send up, a buffer to send left, a buffer to send diagonally
* up and left and a buffer to send right. Each of these buffers
* is actually stored in one buffer BUF where BUF(ISTR1+1) starts
* the first buffer, BUF(ISTR2+1) starts the second, etc.. After
* the values are stored, if there are any values that a node
* needs, they will be sent and received. Then the next major
* loop passes over the data and searches for two consecutive
* small subdiagonals.
*
* Notes:
*
* This routine does a global maximum and must be called by all
* processes.
*
*
* Implemented by: G. Henry, November 17, 1996
*
* =====================================================================
*
* .. Parameters ..
INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
$ LLD_, MB_, M_, NB_, N_, RSRC_
PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
$ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
$ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
* ..
* .. Local Scalars ..
INTEGER CONTXT, DOWN, HBL, IBUF1, IBUF2, IBUF3, IBUF4,
$ IBUF5, ICOL1, II, IRCV1, IRCV2, IRCV3, IRCV4,
$ IRCV5, IROW1, ISRC, ISTR1, ISTR2, ISTR3, ISTR4,
$ ISTR5, JJ, JSRC, LDA, LEFT, MODKM1, MYCOL,
$ MYROW, NPCOL, NPROW, NUM, RIGHT, UP
DOUBLE PRECISION H00, H10, H11, H12, H21, H22, H33S, H44S, S,
$ TST1, ULP, V1, V2, V3
* ..
* .. External Functions ..
INTEGER ILCM
DOUBLE PRECISION PDLAMCH
EXTERNAL ILCM, PDLAMCH
* ..
* .. External Subroutines ..
EXTERNAL BLACS_GRIDINFO, DGERV2D, DGESD2D, IGAMX2D,
$ INFOG2L, PXERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MOD
* ..
* .. Executable Statements ..
*
HBL = DESCA( MB_ )
CONTXT = DESCA( CTXT_ )
LDA = DESCA( LLD_ )
ULP = PDLAMCH( CONTXT, 'PRECISION' )
CALL BLACS_GRIDINFO( CONTXT, NPROW, NPCOL, MYROW, MYCOL )
LEFT = MOD( MYCOL+NPCOL-1, NPCOL )
RIGHT = MOD( MYCOL+1, NPCOL )
UP = MOD( MYROW+NPROW-1, NPROW )
DOWN = MOD( MYROW+1, NPROW )
NUM = NPROW*NPCOL
*
* BUFFER1 starts at BUF(ISTR1+1) and will contain IBUF1 elements
* BUFFER2 starts at BUF(ISTR2+1) and will contain IBUF2 elements
* BUFFER3 starts at BUF(ISTR3+1) and will contain IBUF3 elements
* BUFFER4 starts at BUF(ISTR4+1) and will contain IBUF4 elements
* BUFFER5 starts at BUF(ISTR5+1) and will contain IBUF5 elements
*
ISTR1 = 0
ISTR2 = ( ( I-L-1 ) / HBL )
IF( ISTR2*HBL.LT.( I-L-1 ) )
$ ISTR2 = ISTR2 + 1
II = ISTR2 / ILCM( NPROW, NPCOL )
IF( II*ILCM( NPROW, NPCOL ).LT.ISTR2 ) THEN
ISTR2 = II + 1
ELSE
ISTR2 = II
END IF
IF( LWORK.LT.7*ISTR2 ) THEN
CALL PXERBLA( CONTXT, 'PDLACONSB', 10 )
RETURN
END IF
ISTR3 = 3*ISTR2
ISTR4 = ISTR3 + ISTR2
ISTR5 = ISTR3 + ISTR3
CALL INFOG2L( I-2, I-2, DESCA, NPROW, NPCOL, MYROW, MYCOL, IROW1,
$ ICOL1, II, JJ )
MODKM1 = MOD( I-3+HBL, HBL )
*
* Copy our relevant pieces of triadiagonal that we owe into
* 5 buffers to send to whomever owns H(M,M) as M moves diagonally
* up the tridiagonal
*
IBUF1 = 0
IBUF2 = 0
IBUF3 = 0
IBUF4 = 0
IBUF5 = 0
IRCV1 = 0
IRCV2 = 0
IRCV3 = 0
IRCV4 = 0
IRCV5 = 0
DO 10 M = I - 2, L, -1
IF( ( MODKM1.EQ.0 ) .AND. ( DOWN.EQ.II ) .AND.
$ ( RIGHT.EQ.JJ ) .AND. ( M.GT.L ) ) THEN
*
* We must pack H(M-1,M-1) and send it diagonal down
*
IF( ( DOWN.NE.MYROW ) .OR. ( RIGHT.NE.MYCOL ) ) THEN
CALL INFOG2L( M-1, M-1, DESCA, NPROW, NPCOL, MYROW,
$ MYCOL, IROW1, ICOL1, ISRC, JSRC )
IBUF1 = IBUF1 + 1
BUF( ISTR1+IBUF1 ) = A( ( ICOL1-1 )*LDA+IROW1 )
END IF
END IF
IF( ( MODKM1.EQ.0 ) .AND. ( MYROW.EQ.II ) .AND.
$ ( RIGHT.EQ.JJ ) .AND. ( M.GT.L ) ) THEN
*
* We must pack H(M ,M-1) and send it right
*
IF( NPCOL.GT.1 ) THEN
CALL INFOG2L( M, M-1, DESCA, NPROW, NPCOL, MYROW, MYCOL,
$ IROW1, ICOL1, ISRC, JSRC )
IBUF5 = IBUF5 + 1
BUF( ISTR5+IBUF5 ) = A( ( ICOL1-1 )*LDA+IROW1 )
END IF
END IF
IF( ( MODKM1.EQ.HBL-1 ) .AND. ( UP.EQ.II ) .AND.
$ ( MYCOL.EQ.JJ ) ) THEN
*
* We must pack H(M+1,M) and send it up
*
IF( NPROW.GT.1 ) THEN
CALL INFOG2L( M+1, M, DESCA, NPROW, NPCOL, MYROW, MYCOL,
$ IROW1, ICOL1, ISRC, JSRC )
IBUF2 = IBUF2 + 1
BUF( ISTR2+IBUF2 ) = A( ( ICOL1-1 )*LDA+IROW1 )
END IF
END IF
IF( ( MODKM1.EQ.HBL-1 ) .AND. ( MYROW.EQ.II ) .AND.
$ ( LEFT.EQ.JJ ) ) THEN
*
* We must pack H(M ,M+1) and send it left
*
IF( NPCOL.GT.1 ) THEN
CALL INFOG2L( M, M+1, DESCA, NPROW, NPCOL, MYROW, MYCOL,
$ IROW1, ICOL1, ISRC, JSRC )
IBUF3 = IBUF3 + 1
BUF( ISTR3+IBUF3 ) = A( ( ICOL1-1 )*LDA+IROW1 )
END IF
END IF
IF( ( MODKM1.EQ.HBL-1 ) .AND. ( UP.EQ.II ) .AND.
$ ( LEFT.EQ.JJ ) ) THEN
*
* We must pack H(M+1,M+1) & H(M+2,M+1) and send it
* diagonally up
*
IF( ( UP.NE.MYROW ) .OR. ( LEFT.NE.MYCOL ) ) THEN
CALL INFOG2L( M+1, M+1, DESCA, NPROW, NPCOL, MYROW,
$ MYCOL, IROW1, ICOL1, ISRC, JSRC )
IBUF4 = IBUF4 + 2
BUF( ISTR4+IBUF4-1 ) = A( ( ICOL1-1 )*LDA+IROW1 )
BUF( ISTR4+IBUF4 ) = A( ( ICOL1-1 )*LDA+IROW1+1 )
END IF
END IF
IF( ( MODKM1.EQ.HBL-2 ) .AND. ( UP.EQ.II ) .AND.
$ ( MYCOL.EQ.JJ ) ) THEN
*
* We must pack H(M+2,M+1) and send it up
*
IF( NPROW.GT.1 ) THEN
CALL INFOG2L( M+2, M+1, DESCA, NPROW, NPCOL, MYROW,
$ MYCOL, IROW1, ICOL1, ISRC, JSRC )
IBUF2 = IBUF2 + 1
BUF( ISTR2+IBUF2 ) = A( ( ICOL1-1 )*LDA+IROW1 )
END IF
END IF
*
* Add up the receives
*
IF( ( MYROW.EQ.II ) .AND. ( MYCOL.EQ.JJ ) ) THEN
IF( ( MODKM1.EQ.0 ) .AND. ( M.GT.L ) .AND.
$ ( ( NPROW.GT.1 ) .OR. ( NPCOL.GT.1 ) ) ) THEN
*
* We must receive H(M-1,M-1) from diagonal up
*
IRCV1 = IRCV1 + 1
END IF
IF( ( MODKM1.EQ.0 ) .AND. ( NPCOL.GT.1 ) .AND. ( M.GT.L ) )
$ THEN
*
* We must receive H(M ,M-1) from left
*
IRCV5 = IRCV5 + 1
END IF
IF( ( MODKM1.EQ.HBL-1 ) .AND. ( NPROW.GT.1 ) ) THEN
*
* We must receive H(M+1,M ) from down
*
IRCV2 = IRCV2 + 1
END IF
IF( ( MODKM1.EQ.HBL-1 ) .AND. ( NPCOL.GT.1 ) ) THEN
*
* We must receive H(M ,M+1) from right
*
IRCV3 = IRCV3 + 1
END IF
IF( ( MODKM1.EQ.HBL-1 ) .AND.
$ ( ( NPROW.GT.1 ) .OR. ( NPCOL.GT.1 ) ) ) THEN
*
* We must receive H(M+1:M+2,M+1) from diagonal down
*
IRCV4 = IRCV4 + 2
END IF
IF( ( MODKM1.EQ.HBL-2 ) .AND. ( NPROW.GT.1 ) ) THEN
*
* We must receive H(M+2,M+1) from down
*
IRCV2 = IRCV2 + 1
END IF
END IF
*
* Possibly change owners (occurs only when MOD(M-1,HBL) = 0)
*
IF( MODKM1.EQ.0 ) THEN
II = II - 1
JJ = JJ - 1
IF( II.LT.0 )
$ II = NPROW - 1
IF( JJ.LT.0 )
$ JJ = NPCOL - 1
END IF
MODKM1 = MODKM1 - 1
IF( MODKM1.LT.0 )
$ MODKM1 = HBL - 1
10 CONTINUE
*
*
* Send data on to the appropriate node if there is any data to send
*
IF( IBUF1.GT.0 ) THEN
CALL DGESD2D( CONTXT, IBUF1, 1, BUF( ISTR1+1 ), IBUF1, DOWN,
$ RIGHT )
END IF
IF( IBUF2.GT.0 ) THEN
CALL DGESD2D( CONTXT, IBUF2, 1, BUF( ISTR2+1 ), IBUF2, UP,
$ MYCOL )
END IF
IF( IBUF3.GT.0 ) THEN
CALL DGESD2D( CONTXT, IBUF3, 1, BUF( ISTR3+1 ), IBUF3, MYROW,
$ LEFT )
END IF
IF( IBUF4.GT.0 ) THEN
CALL DGESD2D( CONTXT, IBUF4, 1, BUF( ISTR4+1 ), IBUF4, UP,
$ LEFT )
END IF
IF( IBUF5.GT.0 ) THEN
CALL DGESD2D( CONTXT, IBUF5, 1, BUF( ISTR5+1 ), IBUF5, MYROW,
$ RIGHT )
END IF
*
* Receive appropriate data if there is any
*
IF( IRCV1.GT.0 ) THEN
CALL DGERV2D( CONTXT, IRCV1, 1, BUF( ISTR1+1 ), IRCV1, UP,
$ LEFT )
END IF
IF( IRCV2.GT.0 ) THEN
CALL DGERV2D( CONTXT, IRCV2, 1, BUF( ISTR2+1 ), IRCV2, DOWN,
$ MYCOL )
END IF
IF( IRCV3.GT.0 ) THEN
CALL DGERV2D( CONTXT, IRCV3, 1, BUF( ISTR3+1 ), IRCV3, MYROW,
$ RIGHT )
END IF
IF( IRCV4.GT.0 ) THEN
CALL DGERV2D( CONTXT, IRCV4, 1, BUF( ISTR4+1 ), IRCV4, DOWN,
$ RIGHT )
END IF
IF( IRCV5.GT.0 ) THEN
CALL DGERV2D( CONTXT, IRCV5, 1, BUF( ISTR5+1 ), IRCV5, MYROW,
$ LEFT )
END IF
*
* Start main loop
*
IBUF1 = 0
IBUF2 = 0
IBUF3 = 0
IBUF4 = 0
IBUF5 = 0
CALL INFOG2L( I-2, I-2, DESCA, NPROW, NPCOL, MYROW, MYCOL, IROW1,
$ ICOL1, II, JJ )
MODKM1 = MOD( I-3+HBL, HBL )
IF( ( MYROW.EQ.II ) .AND. ( MYCOL.EQ.JJ ) .AND.
$ ( MODKM1.NE.HBL-1 ) ) THEN
CALL INFOG2L( I-2, I-1, DESCA, NPROW, NPCOL, MYROW, MYCOL,
$ IROW1, ICOL1, ISRC, JSRC )
END IF
*
* Look for two consecutive small subdiagonal elements.
*
DO 20 M = I - 2, L, -1
*
* Determine the effect of starting the double-shift QR
* iteration at row M, and see if this would make H(M,M-1)
* negligible.
*
IF( ( MYROW.EQ.II ) .AND. ( MYCOL.EQ.JJ ) ) THEN
IF( MODKM1.EQ.0 ) THEN
H22 = A( ( ICOL1-1 )*LDA+IROW1+1 )
H11 = A( ( ICOL1-2 )*LDA+IROW1 )
V3 = A( ( ICOL1-1 )*LDA+IROW1+2 )
H21 = A( ( ICOL1-2 )*LDA+IROW1+1 )
H12 = A( ( ICOL1-1 )*LDA+IROW1 )
IF( M.GT.L ) THEN
IF( NUM.GT.1 ) THEN
IBUF1 = IBUF1 + 1
H00 = BUF( ISTR1+IBUF1 )
ELSE
H00 = A( ( ICOL1-3 )*LDA+IROW1-1 )
END IF
IF( NPCOL.GT.1 ) THEN
IBUF5 = IBUF5 + 1
H10 = BUF( ISTR5+IBUF5 )
ELSE
H10 = A( ( ICOL1-3 )*LDA+IROW1 )
END IF
END IF
END IF
IF( MODKM1.EQ.HBL-1 ) THEN
CALL INFOG2L( M, M, DESCA, NPROW, NPCOL, MYROW, MYCOL,
$ IROW1, ICOL1, ISRC, JSRC )
H11 = A( ( ICOL1-1 )*LDA+IROW1 )
IF( NUM.GT.1 ) THEN
IBUF4 = IBUF4 + 2
H22 = BUF( ISTR4+IBUF4-1 )
V3 = BUF( ISTR4+IBUF4 )
ELSE
H22 = A( ICOL1*LDA+IROW1+1 )
V3 = A( ( ICOL1+1 )*LDA+IROW1+1 )
END IF
IF( NPROW.GT.1 ) THEN
IBUF2 = IBUF2 + 1
H21 = BUF( ISTR2+IBUF2 )
ELSE
H21 = A( ( ICOL1-1 )*LDA+IROW1+1 )
END IF
IF( NPCOL.GT.1 ) THEN
IBUF3 = IBUF3 + 1
H12 = BUF( ISTR3+IBUF3 )
ELSE
H12 = A( ICOL1*LDA+IROW1 )
END IF
IF( M.GT.L ) THEN
H00 = A( ( ICOL1-2 )*LDA+IROW1-1 )
H10 = A( ( ICOL1-2 )*LDA+IROW1 )
END IF
*
* Adjust ICOL1 for next iteration where MODKM1=HBL-2
*
ICOL1 = ICOL1 + 1
END IF
IF( MODKM1.EQ.HBL-2 ) THEN
H22 = A( ( ICOL1-1 )*LDA+IROW1+1 )
H11 = A( ( ICOL1-2 )*LDA+IROW1 )
IF( NPROW.GT.1 ) THEN
IBUF2 = IBUF2 + 1
V3 = BUF( ISTR2+IBUF2 )
ELSE
V3 = A( ( ICOL1-1 )*LDA+IROW1+2 )
END IF
H21 = A( ( ICOL1-2 )*LDA+IROW1+1 )
H12 = A( ( ICOL1-1 )*LDA+IROW1 )
IF( M.GT.L ) THEN
H00 = A( ( ICOL1-3 )*LDA+IROW1-1 )
H10 = A( ( ICOL1-3 )*LDA+IROW1 )
END IF
END IF
IF( ( MODKM1.LT.HBL-2 ) .AND. ( MODKM1.GT.0 ) ) THEN
H22 = A( ( ICOL1-1 )*LDA+IROW1+1 )
H11 = A( ( ICOL1-2 )*LDA+IROW1 )
V3 = A( ( ICOL1-1 )*LDA+IROW1+2 )
H21 = A( ( ICOL1-2 )*LDA+IROW1+1 )
H12 = A( ( ICOL1-1 )*LDA+IROW1 )
IF( M.GT.L ) THEN
H00 = A( ( ICOL1-3 )*LDA+IROW1-1 )
H10 = A( ( ICOL1-3 )*LDA+IROW1 )
END IF
END IF
H44S = H44 - H11
H33S = H33 - H11
V1 = ( H33S*H44S-H43H34 ) / H21 + H12
V2 = H22 - H11 - H33S - H44S
S = ABS( V1 ) + ABS( V2 ) + ABS( V3 )
V1 = V1 / S
V2 = V2 / S
V3 = V3 / S
IF( M.EQ.L )
$ GO TO 30
TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) )
IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) ).LE.ULP*TST1 )
$ GO TO 30
*
* Slide indices diagonally up one for next iteration
*
IROW1 = IROW1 - 1
ICOL1 = ICOL1 - 1
END IF
IF( M.EQ.L ) THEN
*
* Stop regardless of which node we are
*
GO TO 30
END IF
*
* Possibly change owners if on border
*
IF( MODKM1.EQ.0 ) THEN
II = II - 1
JJ = JJ - 1
IF( II.LT.0 )
$ II = NPROW - 1
IF( JJ.LT.0 )
$ JJ = NPCOL - 1
END IF
MODKM1 = MODKM1 - 1
IF( MODKM1.LT.0 )
$ MODKM1 = HBL - 1
20 CONTINUE
30 CONTINUE
*
CALL IGAMX2D( CONTXT, 'ALL', ' ', 1, 1, M, 1, L, L, -1, -1, -1 )
*
RETURN
*
* End of PDLACONSB
*
END
|