|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PCLATTRS( UPLO, TRANS, DIAG, NORMIN, N, A, IA, JA, 00002 $ DESCA, X, IX, JX, DESCX, SCALE, CNORM, INFO ) 00003 * 00004 * -- ScaLAPACK routine (version 1.7) -- 00005 * University of Tennessee, Knoxville, Oak Ridge National Laboratory, 00006 * and University of California, Berkeley. 00007 * July 31, 2001 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER DIAG, NORMIN, TRANS, UPLO 00011 INTEGER IA, INFO, IX, JA, JX, N 00012 REAL SCALE 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER DESCA( * ), DESCX( * ) 00016 REAL CNORM( * ) 00017 COMPLEX A( * ), X( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * PCLATTRS solves one of the triangular systems 00024 * 00025 * A * x = s*b, A**T * x = s*b, or A**H * x = s*b, 00026 * 00027 * with scaling to prevent overflow. Here A is an upper or lower 00028 * triangular matrix, A**T denotes the transpose of A, A**H denotes the 00029 * conjugate transpose of A, x and b are n-element vectors, and s is a 00030 * scaling factor, usually less than or equal to 1, chosen so that the 00031 * components of x will be less than the overflow threshold. If the 00032 * unscaled problem will not cause overflow, the Level 2 PBLAS routine 00033 * PCTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j) 00034 * then s is set to 0 and a non-trivial solution to A*x = 0 is returned. 00035 * 00036 * This is very slow relative to PCTRSV. This should only be used 00037 * when scaling is necessary to control overflow, or when it is modified 00038 * to scale better. 00039 * Notes 00040 * 00041 * ===== 00042 * 00043 * Each global data object is described by an associated description 00044 * vector. This vector stores the information required to establish 00045 * the mapping between an object element and its corresponding process 00046 * and memory location. 00047 * 00048 * Let A be a generic term for any 2D block cyclicly distributed array. 00049 * Such a global array has an associated description vector DESCA. 00050 * In the following comments, the character _ should be read as 00051 * "of the global array". 00052 * 00053 * NOTATION STORED IN EXPLANATION 00054 * --------------- -------------- -------------------------------------- 00055 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 00056 * DTYPE_A = 1. 00057 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 00058 * the BLACS process grid A is distribu- 00059 * ted over. The context itself is glo- 00060 * bal, but the handle (the integer 00061 * value) may vary. 00062 * M_A (global) DESCA( M_ ) The number of rows in the global 00063 * array A. 00064 * N_A (global) DESCA( N_ ) The number of columns in the global 00065 * array A. 00066 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 00067 * the rows of the array. 00068 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 00069 * the columns of the array. 00070 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 00071 * row of the array A is distributed. 00072 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the 00073 * first column of the array A is 00074 * distributed. 00075 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 00076 * array. LLD_A >= MAX(1,LOCr(M_A)). 00077 * 00078 * Let K be the number of rows or columns of a distributed matrix, 00079 * and assume that its process grid has dimension r x c. 00080 * LOCr( K ) denotes the number of elements of K that a process 00081 * would receive if K were distributed over the r processes of its 00082 * process column. 00083 * Similarly, LOCc( K ) denotes the number of elements of K that a 00084 * process would receive if K were distributed over the c processes of 00085 * its process row. 00086 * The values of LOCr() and LOCc() may be determined via a call to the 00087 * ScaLAPACK tool function, NUMROC: 00088 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 00089 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 00090 * An upper bound for these quantities may be computed by: 00091 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 00092 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 00093 * 00094 * Arguments 00095 * ========= 00096 * 00097 * UPLO (global input) CHARACTER*1 00098 * Specifies whether the matrix A is upper or lower triangular. 00099 * = 'U': Upper triangular 00100 * = 'L': Lower triangular 00101 * 00102 * TRANS (global input) CHARACTER*1 00103 * Specifies the operation applied to A. 00104 * = 'N': Solve A * x = s*b (No transpose) 00105 * = 'T': Solve A**T * x = s*b (Transpose) 00106 * = 'C': Solve A**H * x = s*b (Conjugate transpose) 00107 * 00108 * DIAG (global input) CHARACTER*1 00109 * Specifies whether or not the matrix A is unit triangular. 00110 * = 'N': Non-unit triangular 00111 * = 'U': Unit triangular 00112 * 00113 * NORMIN (global input) CHARACTER*1 00114 * Specifies whether CNORM has been set or not. 00115 * = 'Y': CNORM contains the column norms on entry 00116 * = 'N': CNORM is not set on entry. On exit, the norms will 00117 * be computed and stored in CNORM. 00118 * 00119 * N (global input) INTEGER 00120 * The order of the matrix A. N >= 0. 00121 * 00122 * A (local input) COMPLEX array, dimension (DESCA(LLD_),*) 00123 * The triangular matrix A. If UPLO = 'U', the leading n by n 00124 * upper triangular part of the array A contains the upper 00125 * triangular matrix, and the strictly lower triangular part of 00126 * A is not referenced. If UPLO = 'L', the leading n by n lower 00127 * triangular part of the array A contains the lower triangular 00128 * matrix, and the strictly upper triangular part of A is not 00129 * referenced. If DIAG = 'U', the diagonal elements of A are 00130 * also not referenced and are assumed to be 1. 00131 * 00132 * IA (global input) pointer to INTEGER 00133 * The global row index of the submatrix of the distributed 00134 * matrix A to operate on. 00135 * 00136 * JA (global input) pointer to INTEGER 00137 * The global column index of the submatrix of the distributed 00138 * matrix A to operate on. 00139 * 00140 * DESCA (global and local input) INTEGER array of dimension DLEN_. 00141 * The array descriptor for the distributed matrix A. 00142 * 00143 * X (local input/output) COMPLEX array, 00144 * dimension (DESCX(LLD_),*) 00145 * On entry, the right hand side b of the triangular system. 00146 * On exit, X is overwritten by the solution vector x. 00147 * 00148 * IX (global input) pointer to INTEGER 00149 * The global row index of the submatrix of the distributed 00150 * matrix X to operate on. 00151 * 00152 * JX (global input) pointer to INTEGER 00153 * The global column index of the submatrix of the distributed 00154 * matrix X to operate on. 00155 * 00156 * DESCX (global and local input) INTEGER array of dimension DLEN_. 00157 * The array descriptor for the distributed matrix X. 00158 * 00159 * SCALE (global output) REAL 00160 * The scaling factor s for the triangular system 00161 * A * x = s*b, A**T * x = s*b, or A**H * x = s*b. 00162 * If SCALE = 0, the matrix A is singular or badly scaled, and 00163 * the vector x is an exact or approximate solution to A*x = 0. 00164 * 00165 * CNORM (global input or global output) REAL array, 00166 * dimension (N) 00167 * If NORMIN = 'Y', CNORM is an input argument and CNORM(j) 00168 * contains the norm of the off-diagonal part of the j-th column 00169 * of A. If TRANS = 'N', CNORM(j) must be greater than or equal 00170 * to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) 00171 * must be greater than or equal to the 1-norm. 00172 * 00173 * If NORMIN = 'N', CNORM is an output argument and CNORM(j) 00174 * returns the 1-norm of the offdiagonal part of the j-th column 00175 * of A. 00176 * 00177 * INFO (global output) INTEGER 00178 * = 0: successful exit 00179 * < 0: if INFO = -k, the k-th argument had an illegal value 00180 * 00181 * Further Details 00182 * ======= ======= 00183 * 00184 * A rough bound on x is computed; if that is less than overflow, PCTRSV 00185 * is called, otherwise, specific code is used which checks for possible 00186 * overflow or divide-by-zero at every operation. 00187 * 00188 * A columnwise scheme is used for solving A*x = b. The basic algorithm 00189 * if A is lower triangular is 00190 * 00191 * x[1:n] := b[1:n] 00192 * for j = 1, ..., n 00193 * x(j) := x(j) / A(j,j) 00194 * x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] 00195 * end 00196 * 00197 * Define bounds on the components of x after j iterations of the loop: 00198 * M(j) = bound on x[1:j] 00199 * G(j) = bound on x[j+1:n] 00200 * Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. 00201 * 00202 * Then for iteration j+1 we have 00203 * M(j+1) <= G(j) / | A(j+1,j+1) | 00204 * G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | 00205 * <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) 00206 * 00207 * where CNORM(j+1) is greater than or equal to the infinity-norm of 00208 * column j+1 of A, not counting the diagonal. Hence 00209 * 00210 * G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) 00211 * 1<=i<=j 00212 * and 00213 * 00214 * |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) 00215 * 1<=i< j 00216 * 00217 * Since |x(j)| <= M(j), we use the Level 2 PBLAS routine PCTRSV if the 00218 * reciprocal of the largest M(j), j=1,..,n, is larger than 00219 * max(underflow, 1/overflow). 00220 * 00221 * The bound on x(j) is also used to determine when a step in the 00222 * columnwise method can be performed without fear of overflow. If 00223 * the computed bound is greater than a large constant, x is scaled to 00224 * prevent overflow, but if the bound overflows, x is set to 0, x(j) to 00225 * 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. 00226 * 00227 * Similarly, a row-wise scheme is used to solve A**T *x = b or 00228 * A**H *x = b. The basic algorithm for A upper triangular is 00229 * 00230 * for j = 1, ..., n 00231 * x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) 00232 * end 00233 * 00234 * We simultaneously compute two bounds 00235 * G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j 00236 * M(j) = bound on x(i), 1<=i<=j 00237 * 00238 * The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we 00239 * add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. 00240 * Then the bound on x(j) is 00241 * 00242 * M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | 00243 * 00244 * <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) 00245 * 1<=i<=j 00246 * 00247 * and we can safely call PCTRSV if 1/M(n) and 1/G(n) are both greater 00248 * than max(underflow, 1/overflow). 00249 * 00250 * Last modified by: Mark R. Fahey, August 2000 00251 * 00252 * ===================================================================== 00253 * 00254 * .. Parameters .. 00255 REAL ZERO, HALF, ONE, TWO 00256 PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0, 00257 $ TWO = 2.0E+0 ) 00258 COMPLEX CZERO, CONE 00259 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00260 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00261 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, 00262 $ MB_, NB_, RSRC_, CSRC_, LLD_ 00263 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00264 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00265 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00266 * .. 00267 * .. Local Scalars .. 00268 LOGICAL NOTRAN, NOUNIT, UPPER 00269 INTEGER CONTXT, CSRC, I, ICOL, ICOLX, IMAX, IROW, 00270 $ IROWX, ITMP1, ITMP1X, ITMP2, ITMP2X, J, JFIRST, 00271 $ JINC, JLAST, LDA, LDX, MB, MYCOL, MYROW, NB, 00272 $ NPCOL, NPROW, RSRC 00273 REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL, 00274 $ XBND, XJ, XMAX 00275 COMPLEX CSUMJ, TJJS, USCAL, XJTMP, ZDUM 00276 * .. 00277 * .. External Functions .. 00278 LOGICAL LSAME 00279 INTEGER ISAMAX 00280 REAL PSLAMCH 00281 COMPLEX CLADIV 00282 EXTERNAL LSAME, ISAMAX, PSLAMCH, CLADIV 00283 * .. 00284 * .. External Subroutines .. 00285 EXTERNAL BLACS_GRIDINFO, SGSUM2D, SSCAL, INFOG2L, 00286 $ PSCASUM, PSLABAD, PXERBLA, PCAMAX, PCAXPY, 00287 $ PCDOTC, PCDOTU, PCSSCAL, PCLASET, PCSCAL, 00288 $ PCTRSV, CGEBR2D, CGEBS2D 00289 * .. 00290 * .. Intrinsic Functions .. 00291 INTRINSIC ABS, REAL, CMPLX, CONJG, AIMAG, MAX, MIN 00292 * .. 00293 * .. Statement Functions .. 00294 REAL CABS1, CABS2 00295 * .. 00296 * .. Statement Function definitions .. 00297 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 00298 CABS2( ZDUM ) = ABS( REAL( ZDUM ) / 2.E0 ) + 00299 $ ABS( AIMAG( ZDUM ) / 2.E0 ) 00300 * .. 00301 * .. Executable Statements .. 00302 * 00303 INFO = 0 00304 UPPER = LSAME( UPLO, 'U' ) 00305 NOTRAN = LSAME( TRANS, 'N' ) 00306 NOUNIT = LSAME( DIAG, 'N' ) 00307 * 00308 CONTXT = DESCA( CTXT_ ) 00309 RSRC = DESCA( RSRC_ ) 00310 CSRC = DESCA( CSRC_ ) 00311 MB = DESCA( MB_ ) 00312 NB = DESCA( NB_ ) 00313 LDA = DESCA( LLD_ ) 00314 LDX = DESCX( LLD_ ) 00315 * 00316 * Test the input parameters. 00317 * 00318 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00319 INFO = -1 00320 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 00321 $ LSAME( TRANS, 'C' ) ) THEN 00322 INFO = -2 00323 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 00324 INFO = -3 00325 ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. 00326 $ LSAME( NORMIN, 'N' ) ) THEN 00327 INFO = -4 00328 ELSE IF( N.LT.0 ) THEN 00329 INFO = -5 00330 END IF 00331 * 00332 CALL BLACS_GRIDINFO( CONTXT, NPROW, NPCOL, MYROW, MYCOL ) 00333 * 00334 IF( INFO.NE.0 ) THEN 00335 CALL PXERBLA( CONTXT, 'PCLATTRS', -INFO ) 00336 RETURN 00337 END IF 00338 * 00339 * Quick return if possible 00340 * 00341 IF( N.EQ.0 ) 00342 $ RETURN 00343 * 00344 * Determine machine dependent parameters to control overflow. 00345 * 00346 SMLNUM = PSLAMCH( CONTXT, 'Safe minimum' ) 00347 BIGNUM = ONE / SMLNUM 00348 CALL PSLABAD( CONTXT, SMLNUM, BIGNUM ) 00349 SMLNUM = SMLNUM / PSLAMCH( CONTXT, 'Precision' ) 00350 BIGNUM = ONE / SMLNUM 00351 SCALE = ONE 00352 * 00353 * 00354 IF( LSAME( NORMIN, 'N' ) ) THEN 00355 * 00356 * Compute the 1-norm of each column, not including the diagonal. 00357 * 00358 IF( UPPER ) THEN 00359 * 00360 * A is upper triangular. 00361 * 00362 CNORM( 1 ) = ZERO 00363 DO 10 J = 2, N 00364 CALL PSCASUM( J-1, CNORM( J ), A, IA, JA+J-1, DESCA, 1 ) 00365 10 CONTINUE 00366 ELSE 00367 * 00368 * A is lower triangular. 00369 * 00370 DO 20 J = 1, N - 1 00371 CALL PSCASUM( N-J, CNORM( J ), A, IA+J, JA+J-1, DESCA, 00372 $ 1 ) 00373 20 CONTINUE 00374 CNORM( N ) = ZERO 00375 END IF 00376 CALL SGSUM2D( CONTXT, 'Row', ' ', N, 1, CNORM, 1, -1, -1 ) 00377 END IF 00378 * 00379 * Scale the column norms by TSCAL if the maximum element in CNORM is 00380 * greater than BIGNUM/2. 00381 * 00382 IMAX = ISAMAX( N, CNORM, 1 ) 00383 TMAX = CNORM( IMAX ) 00384 IF( TMAX.LE.BIGNUM*HALF ) THEN 00385 TSCAL = ONE 00386 ELSE 00387 TSCAL = HALF / ( SMLNUM*TMAX ) 00388 CALL SSCAL( N, TSCAL, CNORM, 1 ) 00389 END IF 00390 * 00391 * Compute a bound on the computed solution vector to see if the 00392 * Level 2 PBLAS routine PCTRSV can be used. 00393 * 00394 XMAX = ZERO 00395 CALL PCAMAX( N, ZDUM, IMAX, X, IX, JX, DESCX, 1 ) 00396 XMAX = CABS2( ZDUM ) 00397 CALL SGSUM2D( CONTXT, 'Row', ' ', 1, 1, XMAX, 1, -1, -1 ) 00398 XBND = XMAX 00399 * 00400 IF( NOTRAN ) THEN 00401 * 00402 * Compute the growth in A * x = b. 00403 * 00404 IF( UPPER ) THEN 00405 JFIRST = N 00406 JLAST = 1 00407 JINC = -1 00408 ELSE 00409 JFIRST = 1 00410 JLAST = N 00411 JINC = 1 00412 END IF 00413 * 00414 IF( TSCAL.NE.ONE ) THEN 00415 GROW = ZERO 00416 GO TO 50 00417 END IF 00418 * 00419 IF( NOUNIT ) THEN 00420 * 00421 * A is non-unit triangular. 00422 * 00423 * Compute GROW = 1/G(j) and XBND = 1/M(j). 00424 * Initially, G(0) = max{x(i), i=1,...,n}. 00425 * 00426 GROW = HALF / MAX( XBND, SMLNUM ) 00427 XBND = GROW 00428 DO 30 J = JFIRST, JLAST, JINC 00429 * 00430 * Exit the loop if the growth factor is too small. 00431 * 00432 IF( GROW.LE.SMLNUM ) 00433 $ GO TO 50 00434 * 00435 * TJJS = A( J, J ) 00436 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, MYROW, 00437 $ MYCOL, IROW, ICOL, ITMP1, ITMP2 ) 00438 IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN 00439 TJJS = A( ( ICOL-1 )*LDA+IROW ) 00440 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1 ) 00441 ELSE 00442 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1, 00443 $ ITMP1, ITMP2 ) 00444 END IF 00445 TJJ = CABS1( TJJS ) 00446 * 00447 IF( TJJ.GE.SMLNUM ) THEN 00448 * 00449 * M(j) = G(j-1) / abs(A(j,j)) 00450 * 00451 XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) 00452 ELSE 00453 * 00454 * M(j) could overflow, set XBND to 0. 00455 * 00456 XBND = ZERO 00457 END IF 00458 * 00459 IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN 00460 * 00461 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) 00462 * 00463 GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) 00464 ELSE 00465 * 00466 * G(j) could overflow, set GROW to 0. 00467 * 00468 GROW = ZERO 00469 END IF 00470 30 CONTINUE 00471 GROW = XBND 00472 ELSE 00473 * 00474 * A is unit triangular. 00475 * 00476 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. 00477 * 00478 GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) ) 00479 DO 40 J = JFIRST, JLAST, JINC 00480 * 00481 * Exit the loop if the growth factor is too small. 00482 * 00483 IF( GROW.LE.SMLNUM ) 00484 $ GO TO 50 00485 * 00486 * G(j) = G(j-1)*( 1 + CNORM(j) ) 00487 * 00488 GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) 00489 40 CONTINUE 00490 END IF 00491 50 CONTINUE 00492 * 00493 ELSE 00494 * 00495 * Compute the growth in A**T * x = b or A**H * x = b. 00496 * 00497 IF( UPPER ) THEN 00498 JFIRST = 1 00499 JLAST = N 00500 JINC = 1 00501 ELSE 00502 JFIRST = N 00503 JLAST = 1 00504 JINC = -1 00505 END IF 00506 * 00507 IF( TSCAL.NE.ONE ) THEN 00508 GROW = ZERO 00509 GO TO 80 00510 END IF 00511 * 00512 IF( NOUNIT ) THEN 00513 * 00514 * A is non-unit triangular. 00515 * 00516 * Compute GROW = 1/G(j) and XBND = 1/M(j). 00517 * Initially, M(0) = max{x(i), i=1,...,n}. 00518 * 00519 GROW = HALF / MAX( XBND, SMLNUM ) 00520 XBND = GROW 00521 DO 60 J = JFIRST, JLAST, JINC 00522 * 00523 * Exit the loop if the growth factor is too small. 00524 * 00525 IF( GROW.LE.SMLNUM ) 00526 $ GO TO 80 00527 * 00528 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) 00529 * 00530 XJ = ONE + CNORM( J ) 00531 GROW = MIN( GROW, XBND / XJ ) 00532 * 00533 * TJJS = A( J, J ) 00534 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, MYROW, 00535 $ MYCOL, IROW, ICOL, ITMP1, ITMP2 ) 00536 IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN 00537 TJJS = A( ( ICOL-1 )*LDA+IROW ) 00538 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1 ) 00539 ELSE 00540 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1, 00541 $ ITMP1, ITMP2 ) 00542 END IF 00543 TJJ = CABS1( TJJS ) 00544 * 00545 IF( TJJ.GE.SMLNUM ) THEN 00546 * 00547 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) 00548 * 00549 IF( XJ.GT.TJJ ) 00550 $ XBND = XBND*( TJJ / XJ ) 00551 ELSE 00552 * 00553 * M(j) could overflow, set XBND to 0. 00554 * 00555 XBND = ZERO 00556 END IF 00557 60 CONTINUE 00558 GROW = MIN( GROW, XBND ) 00559 ELSE 00560 * 00561 * A is unit triangular. 00562 * 00563 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. 00564 * 00565 GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) ) 00566 DO 70 J = JFIRST, JLAST, JINC 00567 * 00568 * Exit the loop if the growth factor is too small. 00569 * 00570 IF( GROW.LE.SMLNUM ) 00571 $ GO TO 80 00572 * 00573 * G(j) = ( 1 + CNORM(j) )*G(j-1) 00574 * 00575 XJ = ONE + CNORM( J ) 00576 GROW = GROW / XJ 00577 70 CONTINUE 00578 END IF 00579 80 CONTINUE 00580 END IF 00581 * 00582 IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN 00583 * 00584 * Use the Level 2 PBLAS solve if the reciprocal of the bound on 00585 * elements of X is not too small. 00586 * 00587 CALL PCTRSV( UPLO, TRANS, DIAG, N, A, IA, JA, DESCA, X, IX, JX, 00588 $ DESCX, 1 ) 00589 ELSE 00590 * 00591 * Use a Level 1 PBLAS solve, scaling intermediate results. 00592 * 00593 IF( XMAX.GT.BIGNUM*HALF ) THEN 00594 * 00595 * Scale X so that its components are less than or equal to 00596 * BIGNUM in absolute value. 00597 * 00598 SCALE = ( BIGNUM*HALF ) / XMAX 00599 CALL PCSSCAL( N, SCALE, X, IX, JX, DESCX, 1 ) 00600 XMAX = BIGNUM 00601 ELSE 00602 XMAX = XMAX*TWO 00603 END IF 00604 * 00605 IF( NOTRAN ) THEN 00606 * 00607 * Solve A * x = b 00608 * 00609 DO 100 J = JFIRST, JLAST, JINC 00610 * 00611 * Compute x(j) = b(j) / A(j,j), scaling x if necessary. 00612 * 00613 * XJ = CABS1( X( J ) ) 00614 CALL INFOG2L( IX+J-1, JX, DESCX, NPROW, NPCOL, MYROW, 00615 $ MYCOL, IROWX, ICOLX, ITMP1X, ITMP2X ) 00616 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) THEN 00617 XJTMP = X( IROWX ) 00618 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1 ) 00619 ELSE 00620 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1, 00621 $ ITMP1X, ITMP2X ) 00622 END IF 00623 XJ = CABS1( XJTMP ) 00624 IF( NOUNIT ) THEN 00625 * TJJS = A( J, J )*TSCAL 00626 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, 00627 $ MYROW, MYCOL, IROW, ICOL, ITMP1, ITMP2 ) 00628 IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN 00629 TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL 00630 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1 ) 00631 ELSE 00632 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1, 00633 $ ITMP1, ITMP2 ) 00634 END IF 00635 ELSE 00636 TJJS = TSCAL 00637 IF( TSCAL.EQ.ONE ) 00638 $ GO TO 90 00639 END IF 00640 TJJ = CABS1( TJJS ) 00641 IF( TJJ.GT.SMLNUM ) THEN 00642 * 00643 * abs(A(j,j)) > SMLNUM: 00644 * 00645 IF( TJJ.LT.ONE ) THEN 00646 IF( XJ.GT.TJJ*BIGNUM ) THEN 00647 * 00648 * Scale x by 1/b(j). 00649 * 00650 REC = ONE / XJ 00651 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 00652 XJTMP = XJTMP*REC 00653 SCALE = SCALE*REC 00654 XMAX = XMAX*REC 00655 END IF 00656 END IF 00657 * X( J ) = CLADIV( X( J ), TJJS ) 00658 * XJ = CABS1( X( J ) ) 00659 XJTMP = CLADIV( XJTMP, TJJS ) 00660 XJ = CABS1( XJTMP ) 00661 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 00662 $ THEN 00663 X( IROWX ) = XJTMP 00664 END IF 00665 ELSE IF( TJJ.GT.ZERO ) THEN 00666 * 00667 * 0 < abs(A(j,j)) <= SMLNUM: 00668 * 00669 IF( XJ.GT.TJJ*BIGNUM ) THEN 00670 * 00671 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM 00672 * to avoid overflow when dividing by A(j,j). 00673 * 00674 REC = ( TJJ*BIGNUM ) / XJ 00675 IF( CNORM( J ).GT.ONE ) THEN 00676 * 00677 * Scale by 1/CNORM(j) to avoid overflow when 00678 * multiplying x(j) times column j. 00679 * 00680 REC = REC / CNORM( J ) 00681 END IF 00682 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 00683 XJTMP = XJTMP*REC 00684 SCALE = SCALE*REC 00685 XMAX = XMAX*REC 00686 END IF 00687 * X( J ) = CLADIV( X( J ), TJJS ) 00688 * XJ = CABS1( X( J ) ) 00689 XJTMP = CLADIV( XJTMP, TJJS ) 00690 XJ = CABS1( XJTMP ) 00691 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 00692 $ THEN 00693 X( IROWX ) = XJTMP 00694 END IF 00695 ELSE 00696 * 00697 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 00698 * scale = 0, and compute a solution to A*x = 0. 00699 * 00700 CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX, 00701 $ DESCX ) 00702 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 00703 $ THEN 00704 X( IROWX ) = CONE 00705 END IF 00706 XJTMP = CONE 00707 XJ = ONE 00708 SCALE = ZERO 00709 XMAX = ZERO 00710 END IF 00711 90 CONTINUE 00712 * 00713 * Scale x if necessary to avoid overflow when adding a 00714 * multiple of column j of A. 00715 * 00716 IF( XJ.GT.ONE ) THEN 00717 REC = ONE / XJ 00718 IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN 00719 * 00720 * Scale x by 1/(2*abs(x(j))). 00721 * 00722 REC = REC*HALF 00723 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 00724 XJTMP = XJTMP*REC 00725 SCALE = SCALE*REC 00726 END IF 00727 ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN 00728 * 00729 * Scale x by 1/2. 00730 * 00731 CALL PCSSCAL( N, HALF, X, IX, JX, DESCX, 1 ) 00732 XJTMP = XJTMP*HALF 00733 SCALE = SCALE*HALF 00734 END IF 00735 * 00736 IF( UPPER ) THEN 00737 IF( J.GT.1 ) THEN 00738 * 00739 * Compute the update 00740 * x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) 00741 * 00742 ZDUM = -XJTMP*TSCAL 00743 CALL PCAXPY( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1, X, 00744 $ IX, JX, DESCX, 1 ) 00745 CALL PCAMAX( J-1, ZDUM, IMAX, X, IX, JX, DESCX, 1 ) 00746 XMAX = CABS1( ZDUM ) 00747 CALL SGSUM2D( CONTXT, 'Row', ' ', 1, 1, XMAX, 1, 00748 $ -1, -1 ) 00749 END IF 00750 ELSE 00751 IF( J.LT.N ) THEN 00752 * 00753 * Compute the update 00754 * x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) 00755 * 00756 ZDUM = -XJTMP*TSCAL 00757 CALL PCAXPY( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1, 00758 $ X, IX+J, JX, DESCX, 1 ) 00759 CALL PCAMAX( N-J, ZDUM, I, X, IX+J, JX, DESCX, 1 ) 00760 XMAX = CABS1( ZDUM ) 00761 CALL SGSUM2D( CONTXT, 'Row', ' ', 1, 1, XMAX, 1, 00762 $ -1, -1 ) 00763 END IF 00764 END IF 00765 100 CONTINUE 00766 * 00767 ELSE IF( LSAME( TRANS, 'T' ) ) THEN 00768 * 00769 * Solve A**T * x = b 00770 * 00771 DO 120 J = JFIRST, JLAST, JINC 00772 * 00773 * Compute x(j) = b(j) - sum A(k,j)*x(k). 00774 * k<>j 00775 * 00776 * XJ = CABS1( X( J ) ) 00777 CALL INFOG2L( IX+J-1, JX, DESCX, NPROW, NPCOL, MYROW, 00778 $ MYCOL, IROWX, ICOLX, ITMP1X, ITMP2X ) 00779 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) THEN 00780 XJTMP = X( IROWX ) 00781 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1 ) 00782 ELSE 00783 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1, 00784 $ ITMP1X, ITMP2X ) 00785 END IF 00786 XJ = CABS1( XJTMP ) 00787 USCAL = CMPLX( TSCAL ) 00788 REC = ONE / MAX( XMAX, ONE ) 00789 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN 00790 * 00791 * If x(j) could overflow, scale x by 1/(2*XMAX). 00792 * 00793 REC = REC*HALF 00794 IF( NOUNIT ) THEN 00795 * TJJS = A( J, J )*TSCAL 00796 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, 00797 $ MYROW, MYCOL, IROW, ICOL, ITMP1, 00798 $ ITMP2 ) 00799 IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) 00800 $ THEN 00801 TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL 00802 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 00803 $ 1 ) 00804 ELSE 00805 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1, 00806 $ ITMP1, ITMP2 ) 00807 END IF 00808 ELSE 00809 TJJS = TSCAL 00810 END IF 00811 TJJ = CABS1( TJJS ) 00812 IF( TJJ.GT.ONE ) THEN 00813 * 00814 * Divide by A(j,j) when scaling x if A(j,j) > 1. 00815 * 00816 REC = MIN( ONE, REC*TJJ ) 00817 USCAL = CLADIV( USCAL, TJJS ) 00818 END IF 00819 IF( REC.LT.ONE ) THEN 00820 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 00821 XJTMP = XJTMP*REC 00822 SCALE = SCALE*REC 00823 XMAX = XMAX*REC 00824 END IF 00825 END IF 00826 * 00827 CSUMJ = CZERO 00828 IF( USCAL.EQ.CONE ) THEN 00829 * 00830 * If the scaling needed for A in the dot product is 1, 00831 * call PCDOTU to perform the dot product. 00832 * 00833 IF( UPPER ) THEN 00834 CALL PCDOTU( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1, 00835 $ X, IX, JX, DESCX, 1 ) 00836 ELSE IF( J.LT.N ) THEN 00837 CALL PCDOTU( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1, 00838 $ X, IX+J, JX, DESCX, 1 ) 00839 END IF 00840 IF( MYCOL.EQ.ITMP2X ) THEN 00841 CALL CGEBS2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1 ) 00842 ELSE 00843 CALL CGEBR2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1, 00844 $ MYROW, ITMP2X ) 00845 END IF 00846 ELSE 00847 * 00848 * Otherwise, scale column of A by USCAL before dot 00849 * product. Below is not the best way to do it. 00850 * 00851 IF( UPPER ) THEN 00852 * DO 130 I = 1, J - 1 00853 * CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I ) 00854 * 130 CONTINUE 00855 ZDUM = CONJG( USCAL ) 00856 CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 ) 00857 CALL PCDOTU( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1, 00858 $ X, IX, JX, DESCX, 1 ) 00859 ZDUM = CLADIV( ZDUM, USCAL ) 00860 CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 ) 00861 ELSE IF( J.LT.N ) THEN 00862 * DO 140 I = J + 1, N 00863 * CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I ) 00864 * 140 CONTINUE 00865 ZDUM = CONJG( USCAL ) 00866 CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 ) 00867 CALL PCDOTU( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1, 00868 $ X, IX+J, JX, DESCX, 1 ) 00869 ZDUM = CLADIV( ZDUM, USCAL ) 00870 CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 ) 00871 END IF 00872 IF( MYCOL.EQ.ITMP2X ) THEN 00873 CALL CGEBS2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1 ) 00874 ELSE 00875 CALL CGEBR2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1, 00876 $ MYROW, ITMP2X ) 00877 END IF 00878 END IF 00879 * 00880 IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN 00881 * 00882 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) 00883 * was not used to scale the dotproduct. 00884 * 00885 * X( J ) = X( J ) - CSUMJ 00886 * XJ = CABS1( X( J ) ) 00887 XJTMP = XJTMP - CSUMJ 00888 XJ = CABS1( XJTMP ) 00889 * IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 00890 * $ X( IROWX ) = XJTMP 00891 IF( NOUNIT ) THEN 00892 * TJJS = A( J, J )*TSCAL 00893 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, 00894 $ MYROW, MYCOL, IROW, ICOL, ITMP1, 00895 $ ITMP2 ) 00896 IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) 00897 $ THEN 00898 TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL 00899 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 00900 $ 1 ) 00901 ELSE 00902 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1, 00903 $ ITMP1, ITMP2 ) 00904 END IF 00905 ELSE 00906 TJJS = TSCAL 00907 IF( TSCAL.EQ.ONE ) 00908 $ GO TO 110 00909 END IF 00910 * 00911 * Compute x(j) = x(j) / A(j,j), scaling if necessary. 00912 * 00913 TJJ = CABS1( TJJS ) 00914 IF( TJJ.GT.SMLNUM ) THEN 00915 * 00916 * abs(A(j,j)) > SMLNUM: 00917 * 00918 IF( TJJ.LT.ONE ) THEN 00919 IF( XJ.GT.TJJ*BIGNUM ) THEN 00920 * 00921 * Scale X by 1/abs(x(j)). 00922 * 00923 REC = ONE / XJ 00924 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 00925 XJTMP = XJTMP*REC 00926 SCALE = SCALE*REC 00927 XMAX = XMAX*REC 00928 END IF 00929 END IF 00930 * X( J ) = CLADIV( X( J ), TJJS ) 00931 XJTMP = CLADIV( XJTMP, TJJS ) 00932 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 00933 $ THEN 00934 X( IROWX ) = XJTMP 00935 END IF 00936 ELSE IF( TJJ.GT.ZERO ) THEN 00937 * 00938 * 0 < abs(A(j,j)) <= SMLNUM: 00939 * 00940 IF( XJ.GT.TJJ*BIGNUM ) THEN 00941 * 00942 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. 00943 * 00944 REC = ( TJJ*BIGNUM ) / XJ 00945 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 00946 XJTMP = XJTMP*REC 00947 SCALE = SCALE*REC 00948 XMAX = XMAX*REC 00949 END IF 00950 * X( J ) = CLADIV( X( J ), TJJS ) 00951 XJTMP = CLADIV( XJTMP, TJJS ) 00952 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 00953 $ THEN 00954 X( IROWX ) = XJTMP 00955 END IF 00956 ELSE 00957 * 00958 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 00959 * scale = 0 and compute a solution to A**T *x = 0. 00960 * 00961 CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX, 00962 $ DESCX ) 00963 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 00964 $ THEN 00965 X( IROWX ) = CONE 00966 END IF 00967 XJTMP = CONE 00968 SCALE = ZERO 00969 XMAX = ZERO 00970 END IF 00971 110 CONTINUE 00972 ELSE 00973 * 00974 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot 00975 * product has already been divided by 1/A(j,j). 00976 * 00977 * X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ 00978 XJTMP = CLADIV( XJTMP, TJJS ) - CSUMJ 00979 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 00980 $ THEN 00981 X( IROWX ) = XJTMP 00982 END IF 00983 END IF 00984 XMAX = MAX( XMAX, CABS1( XJTMP ) ) 00985 120 CONTINUE 00986 * 00987 ELSE 00988 * 00989 * Solve A**H * x = b 00990 * 00991 DO 140 J = JFIRST, JLAST, JINC 00992 * 00993 * Compute x(j) = b(j) - sum A(k,j)*x(k). 00994 * k<>j 00995 * 00996 CALL INFOG2L( IX+J-1, JX, DESCX, NPROW, NPCOL, MYROW, 00997 $ MYCOL, IROWX, ICOLX, ITMP1X, ITMP2X ) 00998 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) THEN 00999 XJTMP = X( IROWX ) 01000 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1 ) 01001 ELSE 01002 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, XJTMP, 1, 01003 $ ITMP1X, ITMP2X ) 01004 END IF 01005 XJ = CABS1( XJTMP ) 01006 USCAL = TSCAL 01007 REC = ONE / MAX( XMAX, ONE ) 01008 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN 01009 * 01010 * If x(j) could overflow, scale x by 1/(2*XMAX). 01011 * 01012 REC = REC*HALF 01013 IF( NOUNIT ) THEN 01014 * TJJS = CONJG( A( J, J ) )*TSCAL 01015 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, 01016 $ MYROW, MYCOL, IROW, ICOL, ITMP1, 01017 $ ITMP2 ) 01018 IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) 01019 $ THEN 01020 TJJS = CONJG( A( ( ICOL-1 )*LDA+IROW ) )*TSCAL 01021 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 01022 $ 1 ) 01023 ELSE 01024 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1, 01025 $ ITMP1, ITMP2 ) 01026 END IF 01027 ELSE 01028 TJJS = TSCAL 01029 END IF 01030 TJJ = CABS1( TJJS ) 01031 IF( TJJ.GT.ONE ) THEN 01032 * 01033 * Divide by A(j,j) when scaling x if A(j,j) > 1. 01034 * 01035 REC = MIN( ONE, REC*TJJ ) 01036 USCAL = CLADIV( USCAL, TJJS ) 01037 END IF 01038 IF( REC.LT.ONE ) THEN 01039 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 01040 XJTMP = XJTMP*REC 01041 SCALE = SCALE*REC 01042 XMAX = XMAX*REC 01043 END IF 01044 END IF 01045 * 01046 CSUMJ = CZERO 01047 IF( USCAL.EQ.CONE ) THEN 01048 * 01049 * If the scaling needed for A in the dot product is 1, 01050 * call PCDOTC to perform the dot product. 01051 * 01052 IF( UPPER ) THEN 01053 CALL PCDOTC( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1, 01054 $ X, IX, JX, DESCX, 1 ) 01055 ELSE IF( J.LT.N ) THEN 01056 CALL PCDOTC( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1, 01057 $ X, IX+J, JX, DESCX, 1 ) 01058 END IF 01059 IF( MYCOL.EQ.ITMP2X ) THEN 01060 CALL CGEBS2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1 ) 01061 ELSE 01062 CALL CGEBR2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1, 01063 $ MYROW, ITMP2X ) 01064 END IF 01065 ELSE 01066 * 01067 * Otherwise, scale column of A by USCAL before dot 01068 * product. Below is not the best way to do it. 01069 * 01070 IF( UPPER ) THEN 01071 * DO 180 I = 1, J - 1 01072 * CSUMJ = CSUMJ + ( CONJG( A( I, J ) )*USCAL )* 01073 * $ X( I ) 01074 * 180 CONTINUE 01075 ZDUM = CONJG( USCAL ) 01076 CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 ) 01077 CALL PCDOTC( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1, 01078 $ X, IX, JX, DESCX, 1 ) 01079 ZDUM = CLADIV( CONE, ZDUM ) 01080 CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 ) 01081 ELSE IF( J.LT.N ) THEN 01082 * DO 190 I = J + 1, N 01083 * CSUMJ = CSUMJ + ( CONJG( A( I, J ) )*USCAL )* 01084 * $ X( I ) 01085 * 190 CONTINUE 01086 ZDUM = CONJG( USCAL ) 01087 CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 ) 01088 CALL PCDOTC( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1, 01089 $ X, IX+J, JX, DESCX, 1 ) 01090 ZDUM = CLADIV( CONE, ZDUM ) 01091 CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 ) 01092 END IF 01093 IF( MYCOL.EQ.ITMP2X ) THEN 01094 CALL CGEBS2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1 ) 01095 ELSE 01096 CALL CGEBR2D( CONTXT, 'Row', ' ', 1, 1, CSUMJ, 1, 01097 $ MYROW, ITMP2X ) 01098 END IF 01099 END IF 01100 * 01101 IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN 01102 * 01103 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) 01104 * was not used to scale the dotproduct. 01105 * 01106 * X( J ) = X( J ) - CSUMJ 01107 * XJ = CABS1( X( J ) ) 01108 XJTMP = XJTMP - CSUMJ 01109 XJ = CABS1( XJTMP ) 01110 * IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 01111 * $ X( IROWX ) = XJTMP 01112 IF( NOUNIT ) THEN 01113 * TJJS = CONJG( A( J, J ) )*TSCAL 01114 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL, 01115 $ MYROW, MYCOL, IROW, ICOL, ITMP1, 01116 $ ITMP2 ) 01117 IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) 01118 $ THEN 01119 TJJS = CONJG( A( ( ICOL-1 )*LDA+IROW ) )*TSCAL 01120 CALL CGEBS2D( CONTXT, 'All', ' ', 1, 1, TJJS, 01121 $ 1 ) 01122 ELSE 01123 CALL CGEBR2D( CONTXT, 'All', ' ', 1, 1, TJJS, 1, 01124 $ ITMP1, ITMP2 ) 01125 END IF 01126 ELSE 01127 TJJS = TSCAL 01128 IF( TSCAL.EQ.ONE ) 01129 $ GO TO 130 01130 END IF 01131 * 01132 * Compute x(j) = x(j) / A(j,j), scaling if necessary. 01133 * 01134 TJJ = CABS1( TJJS ) 01135 IF( TJJ.GT.SMLNUM ) THEN 01136 * 01137 * abs(A(j,j)) > SMLNUM: 01138 * 01139 IF( TJJ.LT.ONE ) THEN 01140 IF( XJ.GT.TJJ*BIGNUM ) THEN 01141 * 01142 * Scale X by 1/abs(x(j)). 01143 * 01144 REC = ONE / XJ 01145 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 01146 XJTMP = XJTMP*REC 01147 SCALE = SCALE*REC 01148 XMAX = XMAX*REC 01149 END IF 01150 END IF 01151 * X( J ) = CLADIV( X( J ), TJJS ) 01152 XJTMP = CLADIV( XJTMP, TJJS ) 01153 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 01154 $ X( IROWX ) = XJTMP 01155 ELSE IF( TJJ.GT.ZERO ) THEN 01156 * 01157 * 0 < abs(A(j,j)) <= SMLNUM: 01158 * 01159 IF( XJ.GT.TJJ*BIGNUM ) THEN 01160 * 01161 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. 01162 * 01163 REC = ( TJJ*BIGNUM ) / XJ 01164 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 ) 01165 XJTMP = XJTMP*REC 01166 SCALE = SCALE*REC 01167 XMAX = XMAX*REC 01168 END IF 01169 * X( J ) = CLADIV( X( J ), TJJS ) 01170 XJTMP = CLADIV( XJTMP, TJJS ) 01171 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 01172 $ X( IROWX ) = XJTMP 01173 ELSE 01174 * 01175 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 01176 * scale = 0 and compute a solution to A**H *x = 0. 01177 * 01178 CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX, 01179 $ DESCX ) 01180 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 01181 $ X( IROWX ) = CONE 01182 XJTMP = CONE 01183 SCALE = ZERO 01184 XMAX = ZERO 01185 END IF 01186 130 CONTINUE 01187 ELSE 01188 * 01189 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot 01190 * product has already been divided by 1/A(j,j). 01191 * 01192 * X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ 01193 XJTMP = CLADIV( XJTMP, TJJS ) - CSUMJ 01194 IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) 01195 $ X( IROWX ) = XJTMP 01196 END IF 01197 XMAX = MAX( XMAX, CABS1( XJTMP ) ) 01198 140 CONTINUE 01199 END IF 01200 SCALE = SCALE / TSCAL 01201 END IF 01202 * 01203 * Scale the column norms by 1/TSCAL for return. 01204 * 01205 IF( TSCAL.NE.ONE ) THEN 01206 CALL SSCAL( N, ONE / TSCAL, CNORM, 1 ) 01207 END IF 01208 * 01209 RETURN 01210 * 01211 * End of PCLATTRS 01212 * 01213 END
