ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzlattrs.f
Go to the documentation of this file.
1  SUBROUTINE pzlattrs( UPLO, TRANS, DIAG, NORMIN, N, A, IA, JA,
2  $ DESCA, X, IX, JX, DESCX, SCALE, CNORM, INFO )
3 *
4 * -- ScaLAPACK routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * July 31, 2001
8 *
9 * .. Scalar Arguments ..
10  CHARACTER DIAG, NORMIN, TRANS, UPLO
11  INTEGER IA, INFO, IX, JA, JX, N
12  DOUBLE PRECISION SCALE
13 * ..
14 * .. Array Arguments ..
15  INTEGER DESCA( * ), DESCX( * )
16  DOUBLE PRECISION CNORM( * )
17  COMPLEX*16 A( * ), X( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PZLATTRS solves one of the triangular systems
24 *
25 * A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
26 *
27 * with scaling to prevent overflow. Here A is an upper or lower
28 * triangular matrix, A**T denotes the transpose of A, A**H denotes the
29 * conjugate transpose of A, x and b are n-element vectors, and s is a
30 * scaling factor, usually less than or equal to 1, chosen so that the
31 * components of x will be less than the overflow threshold. If the
32 * unscaled problem will not cause overflow, the Level 2 PBLAS routine
33 * PZTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j)
34 * then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
35 *
36 * This is very slow relative to PZTRSV. This should only be used
37 * when scaling is necessary to control overflow, or when it is modified
38 * to scale better.
39 * Notes
40 *
41 * =====
42 *
43 * Each global data object is described by an associated description
44 * vector. This vector stores the information required to establish
45 * the mapping between an object element and its corresponding process
46 * and memory location.
47 *
48 * Let A be a generic term for any 2D block cyclicly distributed array.
49 * Such a global array has an associated description vector DESCA.
50 * In the following comments, the character _ should be read as
51 * "of the global array".
52 *
53 * NOTATION STORED IN EXPLANATION
54 * --------------- -------------- --------------------------------------
55 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
56 * DTYPE_A = 1.
57 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
58 * the BLACS process grid A is distribu-
59 * ted over. The context itself is glo-
60 * bal, but the handle (the integer
61 * value) may vary.
62 * M_A (global) DESCA( M_ ) The number of rows in the global
63 * array A.
64 * N_A (global) DESCA( N_ ) The number of columns in the global
65 * array A.
66 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
67 * the rows of the array.
68 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
69 * the columns of the array.
70 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
71 * row of the array A is distributed.
72 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
73 * first column of the array A is
74 * distributed.
75 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
76 * array. LLD_A >= MAX(1,LOCr(M_A)).
77 *
78 * Let K be the number of rows or columns of a distributed matrix,
79 * and assume that its process grid has dimension r x c.
80 * LOCr( K ) denotes the number of elements of K that a process
81 * would receive if K were distributed over the r processes of its
82 * process column.
83 * Similarly, LOCc( K ) denotes the number of elements of K that a
84 * process would receive if K were distributed over the c processes of
85 * its process row.
86 * The values of LOCr() and LOCc() may be determined via a call to the
87 * ScaLAPACK tool function, NUMROC:
88 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
89 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
90 * An upper bound for these quantities may be computed by:
91 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
92 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
93 *
94 * Arguments
95 * =========
96 *
97 * UPLO (global input) CHARACTER*1
98 * Specifies whether the matrix A is upper or lower triangular.
99 * = 'U': Upper triangular
100 * = 'L': Lower triangular
101 *
102 * TRANS (global input) CHARACTER*1
103 * Specifies the operation applied to A.
104 * = 'N': Solve A * x = s*b (No transpose)
105 * = 'T': Solve A**T * x = s*b (Transpose)
106 * = 'C': Solve A**H * x = s*b (Conjugate transpose)
107 *
108 * DIAG (global input) CHARACTER*1
109 * Specifies whether or not the matrix A is unit triangular.
110 * = 'N': Non-unit triangular
111 * = 'U': Unit triangular
112 *
113 * NORMIN (global input) CHARACTER*1
114 * Specifies whether CNORM has been set or not.
115 * = 'Y': CNORM contains the column norms on entry
116 * = 'N': CNORM is not set on entry. On exit, the norms will
117 * be computed and stored in CNORM.
118 *
119 * N (global input) INTEGER
120 * The order of the matrix A. N >= 0.
121 *
122 * A (local input) COMPLEX*16 array, dimension (DESCA(LLD_),*)
123 * The triangular matrix A. If UPLO = 'U', the leading n by n
124 * upper triangular part of the array A contains the upper
125 * triangular matrix, and the strictly lower triangular part of
126 * A is not referenced. If UPLO = 'L', the leading n by n lower
127 * triangular part of the array A contains the lower triangular
128 * matrix, and the strictly upper triangular part of A is not
129 * referenced. If DIAG = 'U', the diagonal elements of A are
130 * also not referenced and are assumed to be 1.
131 *
132 * IA (global input) pointer to INTEGER
133 * The global row index of the submatrix of the distributed
134 * matrix A to operate on.
135 *
136 * JA (global input) pointer to INTEGER
137 * The global column index of the submatrix of the distributed
138 * matrix A to operate on.
139 *
140 * DESCA (global and local input) INTEGER array of dimension DLEN_.
141 * The array descriptor for the distributed matrix A.
142 *
143 * X (local input/output) COMPLEX*16 array,
144 * dimension (DESCX(LLD_),*)
145 * On entry, the right hand side b of the triangular system.
146 * On exit, X is overwritten by the solution vector x.
147 *
148 * IX (global input) pointer to INTEGER
149 * The global row index of the submatrix of the distributed
150 * matrix X to operate on.
151 *
152 * JX (global input) pointer to INTEGER
153 * The global column index of the submatrix of the distributed
154 * matrix X to operate on.
155 *
156 * DESCX (global and local input) INTEGER array of dimension DLEN_.
157 * The array descriptor for the distributed matrix X.
158 *
159 * SCALE (global output) DOUBLE PRECISION
160 * The scaling factor s for the triangular system
161 * A * x = s*b, A**T * x = s*b, or A**H * x = s*b.
162 * If SCALE = 0, the matrix A is singular or badly scaled, and
163 * the vector x is an exact or approximate solution to A*x = 0.
164 *
165 * CNORM (global input or global output) DOUBLE PRECISION array,
166 * dimension (N)
167 * If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
168 * contains the norm of the off-diagonal part of the j-th column
169 * of A. If TRANS = 'N', CNORM(j) must be greater than or equal
170 * to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
171 * must be greater than or equal to the 1-norm.
172 *
173 * If NORMIN = 'N', CNORM is an output argument and CNORM(j)
174 * returns the 1-norm of the offdiagonal part of the j-th column
175 * of A.
176 *
177 * INFO (global output) INTEGER
178 * = 0: successful exit
179 * < 0: if INFO = -k, the k-th argument had an illegal value
180 *
181 * Further Details
182 * ======= =======
183 *
184 * A rough bound on x is computed; if that is less than overflow, PZTRSV
185 * is called, otherwise, specific code is used which checks for possible
186 * overflow or divide-by-zero at every operation.
187 *
188 * A columnwise scheme is used for solving A*x = b. The basic algorithm
189 * if A is lower triangular is
190 *
191 * x[1:n] := b[1:n]
192 * for j = 1, ..., n
193 * x(j) := x(j) / A(j,j)
194 * x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
195 * end
196 *
197 * Define bounds on the components of x after j iterations of the loop:
198 * M(j) = bound on x[1:j]
199 * G(j) = bound on x[j+1:n]
200 * Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
201 *
202 * Then for iteration j+1 we have
203 * M(j+1) <= G(j) / | A(j+1,j+1) |
204 * G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
205 * <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
206 *
207 * where CNORM(j+1) is greater than or equal to the infinity-norm of
208 * column j+1 of A, not counting the diagonal. Hence
209 *
210 * G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
211 * 1<=i<=j
212 * and
213 *
214 * |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
215 * 1<=i< j
216 *
217 * Since |x(j)| <= M(j), we use the Level 2 PBLAS routine PZTRSV if the
218 * reciprocal of the largest M(j), j=1,..,n, is larger than
219 * max(underflow, 1/overflow).
220 *
221 * The bound on x(j) is also used to determine when a step in the
222 * columnwise method can be performed without fear of overflow. If
223 * the computed bound is greater than a large constant, x is scaled to
224 * prevent overflow, but if the bound overflows, x is set to 0, x(j) to
225 * 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
226 *
227 * Similarly, a row-wise scheme is used to solve A**T *x = b or
228 * A**H *x = b. The basic algorithm for A upper triangular is
229 *
230 * for j = 1, ..., n
231 * x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
232 * end
233 *
234 * We simultaneously compute two bounds
235 * G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
236 * M(j) = bound on x(i), 1<=i<=j
237 *
238 * The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
239 * add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
240 * Then the bound on x(j) is
241 *
242 * M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
243 *
244 * <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
245 * 1<=i<=j
246 *
247 * and we can safely call PZTRSV if 1/M(n) and 1/G(n) are both greater
248 * than max(underflow, 1/overflow).
249 *
250 * Last modified by: Mark R. Fahey, August 2000
251 *
252 * =====================================================================
253 *
254 * .. Parameters ..
255  DOUBLE PRECISION ZERO, HALF, ONE, TWO
256  parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0,
257  $ two = 2.0d+0 )
258  COMPLEX*16 CZERO, CONE
259  parameter( czero = ( 0.0d+0, 0.0d+0 ),
260  $ cone = ( 1.0d+0, 0.0d+0 ) )
261  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
262  $ mb_, nb_, rsrc_, csrc_, lld_
263  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
264  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
265  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
266 * ..
267 * .. Local Scalars ..
268  LOGICAL NOTRAN, NOUNIT, UPPER
269  INTEGER CONTXT, CSRC, I, ICOL, ICOLX, IMAX, IROW,
270  $ irowx, itmp1, itmp1x, itmp2, itmp2x, j, jfirst,
271  $ jinc, jlast, lda, ldx, mb, mycol, myrow, nb,
272  $ npcol, nprow, rsrc
273  DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
274  $ xbnd, xj, xmax
275  COMPLEX*16 CSUMJ, TJJS, USCAL, XJTMP, ZDUM
276 * ..
277 * .. External Functions ..
278  LOGICAL LSAME
279  INTEGER IDAMAX
280  DOUBLE PRECISION PDLAMCH
281  COMPLEX*16 ZLADIV
282  EXTERNAL lsame, idamax, pdlamch, zladiv
283 * ..
284 * .. External Subroutines ..
285  EXTERNAL blacs_gridinfo, dgsum2d, dscal, infog2l,
286  $ pdlabad, pdzasum, pxerbla, pzamax, pzaxpy,
287  $ pzdotc, pzdotu, pzdscal, pzlaset, pzscal,
288  $ pztrsv, zgebr2d, zgebs2d
289 * ..
290 * .. Intrinsic Functions ..
291  INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
292 * ..
293 * .. Statement Functions ..
294  DOUBLE PRECISION CABS1, CABS2
295 * ..
296 * .. Statement Function definitions ..
297  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
298  cabs2( zdum ) = abs( dble( zdum ) / 2.d0 ) +
299  $ abs( dimag( zdum ) / 2.d0 )
300 * ..
301 * .. Executable Statements ..
302 *
303  info = 0
304  upper = lsame( uplo, 'U' )
305  notran = lsame( trans, 'N' )
306  nounit = lsame( diag, 'N' )
307 *
308  contxt = desca( ctxt_ )
309  rsrc = desca( rsrc_ )
310  csrc = desca( csrc_ )
311  mb = desca( mb_ )
312  nb = desca( nb_ )
313  lda = desca( lld_ )
314  ldx = descx( lld_ )
315 *
316 * Test the input parameters.
317 *
318  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
319  info = -1
320  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
321  $ lsame( trans, 'C' ) ) THEN
322  info = -2
323  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
324  info = -3
325  ELSE IF( .NOT.lsame( normin, 'Y' ) .AND. .NOT.
326  $ lsame( normin, 'N' ) ) THEN
327  info = -4
328  ELSE IF( n.LT.0 ) THEN
329  info = -5
330  END IF
331 *
332  CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
333 *
334  IF( info.NE.0 ) THEN
335  CALL pxerbla( contxt, 'PZLATTRS', -info )
336  RETURN
337  END IF
338 *
339 * Quick return if possible
340 *
341  IF( n.EQ.0 )
342  $ RETURN
343 *
344 * Determine machine dependent parameters to control overflow.
345 *
346  smlnum = pdlamch( contxt, 'Safe minimum' )
347  bignum = one / smlnum
348  CALL pdlabad( contxt, smlnum, bignum )
349  smlnum = smlnum / pdlamch( contxt, 'Precision' )
350  bignum = one / smlnum
351  scale = one
352 *
353 *
354  IF( lsame( normin, 'N' ) ) THEN
355 *
356 * Compute the 1-norm of each column, not including the diagonal.
357 *
358  IF( upper ) THEN
359 *
360 * A is upper triangular.
361 *
362  cnorm( 1 ) = zero
363  DO 10 j = 2, n
364  CALL pdzasum( j-1, cnorm( j ), a, ia, ja+j-1, desca, 1 )
365  10 CONTINUE
366  ELSE
367 *
368 * A is lower triangular.
369 *
370  DO 20 j = 1, n - 1
371  CALL pdzasum( n-j, cnorm( j ), a, ia+j, ja+j-1, desca,
372  $ 1 )
373  20 CONTINUE
374  cnorm( n ) = zero
375  END IF
376  CALL dgsum2d( contxt, 'Row', ' ', n, 1, cnorm, 1, -1, -1 )
377  END IF
378 *
379 * Scale the column norms by TSCAL if the maximum element in CNORM is
380 * greater than BIGNUM/2.
381 *
382  imax = idamax( n, cnorm, 1 )
383  tmax = cnorm( imax )
384  IF( tmax.LE.bignum*half ) THEN
385  tscal = one
386  ELSE
387  tscal = half / ( smlnum*tmax )
388  CALL dscal( n, tscal, cnorm, 1 )
389  END IF
390 *
391 * Compute a bound on the computed solution vector to see if the
392 * Level 2 PBLAS routine PZTRSV can be used.
393 *
394  xmax = zero
395  CALL pzamax( n, zdum, imax, x, ix, jx, descx, 1 )
396  xmax = cabs2( zdum )
397  CALL dgsum2d( contxt, 'Row', ' ', 1, 1, xmax, 1, -1, -1 )
398  xbnd = xmax
399 *
400  IF( notran ) THEN
401 *
402 * Compute the growth in A * x = b.
403 *
404  IF( upper ) THEN
405  jfirst = n
406  jlast = 1
407  jinc = -1
408  ELSE
409  jfirst = 1
410  jlast = n
411  jinc = 1
412  END IF
413 *
414  IF( tscal.NE.one ) THEN
415  grow = zero
416  GO TO 50
417  END IF
418 *
419  IF( nounit ) THEN
420 *
421 * A is non-unit triangular.
422 *
423 * Compute GROW = 1/G(j) and XBND = 1/M(j).
424 * Initially, G(0) = max{x(i), i=1,...,n}.
425 *
426  grow = half / max( xbnd, smlnum )
427  xbnd = grow
428  DO 30 j = jfirst, jlast, jinc
429 *
430 * Exit the loop if the growth factor is too small.
431 *
432  IF( grow.LE.smlnum )
433  $ GO TO 50
434 *
435 * TJJS = A( J, J )
436  CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol, myrow,
437  $ mycol, irow, icol, itmp1, itmp2 )
438  IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
439  tjjs = a( ( icol-1 )*lda+irow )
440  CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs, 1 )
441  ELSE
442  CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
443  $ itmp1, itmp2 )
444  END IF
445  tjj = cabs1( tjjs )
446 *
447  IF( tjj.GE.smlnum ) THEN
448 *
449 * M(j) = G(j-1) / abs(A(j,j))
450 *
451  xbnd = min( xbnd, min( one, tjj )*grow )
452  ELSE
453 *
454 * M(j) could overflow, set XBND to 0.
455 *
456  xbnd = zero
457  END IF
458 *
459  IF( tjj+cnorm( j ).GE.smlnum ) THEN
460 *
461 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
462 *
463  grow = grow*( tjj / ( tjj+cnorm( j ) ) )
464  ELSE
465 *
466 * G(j) could overflow, set GROW to 0.
467 *
468  grow = zero
469  END IF
470  30 CONTINUE
471  grow = xbnd
472  ELSE
473 *
474 * A is unit triangular.
475 *
476 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
477 *
478  grow = min( one, half / max( xbnd, smlnum ) )
479  DO 40 j = jfirst, jlast, jinc
480 *
481 * Exit the loop if the growth factor is too small.
482 *
483  IF( grow.LE.smlnum )
484  $ GO TO 50
485 *
486 * G(j) = G(j-1)*( 1 + CNORM(j) )
487 *
488  grow = grow*( one / ( one+cnorm( j ) ) )
489  40 CONTINUE
490  END IF
491  50 CONTINUE
492 *
493  ELSE
494 *
495 * Compute the growth in A**T * x = b or A**H * x = b.
496 *
497  IF( upper ) THEN
498  jfirst = 1
499  jlast = n
500  jinc = 1
501  ELSE
502  jfirst = n
503  jlast = 1
504  jinc = -1
505  END IF
506 *
507  IF( tscal.NE.one ) THEN
508  grow = zero
509  GO TO 80
510  END IF
511 *
512  IF( nounit ) THEN
513 *
514 * A is non-unit triangular.
515 *
516 * Compute GROW = 1/G(j) and XBND = 1/M(j).
517 * Initially, M(0) = max{x(i), i=1,...,n}.
518 *
519  grow = half / max( xbnd, smlnum )
520  xbnd = grow
521  DO 60 j = jfirst, jlast, jinc
522 *
523 * Exit the loop if the growth factor is too small.
524 *
525  IF( grow.LE.smlnum )
526  $ GO TO 80
527 *
528 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
529 *
530  xj = one + cnorm( j )
531  grow = min( grow, xbnd / xj )
532 *
533 * TJJS = A( J, J )
534  CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol, myrow,
535  $ mycol, irow, icol, itmp1, itmp2 )
536  IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
537  tjjs = a( ( icol-1 )*lda+irow )
538  CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs, 1 )
539  ELSE
540  CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
541  $ itmp1, itmp2 )
542  END IF
543  tjj = cabs1( tjjs )
544 *
545  IF( tjj.GE.smlnum ) THEN
546 *
547 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
548 *
549  IF( xj.GT.tjj )
550  $ xbnd = xbnd*( tjj / xj )
551  ELSE
552 *
553 * M(j) could overflow, set XBND to 0.
554 *
555  xbnd = zero
556  END IF
557  60 CONTINUE
558  grow = min( grow, xbnd )
559  ELSE
560 *
561 * A is unit triangular.
562 *
563 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
564 *
565  grow = min( one, half / max( xbnd, smlnum ) )
566  DO 70 j = jfirst, jlast, jinc
567 *
568 * Exit the loop if the growth factor is too small.
569 *
570  IF( grow.LE.smlnum )
571  $ GO TO 80
572 *
573 * G(j) = ( 1 + CNORM(j) )*G(j-1)
574 *
575  xj = one + cnorm( j )
576  grow = grow / xj
577  70 CONTINUE
578  END IF
579  80 CONTINUE
580  END IF
581 *
582  IF( ( grow*tscal ).GT.smlnum ) THEN
583 *
584 * Use the Level 2 PBLAS solve if the reciprocal of the bound on
585 * elements of X is not too small.
586 *
587  CALL pztrsv( uplo, trans, diag, n, a, ia, ja, desca, x, ix, jx,
588  $ descx, 1 )
589  ELSE
590 *
591 * Use a Level 1 PBLAS solve, scaling intermediate results.
592 *
593  IF( xmax.GT.bignum*half ) THEN
594 *
595 * Scale X so that its components are less than or equal to
596 * BIGNUM in absolute value.
597 *
598  scale = ( bignum*half ) / xmax
599  CALL pzdscal( n, scale, x, ix, jx, descx, 1 )
600  xmax = bignum
601  ELSE
602  xmax = xmax*two
603  END IF
604 *
605  IF( notran ) THEN
606 *
607 * Solve A * x = b
608 *
609  DO 100 j = jfirst, jlast, jinc
610 *
611 * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
612 *
613 * XJ = CABS1( X( J ) )
614  CALL infog2l( ix+j-1, jx, descx, nprow, npcol, myrow,
615  $ mycol, irowx, icolx, itmp1x, itmp2x )
616  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) ) THEN
617  xjtmp = x( irowx )
618  CALL zgebs2d( contxt, 'All', ' ', 1, 1, xjtmp, 1 )
619  ELSE
620  CALL zgebr2d( contxt, 'All', ' ', 1, 1, xjtmp, 1,
621  $ itmp1x, itmp2x )
622  END IF
623  xj = cabs1( xjtmp )
624  IF( nounit ) THEN
625 * TJJS = A( J, J )*TSCAL
626  CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
627  $ myrow, mycol, irow, icol, itmp1, itmp2 )
628  IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
629  tjjs = a( ( icol-1 )*lda+irow )*tscal
630  CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs, 1 )
631  ELSE
632  CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
633  $ itmp1, itmp2 )
634  END IF
635  ELSE
636  tjjs = tscal
637  IF( tscal.EQ.one )
638  $ GO TO 90
639  END IF
640  tjj = cabs1( tjjs )
641  IF( tjj.GT.smlnum ) THEN
642 *
643 * abs(A(j,j)) > SMLNUM:
644 *
645  IF( tjj.LT.one ) THEN
646  IF( xj.GT.tjj*bignum ) THEN
647 *
648 * Scale x by 1/b(j).
649 *
650  rec = one / xj
651  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
652  xjtmp = xjtmp*rec
653  scale = scale*rec
654  xmax = xmax*rec
655  END IF
656  END IF
657 * X( J ) = ZLADIV( X( J ), TJJS )
658 * XJ = CABS1( X( J ) )
659  xjtmp = zladiv( xjtmp, tjjs )
660  xj = cabs1( xjtmp )
661  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
662  $ THEN
663  x( irowx ) = xjtmp
664  END IF
665  ELSE IF( tjj.GT.zero ) THEN
666 *
667 * 0 < abs(A(j,j)) <= SMLNUM:
668 *
669  IF( xj.GT.tjj*bignum ) THEN
670 *
671 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
672 * to avoid overflow when dividing by A(j,j).
673 *
674  rec = ( tjj*bignum ) / xj
675  IF( cnorm( j ).GT.one ) THEN
676 *
677 * Scale by 1/CNORM(j) to avoid overflow when
678 * multiplying x(j) times column j.
679 *
680  rec = rec / cnorm( j )
681  END IF
682  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
683  xjtmp = xjtmp*rec
684  scale = scale*rec
685  xmax = xmax*rec
686  END IF
687 * X( J ) = ZLADIV( X( J ), TJJS )
688 * XJ = CABS1( X( J ) )
689  xjtmp = zladiv( xjtmp, tjjs )
690  xj = cabs1( xjtmp )
691  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
692  $ THEN
693  x( irowx ) = xjtmp
694  END IF
695  ELSE
696 *
697 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
698 * scale = 0, and compute a solution to A*x = 0.
699 *
700  CALL pzlaset( ' ', n, 1, czero, czero, x, ix, jx,
701  $ descx )
702  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
703  $ THEN
704  x( irowx ) = cone
705  END IF
706  xjtmp = cone
707  xj = one
708  scale = zero
709  xmax = zero
710  END IF
711  90 CONTINUE
712 *
713 * Scale x if necessary to avoid overflow when adding a
714 * multiple of column j of A.
715 *
716  IF( xj.GT.one ) THEN
717  rec = one / xj
718  IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
719 *
720 * Scale x by 1/(2*abs(x(j))).
721 *
722  rec = rec*half
723  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
724  xjtmp = xjtmp*rec
725  scale = scale*rec
726  END IF
727  ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
728 *
729 * Scale x by 1/2.
730 *
731  CALL pzdscal( n, half, x, ix, jx, descx, 1 )
732  xjtmp = xjtmp*half
733  scale = scale*half
734  END IF
735 *
736  IF( upper ) THEN
737  IF( j.GT.1 ) THEN
738 *
739 * Compute the update
740 * x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
741 *
742  zdum = -xjtmp*tscal
743  CALL pzaxpy( j-1, zdum, a, ia, ja+j-1, desca, 1, x,
744  $ ix, jx, descx, 1 )
745  CALL pzamax( j-1, zdum, imax, x, ix, jx, descx, 1 )
746  xmax = cabs1( zdum )
747  CALL dgsum2d( contxt, 'Row', ' ', 1, 1, xmax, 1,
748  $ -1, -1 )
749  END IF
750  ELSE
751  IF( j.LT.n ) THEN
752 *
753 * Compute the update
754 * x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
755 *
756  zdum = -xjtmp*tscal
757  CALL pzaxpy( n-j, zdum, a, ia+j, ja+j-1, desca, 1,
758  $ x, ix+j, jx, descx, 1 )
759  CALL pzamax( n-j, zdum, i, x, ix+j, jx, descx, 1 )
760  xmax = cabs1( zdum )
761  CALL dgsum2d( contxt, 'Row', ' ', 1, 1, xmax, 1,
762  $ -1, -1 )
763  END IF
764  END IF
765  100 CONTINUE
766 *
767  ELSE IF( lsame( trans, 'T' ) ) THEN
768 *
769 * Solve A**T * x = b
770 *
771  DO 120 j = jfirst, jlast, jinc
772 *
773 * Compute x(j) = b(j) - sum A(k,j)*x(k).
774 * k<>j
775 *
776 * XJ = CABS1( X( J ) )
777  CALL infog2l( ix+j-1, jx, descx, nprow, npcol, myrow,
778  $ mycol, irowx, icolx, itmp1x, itmp2x )
779  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) ) THEN
780  xjtmp = x( irowx )
781  CALL zgebs2d( contxt, 'All', ' ', 1, 1, xjtmp, 1 )
782  ELSE
783  CALL zgebr2d( contxt, 'All', ' ', 1, 1, xjtmp, 1,
784  $ itmp1x, itmp2x )
785  END IF
786  xj = cabs1( xjtmp )
787  uscal = dcmplx( tscal )
788  rec = one / max( xmax, one )
789  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
790 *
791 * If x(j) could overflow, scale x by 1/(2*XMAX).
792 *
793  rec = rec*half
794  IF( nounit ) THEN
795 * TJJS = A( J, J )*TSCAL
796  CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
797  $ myrow, mycol, irow, icol, itmp1,
798  $ itmp2 )
799  IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
800  $ THEN
801  tjjs = a( ( icol-1 )*lda+irow )*tscal
802  CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs,
803  $ 1 )
804  ELSE
805  CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
806  $ itmp1, itmp2 )
807  END IF
808  ELSE
809  tjjs = tscal
810  END IF
811  tjj = cabs1( tjjs )
812  IF( tjj.GT.one ) THEN
813 *
814 * Divide by A(j,j) when scaling x if A(j,j) > 1.
815 *
816  rec = min( one, rec*tjj )
817  uscal = zladiv( uscal, tjjs )
818  END IF
819  IF( rec.LT.one ) THEN
820  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
821  xjtmp = xjtmp*rec
822  scale = scale*rec
823  xmax = xmax*rec
824  END IF
825  END IF
826 *
827  csumj = czero
828  IF( uscal.EQ.cone ) THEN
829 *
830 * If the scaling needed for A in the dot product is 1,
831 * call PZDOTU to perform the dot product.
832 *
833  IF( upper ) THEN
834  CALL pzdotu( j-1, csumj, a, ia, ja+j-1, desca, 1,
835  $ x, ix, jx, descx, 1 )
836  ELSE IF( j.LT.n ) THEN
837  CALL pzdotu( n-j, csumj, a, ia+j, ja+j-1, desca, 1,
838  $ x, ix+j, jx, descx, 1 )
839  END IF
840  IF( mycol.EQ.itmp2x ) THEN
841  CALL zgebs2d( contxt, 'Row', ' ', 1, 1, csumj, 1 )
842  ELSE
843  CALL zgebr2d( contxt, 'Row', ' ', 1, 1, csumj, 1,
844  $ myrow, itmp2x )
845  END IF
846  ELSE
847 *
848 * Otherwise, scale column of A by USCAL before dot
849 * product. Below is not the best way to do it.
850 *
851  IF( upper ) THEN
852 * DO 130 I = 1, J - 1
853 * CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
854 * 130 CONTINUE
855  zdum = dconjg( uscal )
856  CALL pzscal( j-1, zdum, a, ia, ja+j-1, desca, 1 )
857  CALL pzdotu( j-1, csumj, a, ia, ja+j-1, desca, 1,
858  $ x, ix, jx, descx, 1 )
859  zdum = zladiv( zdum, uscal )
860  CALL pzscal( j-1, zdum, a, ia, ja+j-1, desca, 1 )
861  ELSE IF( j.LT.n ) THEN
862 * DO 140 I = J + 1, N
863 * CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
864 * 140 CONTINUE
865  zdum = dconjg( uscal )
866  CALL pzscal( n-j, zdum, a, ia+j, ja+j-1, desca, 1 )
867  CALL pzdotu( n-j, csumj, a, ia+j, ja+j-1, desca, 1,
868  $ x, ix+j, jx, descx, 1 )
869  zdum = zladiv( zdum, uscal )
870  CALL pzscal( n-j, zdum, a, ia+j, ja+j-1, desca, 1 )
871  END IF
872  IF( mycol.EQ.itmp2x ) THEN
873  CALL zgebs2d( contxt, 'Row', ' ', 1, 1, csumj, 1 )
874  ELSE
875  CALL zgebr2d( contxt, 'Row', ' ', 1, 1, csumj, 1,
876  $ myrow, itmp2x )
877  END IF
878  END IF
879 *
880  IF( uscal.EQ.dcmplx( tscal ) ) THEN
881 *
882 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
883 * was not used to scale the dotproduct.
884 *
885 * X( J ) = X( J ) - CSUMJ
886 * XJ = CABS1( X( J ) )
887  xjtmp = xjtmp - csumj
888  xj = cabs1( xjtmp )
889 * IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
890 * $ X( IROWX ) = XJTMP
891  IF( nounit ) THEN
892 * TJJS = A( J, J )*TSCAL
893  CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
894  $ myrow, mycol, irow, icol, itmp1,
895  $ itmp2 )
896  IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
897  $ THEN
898  tjjs = a( ( icol-1 )*lda+irow )*tscal
899  CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs,
900  $ 1 )
901  ELSE
902  CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
903  $ itmp1, itmp2 )
904  END IF
905  ELSE
906  tjjs = tscal
907  IF( tscal.EQ.one )
908  $ GO TO 110
909  END IF
910 *
911 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
912 *
913  tjj = cabs1( tjjs )
914  IF( tjj.GT.smlnum ) THEN
915 *
916 * abs(A(j,j)) > SMLNUM:
917 *
918  IF( tjj.LT.one ) THEN
919  IF( xj.GT.tjj*bignum ) THEN
920 *
921 * Scale X by 1/abs(x(j)).
922 *
923  rec = one / xj
924  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
925  xjtmp = xjtmp*rec
926  scale = scale*rec
927  xmax = xmax*rec
928  END IF
929  END IF
930 * X( J ) = ZLADIV( X( J ), TJJS )
931  xjtmp = zladiv( xjtmp, tjjs )
932  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
933  $ THEN
934  x( irowx ) = xjtmp
935  END IF
936  ELSE IF( tjj.GT.zero ) THEN
937 *
938 * 0 < abs(A(j,j)) <= SMLNUM:
939 *
940  IF( xj.GT.tjj*bignum ) THEN
941 *
942 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
943 *
944  rec = ( tjj*bignum ) / xj
945  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
946  xjtmp = xjtmp*rec
947  scale = scale*rec
948  xmax = xmax*rec
949  END IF
950 * X( J ) = ZLADIV( X( J ), TJJS )
951  xjtmp = zladiv( xjtmp, tjjs )
952  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
953  $ THEN
954  x( irowx ) = xjtmp
955  END IF
956  ELSE
957 *
958 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
959 * scale = 0 and compute a solution to A**T *x = 0.
960 *
961  CALL pzlaset( ' ', n, 1, czero, czero, x, ix, jx,
962  $ descx )
963  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
964  $ THEN
965  x( irowx ) = cone
966  END IF
967  xjtmp = cone
968  scale = zero
969  xmax = zero
970  END IF
971  110 CONTINUE
972  ELSE
973 *
974 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
975 * product has already been divided by 1/A(j,j).
976 *
977 * X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
978  xjtmp = zladiv( xjtmp, tjjs ) - csumj
979  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
980  $ THEN
981  x( irowx ) = xjtmp
982  END IF
983  END IF
984  xmax = max( xmax, cabs1( xjtmp ) )
985  120 CONTINUE
986 *
987  ELSE
988 *
989 * Solve A**H * x = b
990 *
991  DO 140 j = jfirst, jlast, jinc
992 *
993 * Compute x(j) = b(j) - sum A(k,j)*x(k).
994 * k<>j
995 *
996  CALL infog2l( ix+j-1, jx, descx, nprow, npcol, myrow,
997  $ mycol, irowx, icolx, itmp1x, itmp2x )
998  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) ) THEN
999  xjtmp = x( irowx )
1000  CALL zgebs2d( contxt, 'All', ' ', 1, 1, xjtmp, 1 )
1001  ELSE
1002  CALL zgebr2d( contxt, 'All', ' ', 1, 1, xjtmp, 1,
1003  $ itmp1x, itmp2x )
1004  END IF
1005  xj = cabs1( xjtmp )
1006  uscal = tscal
1007  rec = one / max( xmax, one )
1008  IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
1009 *
1010 * If x(j) could overflow, scale x by 1/(2*XMAX).
1011 *
1012  rec = rec*half
1013  IF( nounit ) THEN
1014 * TJJS = DCONJG( A( J, J ) )*TSCAL
1015  CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
1016  $ myrow, mycol, irow, icol, itmp1,
1017  $ itmp2 )
1018  IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
1019  $ THEN
1020  tjjs = dconjg( a( ( icol-1 )*lda+irow ) )*tscal
1021  CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs,
1022  $ 1 )
1023  ELSE
1024  CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
1025  $ itmp1, itmp2 )
1026  END IF
1027  ELSE
1028  tjjs = tscal
1029  END IF
1030  tjj = cabs1( tjjs )
1031  IF( tjj.GT.one ) THEN
1032 *
1033 * Divide by A(j,j) when scaling x if A(j,j) > 1.
1034 *
1035  rec = min( one, rec*tjj )
1036  uscal = zladiv( uscal, tjjs )
1037  END IF
1038  IF( rec.LT.one ) THEN
1039  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
1040  xjtmp = xjtmp*rec
1041  scale = scale*rec
1042  xmax = xmax*rec
1043  END IF
1044  END IF
1045 *
1046  csumj = czero
1047  IF( uscal.EQ.cone ) THEN
1048 *
1049 * If the scaling needed for A in the dot product is 1,
1050 * call PZDOTC to perform the dot product.
1051 *
1052  IF( upper ) THEN
1053  CALL pzdotc( j-1, csumj, a, ia, ja+j-1, desca, 1,
1054  $ x, ix, jx, descx, 1 )
1055  ELSE IF( j.LT.n ) THEN
1056  CALL pzdotc( n-j, csumj, a, ia+j, ja+j-1, desca, 1,
1057  $ x, ix+j, jx, descx, 1 )
1058  END IF
1059  IF( mycol.EQ.itmp2x ) THEN
1060  CALL zgebs2d( contxt, 'Row', ' ', 1, 1, csumj, 1 )
1061  ELSE
1062  CALL zgebr2d( contxt, 'Row', ' ', 1, 1, csumj, 1,
1063  $ myrow, itmp2x )
1064  END IF
1065  ELSE
1066 *
1067 * Otherwise, scale column of A by USCAL before dot
1068 * product. Below is not the best way to do it.
1069 *
1070  IF( upper ) THEN
1071 * DO 180 I = 1, J - 1
1072 * CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )*
1073 * $ X( I )
1074 * 180 CONTINUE
1075  zdum = dconjg( uscal )
1076  CALL pzscal( j-1, zdum, a, ia, ja+j-1, desca, 1 )
1077  CALL pzdotc( j-1, csumj, a, ia, ja+j-1, desca, 1,
1078  $ x, ix, jx, descx, 1 )
1079  zdum = zladiv( cone, zdum )
1080  CALL pzscal( j-1, zdum, a, ia, ja+j-1, desca, 1 )
1081  ELSE IF( j.LT.n ) THEN
1082 * DO 190 I = J + 1, N
1083 * CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )*
1084 * $ X( I )
1085 * 190 CONTINUE
1086  zdum = dconjg( uscal )
1087  CALL pzscal( n-j, zdum, a, ia+j, ja+j-1, desca, 1 )
1088  CALL pzdotc( n-j, csumj, a, ia+j, ja+j-1, desca, 1,
1089  $ x, ix+j, jx, descx, 1 )
1090  zdum = zladiv( cone, zdum )
1091  CALL pzscal( n-j, zdum, a, ia+j, ja+j-1, desca, 1 )
1092  END IF
1093  IF( mycol.EQ.itmp2x ) THEN
1094  CALL zgebs2d( contxt, 'Row', ' ', 1, 1, csumj, 1 )
1095  ELSE
1096  CALL zgebr2d( contxt, 'Row', ' ', 1, 1, csumj, 1,
1097  $ myrow, itmp2x )
1098  END IF
1099  END IF
1100 *
1101  IF( uscal.EQ.dcmplx( tscal ) ) THEN
1102 *
1103 * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
1104 * was not used to scale the dotproduct.
1105 *
1106 * X( J ) = X( J ) - CSUMJ
1107 * XJ = CABS1( X( J ) )
1108  xjtmp = xjtmp - csumj
1109  xj = cabs1( xjtmp )
1110 * IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
1111 * $ X( IROWX ) = XJTMP
1112  IF( nounit ) THEN
1113 * TJJS = DCONJG( A( J, J ) )*TSCAL
1114  CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
1115  $ myrow, mycol, irow, icol, itmp1,
1116  $ itmp2 )
1117  IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
1118  $ THEN
1119  tjjs = dconjg( a( ( icol-1 )*lda+irow ) )*tscal
1120  CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs,
1121  $ 1 )
1122  ELSE
1123  CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
1124  $ itmp1, itmp2 )
1125  END IF
1126  ELSE
1127  tjjs = tscal
1128  IF( tscal.EQ.one )
1129  $ GO TO 130
1130  END IF
1131 *
1132 * Compute x(j) = x(j) / A(j,j), scaling if necessary.
1133 *
1134  tjj = cabs1( tjjs )
1135  IF( tjj.GT.smlnum ) THEN
1136 *
1137 * abs(A(j,j)) > SMLNUM:
1138 *
1139  IF( tjj.LT.one ) THEN
1140  IF( xj.GT.tjj*bignum ) THEN
1141 *
1142 * Scale X by 1/abs(x(j)).
1143 *
1144  rec = one / xj
1145  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
1146  xjtmp = xjtmp*rec
1147  scale = scale*rec
1148  xmax = xmax*rec
1149  END IF
1150  END IF
1151 * X( J ) = ZLADIV( X( J ), TJJS )
1152  xjtmp = zladiv( xjtmp, tjjs )
1153  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1154  $ x( irowx ) = xjtmp
1155  ELSE IF( tjj.GT.zero ) THEN
1156 *
1157 * 0 < abs(A(j,j)) <= SMLNUM:
1158 *
1159  IF( xj.GT.tjj*bignum ) THEN
1160 *
1161 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
1162 *
1163  rec = ( tjj*bignum ) / xj
1164  CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
1165  xjtmp = xjtmp*rec
1166  scale = scale*rec
1167  xmax = xmax*rec
1168  END IF
1169 * X( J ) = ZLADIV( X( J ), TJJS )
1170  xjtmp = zladiv( xjtmp, tjjs )
1171  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1172  $ x( irowx ) = xjtmp
1173  ELSE
1174 *
1175 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
1176 * scale = 0 and compute a solution to A**H *x = 0.
1177 *
1178  CALL pzlaset( ' ', n, 1, czero, czero, x, ix, jx,
1179  $ descx )
1180  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1181  $ x( irowx ) = cone
1182  xjtmp = cone
1183  scale = zero
1184  xmax = zero
1185  END IF
1186  130 CONTINUE
1187  ELSE
1188 *
1189 * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
1190 * product has already been divided by 1/A(j,j).
1191 *
1192 * X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
1193  xjtmp = zladiv( xjtmp, tjjs ) - csumj
1194  IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1195  $ x( irowx ) = xjtmp
1196  END IF
1197  xmax = max( xmax, cabs1( xjtmp ) )
1198  140 CONTINUE
1199  END IF
1200  scale = scale / tscal
1201  END IF
1202 *
1203 * Scale the column norms by 1/TSCAL for return.
1204 *
1205  IF( tscal.NE.one ) THEN
1206  CALL dscal( n, one / tscal, cnorm, 1 )
1207  END IF
1208 *
1209  RETURN
1210 *
1211 * End of PZLATTRS
1212 *
1213  END
max
#define max(A, B)
Definition: pcgemr.c:180
pzlattrs
subroutine pzlattrs(UPLO, TRANS, DIAG, NORMIN, N, A, IA, JA, DESCA, X, IX, JX, DESCX, SCALE, CNORM, INFO)
Definition: pzlattrs.f:3
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pdlabad
subroutine pdlabad(ICTXT, SMALL, LARGE)
Definition: pdlabad.f:2
pzlaset
subroutine pzlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pzblastst.f:7509
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
min
#define min(A, B)
Definition: pcgemr.c:181