|
SRC\pclattrs.f |
|
| #lines: 1213 size: 44 Kb creation: 18/01/2006 23:36:04 last modification: 08/05/2008 18:37:47 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: 575: 576: 577: 578: 579: 580: 581: 582: 583: 584: 585: 586: 587: 588: 589: 590: 591: 592: 593: 594: 595: 596: 597: 598: 599: 600: 601: 602: 603: 604: 605: 606: 607: 608: 609: 610: 611: 612: 613: 614: 615: 616: 617: 618: 619: 620: 621: 622: 623: 624: 625: 626: 627: 628: 629: 630: 631: 632: 633: 634: 635: 636: 637: 638: 639: 640: 641: 642: 643: 644: 645: 646: 647: 648: 649: 650: 651: 652: 653: 654: 655: 656: 657: 658: 659: 660: 661: 662: 663: 664: 665: 666: 667: 668: 669: 670: 671: 672: 673: 674: 675: 676: 677: 678: 679: 680: 681: 682: 683: 684: 685: 686: 687: 688: 689: 690: 691: 692: 693: 694: 695: 696: 697: 698: 699: 700: 701: 702: 703: 704: 705: 706: 707: 708: 709: 710: 711: 712: 713: 714: 715: 716: 717: 718: 719: 720: 721: 722: 723: 724: 725: 726: 727: 728: 729: 730: 731: 732: 733: 734: 735: 736: 737: 738: 739: 740: 741: 742: 743: 744: 745: 746: 747: 748: 749: 750: 751: 752: 753: 754: 755: 756: 757: 758: 759: 760: 761: 762: 763: 764: 765: 766: 767: 768: 769: 770: 771: 772: 773: 774: 775: 776: 777: 778: 779: 780: 781: 782: 783: 784: 785: 786: 787: 788: 789: 790: 791: 792: 793: 794: 795: 796: 797: 798: 799: 800: 801: 802: 803: 804: 805: 806: 807: 808: 809: 810: 811: 812: 813: 814: 815: 816: 817: 818: 819: 820: 821: 822: 823: 824: 825: 826: 827: 828: 829: 830: 831: 832: 833: 834: 835: 836: 837: 838: 839: 840: 841: 842: 843: 844: 845: 846: 847: 848: 849: 850: 851: 852: 853: 854: 855: 856: 857: 858: 859: 860: 861: 862: 863: 864: 865: 866: 867: 868: 869: 870: 871: 872: 873: 874: 875: 876: 877: 878: 879: 880: 881: 882: 883: 884: 885: 886: 887: 888: 889: 890: 891: 892: 893: 894: 895: 896: 897: 898: 899: 900: 901: 902: 903: 904: 905: 906: 907: 908: 909: 910: 911: 912: 913: 914: 915: 916: 917: 918: 919: 920: 921: 922: 923: 924: 925: 926: 927: 928: 929: 930: 931: 932: 933: 934: 935: 936: 937: 938: 939: 940: 941: 942: 943: 944: 945: 946: 947: 948: 949: 950: 951: 952: 953: 954: 955: 956: 957: 958: 959: 960: 961: 962: 963: 964: 965: 966: 967: 968: 969: 970: 971: 972: 973: 974: 975: 976: 977: 978: 979: 980: 981: 982: 983: 984: 985: 986: 987: 988: 989: 990: 991: 992: 993: 994: 995: 996: 997: 998: 999: 1000: 1001: 1002: 1003: 1004: 1005: 1006: 1007: 1008: 1009: 1010: 1011: 1012: 1013: 1014: 1015: 1016: 1017: 1018: 1019: 1020: 1021: 1022: 1023: 1024: 1025: 1026: 1027: 1028: 1029: 1030: 1031: 1032: 1033: 1034: 1035: 1036: 1037: 1038: 1039: 1040: 1041: 1042: 1043: 1044: 1045: 1046: 1047: 1048: 1049: 1050: 1051: 1052: 1053: 1054: 1055: 1056: 1057: 1058: 1059: 1060: 1061: 1062: 1063: 1064: 1065: 1066: 1067: 1068: 1069: 1070: 1071: 1072: 1073: 1074: 1075: 1076: 1077: 1078: 1079: 1080: 1081: 1082: 1083: 1084: 1085: 1086: 1087: 1088: 1089: 1090: 1091: 1092: 1093: 1094: 1095: 1096: 1097: 1098: 1099: 1100: 1101: 1102: 1103: 1104: 1105: 1106: 1107: 1108: 1109: 1110: 1111: 1112: 1113: 1114: 1115: 1116: 1117: 1118: 1119: 1120: 1121: 1122: 1123: 1124: 1125: 1126: 1127: 1128: 1129: 1130: 1131: 1132: 1133: 1134: 1135: 1136: 1137: 1138: 1139: 1140: 1141: 1142: 1143: 1144: 1145: 1146: 1147: 1148: 1149: 1150: 1151: 1152: 1153: 1154: 1155: 1156: 1157: 1158: 1159: 1160: 1161: 1162: 1163: 1164: 1165: 1166: 1167: 1168: 1169: 1170: 1171: 1172: 1173: 1174: 1175: 1176: 1177: 1178: 1179: 1180: 1181: 1182: 1183: 1184: 1185: 1186: 1187: 1188: 1189: 1190: 1191: 1192: 1193: 1194: 1195: 1196: 1197: 1198: 1199: 1200: 1201: 1202: 1203: 1204: 1205: 1206: 1207: 1208: 1209: 1210: 1211: 1212: 1213: |
SUBROUTINE PCLATTRS( UPLO, TRANS, DIAG, NORMIN, N, A, IA, JA,
$ DESCA, X, IX, JX, DESCX, SCALE, CNORM, INFO )
*
* -- ScaLAPACK routine (version 1.7) --
* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
* and University of California, Berkeley.
* July 31, 2001
*
* .. Scalar Arguments ..
CHARACTER DIAG, NORMIN, TRANS, UPLO
INTEGER IA, INFO, IX, JA, JX, N
REAL SCALE
* ..
* .. Array Arguments ..
INTEGER DESCA( * ), DESCX( * )
REAL CNORM( * )
COMPLEX A( * ), X( * )
* ..
*
* Purpose
* =======
*
* PCLATTRS solves one of the triangular systems
*
* A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
*
* with scaling to prevent overflow. Here A is an upper or lower
* triangular matrix, A**T denotes the transpose of A, A**H denotes the
* conjugate transpose of A, x and b are n-element vectors, and s is a
* scaling factor, usually less than or equal to 1, chosen so that the
* components of x will be less than the overflow threshold. If the
* unscaled problem will not cause overflow, the Level 2 PBLAS routine
* PCTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j)
* then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
*
* This is very slow relative to PCTRSV. This should only be used
* when scaling is necessary to control overflow, or when it is modified
* to scale better.
* 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 r x c.
* LOCr( K ) denotes the number of elements of K that a process
* would receive if K were distributed over the r 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 c 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
* =========
*
* UPLO (global input) CHARACTER*1
* Specifies whether the matrix A is upper or lower triangular.
* = 'U': Upper triangular
* = 'L': Lower triangular
*
* TRANS (global input) CHARACTER*1
* Specifies the operation applied to A.
* = 'N': Solve A * x = s*b (No transpose)
* = 'T': Solve A**T * x = s*b (Transpose)
* = 'C': Solve A**H * x = s*b (Conjugate transpose)
*
* DIAG (global input) CHARACTER*1
* Specifies whether or not the matrix A is unit triangular.
* = 'N': Non-unit triangular
* = 'U': Unit triangular
*
* NORMIN (global input) CHARACTER*1
* Specifies whether CNORM has been set or not.
* = 'Y': CNORM contains the column norms on entry
* = 'N': CNORM is not set on entry. On exit, the norms will
* be computed and stored in CNORM.
*
* N (global input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (local input) COMPLEX array, dimension (DESCA(LLD_),*)
* The triangular matrix A. If UPLO = 'U', the leading n by n
* upper triangular part of the array A contains the upper
* triangular matrix, and the strictly lower triangular part of
* A is not referenced. If UPLO = 'L', the leading n by n lower
* triangular part of the array A contains the lower triangular
* matrix, and the strictly upper triangular part of A is not
* referenced. If DIAG = 'U', the diagonal elements of A are
* also not referenced and are assumed to be 1.
*
* IA (global input) pointer to INTEGER
* The global row index of the submatrix of the distributed
* matrix A to operate on.
*
* JA (global input) pointer to INTEGER
* The global column index of the submatrix of the distributed
* matrix A to operate on.
*
* DESCA (global and local input) INTEGER array of dimension DLEN_.
* The array descriptor for the distributed matrix A.
*
* X (local input/output) COMPLEX array,
* dimension (DESCX(LLD_),*)
* On entry, the right hand side b of the triangular system.
* On exit, X is overwritten by the solution vector x.
*
* IX (global input) pointer to INTEGER
* The global row index of the submatrix of the distributed
* matrix X to operate on.
*
* JX (global input) pointer to INTEGER
* The global column index of the submatrix of the distributed
* matrix X to operate on.
*
* DESCX (global and local input) INTEGER array of dimension DLEN_.
* The array descriptor for the distributed matrix X.
*
* SCALE (global output) REAL
* The scaling factor s for the triangular system
* A * x = s*b, A**T * x = s*b, or A**H * x = s*b.
* If SCALE = 0, the matrix A is singular or badly scaled, and
* the vector x is an exact or approximate solution to A*x = 0.
*
* CNORM (global input or global output) REAL array,
* dimension (N)
* If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
* contains the norm of the off-diagonal part of the j-th column
* of A. If TRANS = 'N', CNORM(j) must be greater than or equal
* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
* must be greater than or equal to the 1-norm.
*
* If NORMIN = 'N', CNORM is an output argument and CNORM(j)
* returns the 1-norm of the offdiagonal part of the j-th column
* of A.
*
* INFO (global output) INTEGER
* = 0: successful exit
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* Further Details
* ======= =======
*
* A rough bound on x is computed; if that is less than overflow, PCTRSV
* is called, otherwise, specific code is used which checks for possible
* overflow or divide-by-zero at every operation.
*
* A columnwise scheme is used for solving A*x = b. The basic algorithm
* if A is lower triangular is
*
* x[1:n] := b[1:n]
* for j = 1, ..., n
* x(j) := x(j) / A(j,j)
* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
* end
*
* Define bounds on the components of x after j iterations of the loop:
* M(j) = bound on x[1:j]
* G(j) = bound on x[j+1:n]
* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
*
* Then for iteration j+1 we have
* M(j+1) <= G(j) / | A(j+1,j+1) |
* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
*
* where CNORM(j+1) is greater than or equal to the infinity-norm of
* column j+1 of A, not counting the diagonal. Hence
*
* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
* 1<=i<=j
* and
*
* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
* 1<=i< j
*
* Since |x(j)| <= M(j), we use the Level 2 PBLAS routine PCTRSV if the
* reciprocal of the largest M(j), j=1,..,n, is larger than
* max(underflow, 1/overflow).
*
* The bound on x(j) is also used to determine when a step in the
* columnwise method can be performed without fear of overflow. If
* the computed bound is greater than a large constant, x is scaled to
* prevent overflow, but if the bound overflows, x is set to 0, x(j) to
* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
*
* Similarly, a row-wise scheme is used to solve A**T *x = b or
* A**H *x = b. The basic algorithm for A upper triangular is
*
* for j = 1, ..., n
* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
* end
*
* We simultaneously compute two bounds
* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
* M(j) = bound on x(i), 1<=i<=j
*
* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
* Then the bound on x(j) is
*
* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
*
* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
* 1<=i<=j
*
* and we can safely call PCTRSV if 1/M(n) and 1/G(n) are both greater
* than max(underflow, 1/overflow).
*
* Last modified by: Mark R. Fahey, August 2000
*
* =====================================================================
*
* .. Parameters ..
REAL ZERO, HALF, ONE, TWO
PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0,
$ TWO = 2.0E+0 )
COMPLEX CZERO, CONE
PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
$ CONE = ( 1.0E+0, 0.0E+0 ) )
INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
$ MB_, NB_, RSRC_, CSRC_, LLD_
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 ..
LOGICAL NOTRAN, NOUNIT, UPPER
INTEGER CONTXT, CSRC, I, ICOL, ICOLX, IMAX, IROW,
$ IROWX, ITMP1, ITMP1X, ITMP2, ITMP2X, J, JFIRST,
$ JINC, JLAST, LDA, LDX, MB, MYCOL, MYROW, NB,
$ NPCOL, NPROW, RSRC
REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
$ XBND, XJ, XMAX
COMPLEX CSUMJ, TJJS, USCAL, XJTMP, ZDUM
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER ISAMAX
REAL PSLAMCH
COMPLEX CLADIV
EXTERNAL LSAME, ISAMAX, PSLAMCH, CLADIV
* ..
* .. External Subroutines ..
EXTERNAL BLACS_GRIDINFO, SGSUM2D, SSCAL, INFOG2L,
$ PSCASUM, PSLABAD, PXERBLA, PCAMAX, PCAXPY,
$ PCDOTC, PCDOTU, PCSSCAL, PCLASET, PCSCAL,
$ PCTRSV, CGEBR2D, CGEBS2D
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, REAL, CMPLX, CONJG, AIMAG, MAX, MIN
* ..
* .. Statement Functions ..
REAL CABS1, CABS2
* ..
* .. Statement Function definitions ..
CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
CABS2( ZDUM ) = ABS( REAL( ZDUM ) / 2.E0 ) +
$ ABS( AIMAG( ZDUM ) / 2.E0 )
* ..
* .. Executable Statements ..
*
INFO = 0
UPPER = LSAME( UPLO, 'U' )
NOTRAN = LSAME( TRANS, 'N' )
NOUNIT = LSAME( DIAG, 'N' )
*
CONTXT = DESCA( CTXT_ )
RSRC = DESCA( RSRC_ )
CSRC = DESCA( CSRC_ )
MB = DESCA( MB_ )
NB = DESCA( NB_ )
LDA = DESCA( LLD_ )
LDX = DESCX( LLD_ )
*
* Test the input parameters.
*
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
$ LSAME( TRANS, 'C' ) ) THEN
INFO = -2
ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
INFO = -3
ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
$ LSAME( NORMIN, 'N' ) ) THEN
INFO = -4
ELSE IF( N.LT.0 ) THEN
INFO = -5
END IF
*
CALL BLACS_GRIDINFO( CONTXT, NPROW, NPCOL, MYROW, MYCOL )
*
IF( INFO.NE.0 ) THEN
CALL PXERBLA( CONTXT, 'PCLATTRS', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
*
* Determine machine dependent parameters to control overflow.
*
SMLNUM = PSLAMCH( CONTXT, 'Safe minimum' )
BIGNUM = ONE / SMLNUM
CALL PSLABAD( CONTXT, SMLNUM, BIGNUM )
SMLNUM = SMLNUM / PSLAMCH( CONTXT, 'Precision' )
BIGNUM = ONE / SMLNUM
SCALE = ONE
*
*
IF( LSAME( NORMIN, 'N' ) ) THEN
*
* Compute the 1-norm of each column, not including the diagonal.
*
IF( UPPER ) THEN
*
* A is upper triangular.
*
CNORM( 1 ) = ZERO
DO 10 J = 2, N
CALL PSCASUM( J-1, CNORM( J ), A, IA, JA+J-1, DESCA, 1 )
10 CONTINUE
ELSE
*
* A is lower triangular.
*
DO 20 J = 1, N - 1
CALL PSCASUM( N-J, CNORM( J ), A, IA+J, JA+J-1, DESCA,
$ 1 )
20 CONTINUE
CNORM( N ) = ZERO
END IF
CALL SGSUM2D( CONTXT, 'Row', ' ', N, 1, CNORM, 1, -1, -1 )
END IF
*
* Scale the column norms by TSCAL if the maximum element in CNORM is
* greater than BIGNUM/2.
*
IMAX = ISAMAX( N, CNORM, 1 )
TMAX = CNORM( IMAX )
IF( TMAX.LE.BIGNUM*HALF ) THEN
TSCAL = ONE
ELSE
TSCAL = HALF / ( SMLNUM*TMAX )
CALL SSCAL( N, TSCAL, CNORM, 1 )
END IF
*
* Compute a bound on the computed solution vector to see if the
* Level 2 PBLAS routine PCTRSV can be used.
*
XMAX = ZERO
CALL PCAMAX( N, ZDUM, IMAX, X, IX, JX, DESCX, 1 )
XMAX = CABS2( ZDUM )
CALL SGSUM2D( CONTXT, 'Row', ' ', 1, 1, XMAX, 1, -1, -1 )
XBND = XMAX
*
IF( NOTRAN ) THEN
*
* Compute the growth in A * x = b.
*
IF( UPPER ) THEN
JFIRST = N
JLAST = 1
JINC = -1
ELSE
JFIRST = 1
JLAST = N
JINC = 1
END IF
*
IF( TSCAL.NE.ONE ) THEN
GROW = ZERO
GO TO 50
END IF
*
IF( NOUNIT ) THEN
*
* A is non-unit triangular.
*
* Compute GROW = 1/G(j) and XBND = 1/M(j).
* Initially, G(0) = max{x(i), i=1,...,n}.
*
GROW = HALF / MAX( XBND, SMLNUM )
XBND = GROW
DO 30 J = JFIRST, JLAST, JINC
*
* Exit the loop if the growth factor is too small.
*
IF( GROW.LE.SMLNUM )
$ GO TO 50
*
* TJJS = A( J, J )
CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, MYROW,
$ MYCOL, IROW, ICOL, ITMP1, ITMP2 )
IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
TJJS = A( ( ICOL-1 )*LDA+IROW )
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1,
$ ITMP1, ITMP2 )
END IF
TJJ = CABS1( TJJS )
*
IF( TJJ.GE.SMLNUM ) THEN
*
* M(j) = G(j-1) / abs(A(j,j))
*
XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
ELSE
*
* M(j) could overflow, set XBND to 0.
*
XBND = ZERO
END IF
*
IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
*
* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
*
GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
ELSE
*
* G(j) could overflow, set GROW to 0.
*
GROW = ZERO
END IF
30 CONTINUE
GROW = XBND
ELSE
*
* A is unit triangular.
*
* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
*
GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
DO 40 J = JFIRST, JLAST, JINC
*
* Exit the loop if the growth factor is too small.
*
IF( GROW.LE.SMLNUM )
$ GO TO 50
*
* G(j) = G(j-1)*( 1 + CNORM(j) )
*
GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
40 CONTINUE
END IF
50 CONTINUE
*
ELSE
*
* Compute the growth in A**T * x = b or A**H * x = b.
*
IF( UPPER ) THEN
JFIRST = 1
JLAST = N
JINC = 1
ELSE
JFIRST = N
JLAST = 1
JINC = -1
END IF
*
IF( TSCAL.NE.ONE ) THEN
GROW = ZERO
GO TO 80
END IF
*
IF( NOUNIT ) THEN
*
* A is non-unit triangular.
*
* Compute GROW = 1/G(j) and XBND = 1/M(j).
* Initially, M(0) = max{x(i), i=1,...,n}.
*
GROW = HALF / MAX( XBND, SMLNUM )
XBND = GROW
DO 60 J = JFIRST, JLAST, JINC
*
* Exit the loop if the growth factor is too small.
*
IF( GROW.LE.SMLNUM )
$ GO TO 80
*
* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
*
XJ = ONE + CNORM( J )
GROW = MIN( GROW, XBND / XJ )
*
* TJJS = A( J, J )
CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, MYROW,
$ MYCOL, IROW, ICOL, ITMP1, ITMP2 )
IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
TJJS = A( ( ICOL-1 )*LDA+IROW )
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1,
$ ITMP1, ITMP2 )
END IF
TJJ = CABS1( TJJS )
*
IF( TJJ.GE.SMLNUM ) THEN
*
* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
*
IF( XJ.GT.TJJ )
$ XBND = XBND*( TJJ / XJ )
ELSE
*
* M(j) could overflow, set XBND to 0.
*
XBND = ZERO
END IF
60 CONTINUE
GROW = MIN( GROW, XBND )
ELSE
*
* A is unit triangular.
*
* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
*
GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
DO 70 J = JFIRST, JLAST, JINC
*
* Exit the loop if the growth factor is too small.
*
IF( GROW.LE.SMLNUM )
$ GO TO 80
*
* G(j) = ( 1 + CNORM(j) )*G(j-1)
*
XJ = ONE + CNORM( J )
GROW = GROW / XJ
70 CONTINUE
END IF
80 CONTINUE
END IF
*
IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
*
* Use the Level 2 PBLAS solve if the reciprocal of the bound on
* elements of X is not too small.
*
CALL PCTRSV( UPLO, TRANS, DIAG, N, A, IA, JA, DESCA, X, IX, JX,
$ DESCX, 1 )
ELSE
*
* Use a Level 1 PBLAS solve, scaling intermediate results.
*
IF( XMAX.GT.BIGNUM*HALF ) THEN
*
* Scale X so that its components are less than or equal to
* BIGNUM in absolute value.
*
SCALE = ( BIGNUM*HALF ) / XMAX
CALL PCSSCAL( N, SCALE, X, IX, JX, DESCX, 1 )
XMAX = BIGNUM
ELSE
XMAX = XMAX*TWO
END IF
*
IF( NOTRAN ) THEN
*
* Solve A * x = b
*
DO 100 J = JFIRST, JLAST, JINC
*
* Compute x(j) = b(j) / A(j,j), scaling x if necessary.
*
* XJ = CABS1( X( J ) )
CALL INFOG2L( IX+J-1, JX, DESCX, NPROW, NPCOL, MYROW,
$ MYCOL, IROWX, ICOLX, ITMP1X, ITMP2X )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) THEN
XJTMP = X( IROWX )
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1,
$ ITMP1X, ITMP2X )
END IF
XJ = CABS1( XJTMP )
IF( NOUNIT ) THEN
* TJJS = A( J, J )*TSCAL
CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
$ MYROW, MYCOL, IROW, ICOL, ITMP1, ITMP2 )
IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1,
$ ITMP1, ITMP2 )
END IF
ELSE
TJJS = TSCAL
IF( TSCAL.EQ.ONE )
$ GO TO 90
END IF
TJJ = CABS1( TJJS )
IF( TJJ.GT.SMLNUM ) THEN
*
* abs(A(j,j)) > SMLNUM:
*
IF( TJJ.LT.ONE ) THEN
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale x by 1/b(j).
*
REC = ONE / XJ
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
END IF
* X( J ) = CLADIV( X( J ), TJJS )
* XJ = CABS1( X( J ) )
XJTMP = CLADIV( XJTMP, TJJS )
XJ = CABS1( XJTMP )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ THEN
X( IROWX ) = XJTMP
END IF
ELSE IF( TJJ.GT.ZERO ) THEN
*
* 0 < abs(A(j,j)) <= SMLNUM:
*
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
* to avoid overflow when dividing by A(j,j).
*
REC = ( TJJ*BIGNUM ) / XJ
IF( CNORM( J ).GT.ONE ) THEN
*
* Scale by 1/CNORM(j) to avoid overflow when
* multiplying x(j) times column j.
*
REC = REC / CNORM( J )
END IF
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
* X( J ) = CLADIV( X( J ), TJJS )
* XJ = CABS1( X( J ) )
XJTMP = CLADIV( XJTMP, TJJS )
XJ = CABS1( XJTMP )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ THEN
X( IROWX ) = XJTMP
END IF
ELSE
*
* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
* scale = 0, and compute a solution to A*x = 0.
*
CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX,
$ DESCX )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ THEN
X( IROWX ) = CONE
END IF
XJTMP = CONE
XJ = ONE
SCALE = ZERO
XMAX = ZERO
END IF
90 CONTINUE
*
* Scale x if necessary to avoid overflow when adding a
* multiple of column j of A.
*
IF( XJ.GT.ONE ) THEN
REC = ONE / XJ
IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
*
* Scale x by 1/(2*abs(x(j))).
*
REC = REC*HALF
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
END IF
ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
*
* Scale x by 1/2.
*
CALL PCSSCAL( N, HALF, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*HALF
SCALE = SCALE*HALF
END IF
*
IF( UPPER ) THEN
IF( J.GT.1 ) THEN
*
* Compute the update
* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
*
ZDUM = -XJTMP*TSCAL
CALL PCAXPY( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1, X,
$ IX, JX, DESCX, 1 )
CALL PCAMAX( J-1, ZDUM, IMAX, X, IX, JX, DESCX, 1 )
XMAX = CABS1( ZDUM )
CALL SGSUM2D( CONTXT, 'Row', ' ', 1, 1, XMAX, 1,
$ -1, -1 )
END IF
ELSE
IF( J.LT.N ) THEN
*
* Compute the update
* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
*
ZDUM = -XJTMP*TSCAL
CALL PCAXPY( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1,
$ X, IX+J, JX, DESCX, 1 )
CALL PCAMAX( N-J, ZDUM, I, X, IX+J, JX, DESCX, 1 )
XMAX = CABS1( ZDUM )
CALL SGSUM2D( CONTXT, 'Row', ' ', 1, 1, XMAX, 1,
$ -1, -1 )
END IF
END IF
100 CONTINUE
*
ELSE IF( LSAME( TRANS, 'T' ) ) THEN
*
* Solve A**T * x = b
*
DO 120 J = JFIRST, JLAST, JINC
*
* Compute x(j) = b(j) - sum A(k,j)*x(k).
* k<>j
*
* XJ = CABS1( X( J ) )
CALL INFOG2L( IX+J-1, JX, DESCX, NPROW, NPCOL, MYROW,
$ MYCOL, IROWX, ICOLX, ITMP1X, ITMP2X )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) THEN
XJTMP = X( IROWX )
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1,
$ ITMP1X, ITMP2X )
END IF
XJ = CABS1( XJTMP )
USCAL = CMPLX( TSCAL )
REC = ONE / MAX( XMAX, ONE )
IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
*
* If x(j) could overflow, scale x by 1/(2*XMAX).
*
REC = REC*HALF
IF( NOUNIT ) THEN
* TJJS = A( J, J )*TSCAL
CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
$ MYROW, MYCOL, IROW, ICOL, ITMP1,
$ ITMP2 )
IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) )
$ THEN
TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS,
$ 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1,
$ ITMP1, ITMP2 )
END IF
ELSE
TJJS = TSCAL
END IF
TJJ = CABS1( TJJS )
IF( TJJ.GT.ONE ) THEN
*
* Divide by A(j,j) when scaling x if A(j,j) > 1.
*
REC = MIN( ONE, REC*TJJ )
USCAL = CLADIV( USCAL, TJJS )
END IF
IF( REC.LT.ONE ) THEN
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
END IF
*
CSUMJ = CZERO
IF( USCAL.EQ.CONE ) THEN
*
* If the scaling needed for A in the dot product is 1,
* call PCDOTU to perform the dot product.
*
IF( UPPER ) THEN
CALL PCDOTU( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1,
$ X, IX, JX, DESCX, 1 )
ELSE IF( J.LT.N ) THEN
CALL PCDOTU( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1,
$ X, IX+J, JX, DESCX, 1 )
END IF
IF( MYCOL.EQ.ITMP2X ) THEN
CALL CGEBS2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1,
$ MYROW, ITMP2X )
END IF
ELSE
*
* Otherwise, scale column of A by USCAL before dot
* product. Below is not the best way to do it.
*
IF( UPPER ) THEN
* DO 130 I = 1, J - 1
* CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
* 130 CONTINUE
ZDUM = CONJG( USCAL )
CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 )
CALL PCDOTU( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1,
$ X, IX, JX, DESCX, 1 )
ZDUM = CLADIV( ZDUM, USCAL )
CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 )
ELSE IF( J.LT.N ) THEN
* DO 140 I = J + 1, N
* CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
* 140 CONTINUE
ZDUM = CONJG( USCAL )
CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 )
CALL PCDOTU( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1,
$ X, IX+J, JX, DESCX, 1 )
ZDUM = CLADIV( ZDUM, USCAL )
CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 )
END IF
IF( MYCOL.EQ.ITMP2X ) THEN
CALL CGEBS2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1,
$ MYROW, ITMP2X )
END IF
END IF
*
IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
*
* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
* was not used to scale the dotproduct.
*
* X( J ) = X( J ) - CSUMJ
* XJ = CABS1( X( J ) )
XJTMP = XJTMP - CSUMJ
XJ = CABS1( XJTMP )
* IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
* $ X( IROWX ) = XJTMP
IF( NOUNIT ) THEN
* TJJS = A( J, J )*TSCAL
CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
$ MYROW, MYCOL, IROW, ICOL, ITMP1,
$ ITMP2 )
IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) )
$ THEN
TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS,
$ 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1,
$ ITMP1, ITMP2 )
END IF
ELSE
TJJS = TSCAL
IF( TSCAL.EQ.ONE )
$ GO TO 110
END IF
*
* Compute x(j) = x(j) / A(j,j), scaling if necessary.
*
TJJ = CABS1( TJJS )
IF( TJJ.GT.SMLNUM ) THEN
*
* abs(A(j,j)) > SMLNUM:
*
IF( TJJ.LT.ONE ) THEN
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale X by 1/abs(x(j)).
*
REC = ONE / XJ
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
END IF
* X( J ) = CLADIV( X( J ), TJJS )
XJTMP = CLADIV( XJTMP, TJJS )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ THEN
X( IROWX ) = XJTMP
END IF
ELSE IF( TJJ.GT.ZERO ) THEN
*
* 0 < abs(A(j,j)) <= SMLNUM:
*
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
*
REC = ( TJJ*BIGNUM ) / XJ
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
* X( J ) = CLADIV( X( J ), TJJS )
XJTMP = CLADIV( XJTMP, TJJS )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ THEN
X( IROWX ) = XJTMP
END IF
ELSE
*
* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
* scale = 0 and compute a solution to A**T *x = 0.
*
CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX,
$ DESCX )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ THEN
X( IROWX ) = CONE
END IF
XJTMP = CONE
SCALE = ZERO
XMAX = ZERO
END IF
110 CONTINUE
ELSE
*
* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
* product has already been divided by 1/A(j,j).
*
* X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
XJTMP = CLADIV( XJTMP, TJJS ) - CSUMJ
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ THEN
X( IROWX ) = XJTMP
END IF
END IF
XMAX = MAX( XMAX, CABS1( XJTMP ) )
120 CONTINUE
*
ELSE
*
* Solve A**H * x = b
*
DO 140 J = JFIRST, JLAST, JINC
*
* Compute x(j) = b(j) - sum A(k,j)*x(k).
* k<>j
*
CALL INFOG2L( IX+J-1, JX, DESCX, NPROW, NPCOL, MYROW,
$ MYCOL, IROWX, ICOLX, ITMP1X, ITMP2X )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) THEN
XJTMP = X( IROWX )
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1,
$ ITMP1X, ITMP2X )
END IF
XJ = CABS1( XJTMP )
USCAL = TSCAL
REC = ONE / MAX( XMAX, ONE )
IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
*
* If x(j) could overflow, scale x by 1/(2*XMAX).
*
REC = REC*HALF
IF( NOUNIT ) THEN
* TJJS = CONJG( A( J, J ) )*TSCAL
CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
$ MYROW, MYCOL, IROW, ICOL, ITMP1,
$ ITMP2 )
IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) )
$ THEN
TJJS = CONJG( A( ( ICOL-1 )*LDA+IROW ) )*TSCAL
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS,
$ 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1,
$ ITMP1, ITMP2 )
END IF
ELSE
TJJS = TSCAL
END IF
TJJ = CABS1( TJJS )
IF( TJJ.GT.ONE ) THEN
*
* Divide by A(j,j) when scaling x if A(j,j) > 1.
*
REC = MIN( ONE, REC*TJJ )
USCAL = CLADIV( USCAL, TJJS )
END IF
IF( REC.LT.ONE ) THEN
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
END IF
*
CSUMJ = CZERO
IF( USCAL.EQ.CONE ) THEN
*
* If the scaling needed for A in the dot product is 1,
* call PCDOTC to perform the dot product.
*
IF( UPPER ) THEN
CALL PCDOTC( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1,
$ X, IX, JX, DESCX, 1 )
ELSE IF( J.LT.N ) THEN
CALL PCDOTC( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1,
$ X, IX+J, JX, DESCX, 1 )
END IF
IF( MYCOL.EQ.ITMP2X ) THEN
CALL CGEBS2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1,
$ MYROW, ITMP2X )
END IF
ELSE
*
* Otherwise, scale column of A by USCAL before dot
* product. Below is not the best way to do it.
*
IF( UPPER ) THEN
* DO 180 I = 1, J - 1
* CSUMJ = CSUMJ + ( CONJG( A( I, J ) )*USCAL )*
* $ X( I )
* 180 CONTINUE
ZDUM = CONJG( USCAL )
CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 )
CALL PCDOTC( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1,
$ X, IX, JX, DESCX, 1 )
ZDUM = CLADIV( CONE, ZDUM )
CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 )
ELSE IF( J.LT.N ) THEN
* DO 190 I = J + 1, N
* CSUMJ = CSUMJ + ( CONJG( A( I, J ) )*USCAL )*
* $ X( I )
* 190 CONTINUE
ZDUM = CONJG( USCAL )
CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 )
CALL PCDOTC( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1,
$ X, IX+J, JX, DESCX, 1 )
ZDUM = CLADIV( CONE, ZDUM )
CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 )
END IF
IF( MYCOL.EQ.ITMP2X ) THEN
CALL CGEBS2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1 )
ELSE
CALL CGEBR2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1,
$ MYROW, ITMP2X )
END IF
END IF
*
IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
*
* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
* was not used to scale the dotproduct.
*
* X( J ) = X( J ) - CSUMJ
* XJ = CABS1( X( J ) )
XJTMP = XJTMP - CSUMJ
XJ = CABS1( XJTMP )
* IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
* $ X( IROWX ) = XJTMP
IF( NOUNIT ) THEN
* TJJS = CONJG( A( J, J ) )*TSCAL
CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
$ MYROW, MYCOL, IROW, ICOL, ITMP1,
$ ITMP2 )
IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) )
$ THEN
TJJS = CONJG( A( ( ICOL-1 )*LDA+IROW ) )*TSCAL
CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS,
$ 1 )
ELSE
CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1,
$ ITMP1, ITMP2 )
END IF
ELSE
TJJS = TSCAL
IF( TSCAL.EQ.ONE )
$ GO TO 130
END IF
*
* Compute x(j) = x(j) / A(j,j), scaling if necessary.
*
TJJ = CABS1( TJJS )
IF( TJJ.GT.SMLNUM ) THEN
*
* abs(A(j,j)) > SMLNUM:
*
IF( TJJ.LT.ONE ) THEN
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale X by 1/abs(x(j)).
*
REC = ONE / XJ
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
END IF
* X( J ) = CLADIV( X( J ), TJJS )
XJTMP = CLADIV( XJTMP, TJJS )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ X( IROWX ) = XJTMP
ELSE IF( TJJ.GT.ZERO ) THEN
*
* 0 < abs(A(j,j)) <= SMLNUM:
*
IF( XJ.GT.TJJ*BIGNUM ) THEN
*
* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
*
REC = ( TJJ*BIGNUM ) / XJ
CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
XJTMP = XJTMP*REC
SCALE = SCALE*REC
XMAX = XMAX*REC
END IF
* X( J ) = CLADIV( X( J ), TJJS )
XJTMP = CLADIV( XJTMP, TJJS )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ X( IROWX ) = XJTMP
ELSE
*
* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
* scale = 0 and compute a solution to A**H *x = 0.
*
CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX,
$ DESCX )
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ X( IROWX ) = CONE
XJTMP = CONE
SCALE = ZERO
XMAX = ZERO
END IF
130 CONTINUE
ELSE
*
* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
* product has already been divided by 1/A(j,j).
*
* X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
XJTMP = CLADIV( XJTMP, TJJS ) - CSUMJ
IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
$ X( IROWX ) = XJTMP
END IF
XMAX = MAX( XMAX, CABS1( XJTMP ) )
140 CONTINUE
END IF
SCALE = SCALE / TSCAL
END IF
*
* Scale the column norms by 1/TSCAL for return.
*
IF( TSCAL.NE.ONE ) THEN
CALL SSCAL( N, ONE / TSCAL, CNORM, 1 )
END IF
*
RETURN
*
* End of PCLATTRS
*
END
|