SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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, zr, zi
275 COMPLEX*16 CSUMJ, TJJS, USCAL, XJTMP, ZDUM
276 DOUBLE PRECISION XMAX( 1 )
277* ..
278* .. External Functions ..
279 LOGICAL LSAME
280 INTEGER IDAMAX
281 DOUBLE PRECISION PDLAMCH
282 EXTERNAL lsame, idamax, pdlamch
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, dladiv
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( 1 ) = zero
395 CALL pzamax( n, zdum, imax, x, ix, jx, descx, 1 )
396 xmax( 1 ) = cabs2( zdum )
397 CALL dgsum2d( contxt, 'Row', ' ', 1, 1, xmax, 1, -1, -1 )
398 xbnd = xmax( 1 )
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( 1 ).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( 1 )
599 CALL pzdscal( n, scale, x, ix, jx, descx, 1 )
600 xmax( 1 ) = bignum
601 ELSE
602 xmax( 1 ) = xmax( 1 )*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( 1 ) = xmax( 1 )*rec
655 END IF
656 END IF
657* X( J ) = ZLADIV( X( J ), TJJS )
658* XJ = CABS1( X( J ) )
659 CALL dladiv( dble( xjtmp ), dimag( xjtmp ),
660 $ dble( tjjs ), dimag( tjjs ), zr, zi )
661 xjtmp = dcmplx( zr, zi )
662 xj = cabs1( xjtmp )
663 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
664 $ THEN
665 x( irowx ) = xjtmp
666 END IF
667 ELSE IF( tjj.GT.zero ) THEN
668*
669* 0 < abs(A(j,j)) <= SMLNUM:
670*
671 IF( xj.GT.tjj*bignum ) THEN
672*
673* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
674* to avoid overflow when dividing by A(j,j).
675*
676 rec = ( tjj*bignum ) / xj
677 IF( cnorm( j ).GT.one ) THEN
678*
679* Scale by 1/CNORM(j) to avoid overflow when
680* multiplying x(j) times column j.
681*
682 rec = rec / cnorm( j )
683 END IF
684 CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
685 xjtmp = xjtmp*rec
686 scale = scale*rec
687 xmax( 1 ) = xmax( 1 )*rec
688 END IF
689* X( J ) = ZLADIV( X( J ), TJJS )
690* XJ = CABS1( X( J ) )
691 CALL dladiv( dble( xjtmp ), dimag( xjtmp ),
692 $ dble( tjjs ), dimag( tjjs ), zr, zi )
693 xjtmp = dcmplx( zr, zi )
694 xj = cabs1( xjtmp )
695 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
696 $ THEN
697 x( irowx ) = xjtmp
698 END IF
699 ELSE
700*
701* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
702* scale = 0, and compute a solution to A*x = 0.
703*
704 CALL pzlaset( ' ', n, 1, czero, czero, x, ix, jx,
705 $ descx )
706 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
707 $ THEN
708 x( irowx ) = cone
709 END IF
710 xjtmp = cone
711 xj = one
712 scale = zero
713 xmax( 1 ) = zero
714 END IF
715 90 CONTINUE
716*
717* Scale x if necessary to avoid overflow when adding a
718* multiple of column j of A.
719*
720 IF( xj.GT.one ) THEN
721 rec = one / xj
722 IF( cnorm( j ).GT.( bignum-xmax( 1 ) )*rec ) THEN
723*
724* Scale x by 1/(2*abs(x(j))).
725*
726 rec = rec*half
727 CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
728 xjtmp = xjtmp*rec
729 scale = scale*rec
730 END IF
731 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax( 1 ) ) ) THEN
732*
733* Scale x by 1/2.
734*
735 CALL pzdscal( n, half, x, ix, jx, descx, 1 )
736 xjtmp = xjtmp*half
737 scale = scale*half
738 END IF
739*
740 IF( upper ) THEN
741 IF( j.GT.1 ) THEN
742*
743* Compute the update
744* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
745*
746 zdum = -xjtmp*tscal
747 CALL pzaxpy( j-1, zdum, a, ia, ja+j-1, desca, 1, x,
748 $ ix, jx, descx, 1 )
749 CALL pzamax( j-1, zdum, imax, x, ix, jx, descx, 1 )
750 xmax( 1 ) = cabs1( zdum )
751 CALL dgsum2d( contxt, 'Row', ' ', 1, 1, xmax, 1,
752 $ -1, -1 )
753 END IF
754 ELSE
755 IF( j.LT.n ) THEN
756*
757* Compute the update
758* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
759*
760 zdum = -xjtmp*tscal
761 CALL pzaxpy( n-j, zdum, a, ia+j, ja+j-1, desca, 1,
762 $ x, ix+j, jx, descx, 1 )
763 CALL pzamax( n-j, zdum, i, x, ix+j, jx, descx, 1 )
764 xmax( 1 ) = cabs1( zdum )
765 CALL dgsum2d( contxt, 'Row', ' ', 1, 1, xmax, 1,
766 $ -1, -1 )
767 END IF
768 END IF
769 100 CONTINUE
770*
771 ELSE IF( lsame( trans, 'T' ) ) THEN
772*
773* Solve A**T * x = b
774*
775 DO 120 j = jfirst, jlast, jinc
776*
777* Compute x(j) = b(j) - sum A(k,j)*x(k).
778* k<>j
779*
780* XJ = CABS1( X( J ) )
781 CALL infog2l( ix+j-1, jx, descx, nprow, npcol, myrow,
782 $ mycol, irowx, icolx, itmp1x, itmp2x )
783 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) ) THEN
784 xjtmp = x( irowx )
785 CALL zgebs2d( contxt, 'All', ' ', 1, 1, xjtmp, 1 )
786 ELSE
787 CALL zgebr2d( contxt, 'All', ' ', 1, 1, xjtmp, 1,
788 $ itmp1x, itmp2x )
789 END IF
790 xj = cabs1( xjtmp )
791 uscal = dcmplx( tscal )
792 rec = one / max( xmax( 1 ), one )
793 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
794*
795* If x(j) could overflow, scale x by 1/(2*XMAX).
796*
797 rec = rec*half
798 IF( nounit ) THEN
799* TJJS = A( J, J )*TSCAL
800 CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
801 $ myrow, mycol, irow, icol, itmp1,
802 $ itmp2 )
803 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
804 $ THEN
805 tjjs = a( ( icol-1 )*lda+irow )*tscal
806 CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs,
807 $ 1 )
808 ELSE
809 CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
810 $ itmp1, itmp2 )
811 END IF
812 ELSE
813 tjjs = tscal
814 END IF
815 tjj = cabs1( tjjs )
816 IF( tjj.GT.one ) THEN
817*
818* Divide by A(j,j) when scaling x if A(j,j) > 1.
819*
820 rec = min( one, rec*tjj )
821 CALL dladiv( dble( uscal ), dimag( uscal ),
822 $ dble( tjjs ), dimag( tjjs ), zr, zi )
823 uscal = dcmplx( zr, zi )
824 END IF
825 IF( rec.LT.one ) THEN
826 CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
827 xjtmp = xjtmp*rec
828 scale = scale*rec
829 xmax( 1 ) = xmax( 1 )*rec
830 END IF
831 END IF
832*
833 csumj = czero
834 IF( uscal.EQ.cone ) THEN
835*
836* If the scaling needed for A in the dot product is 1,
837* call PZDOTU to perform the dot product.
838*
839 IF( upper ) THEN
840 CALL pzdotu( j-1, csumj, a, ia, ja+j-1, desca, 1,
841 $ x, ix, jx, descx, 1 )
842 ELSE IF( j.LT.n ) THEN
843 CALL pzdotu( n-j, csumj, a, ia+j, ja+j-1, desca, 1,
844 $ x, ix+j, jx, descx, 1 )
845 END IF
846 IF( mycol.EQ.itmp2x ) THEN
847 CALL zgebs2d( contxt, 'Row', ' ', 1, 1, csumj, 1 )
848 ELSE
849 CALL zgebr2d( contxt, 'Row', ' ', 1, 1, csumj, 1,
850 $ myrow, itmp2x )
851 END IF
852 ELSE
853*
854* Otherwise, scale column of A by USCAL before dot
855* product. Below is not the best way to do it.
856*
857 IF( upper ) THEN
858* DO 130 I = 1, J - 1
859* CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
860* 130 CONTINUE
861 zdum = dconjg( uscal )
862 CALL pzscal( j-1, zdum, a, ia, ja+j-1, desca, 1 )
863 CALL pzdotu( j-1, csumj, a, ia, ja+j-1, desca, 1,
864 $ x, ix, jx, descx, 1 )
865 CALL dladiv( dble( zdum ), dimag( zdum ),
866 $ dble( uscal ), dimag( uscal ), zr, zi)
867 zdum = dcmplx( zr, zi )
868 CALL pzscal( j-1, zdum, a, ia, ja+j-1, desca, 1 )
869 ELSE IF( j.LT.n ) THEN
870* DO 140 I = J + 1, N
871* CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
872* 140 CONTINUE
873 zdum = dconjg( uscal )
874 CALL pzscal( n-j, zdum, a, ia+j, ja+j-1, desca, 1 )
875 CALL pzdotu( n-j, csumj, a, ia+j, ja+j-1, desca, 1,
876 $ x, ix+j, jx, descx, 1 )
877 CALL dladiv( dble( zdum ), dimag( zdum ),
878 $ dble( uscal ), dimag( uscal ), zr, zi)
879 zdum = dcmplx( zr, zi )
880 CALL pzscal( n-j, zdum, a, ia+j, ja+j-1, desca, 1 )
881 END IF
882 IF( mycol.EQ.itmp2x ) THEN
883 CALL zgebs2d( contxt, 'Row', ' ', 1, 1, csumj, 1 )
884 ELSE
885 CALL zgebr2d( contxt, 'Row', ' ', 1, 1, csumj, 1,
886 $ myrow, itmp2x )
887 END IF
888 END IF
889*
890 IF( uscal.EQ.dcmplx( tscal ) ) THEN
891*
892* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
893* was not used to scale the dotproduct.
894*
895* X( J ) = X( J ) - CSUMJ
896* XJ = CABS1( X( J ) )
897 xjtmp = xjtmp - csumj
898 xj = cabs1( xjtmp )
899* IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
900* $ X( IROWX ) = XJTMP
901 IF( nounit ) THEN
902* TJJS = A( J, J )*TSCAL
903 CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
904 $ myrow, mycol, irow, icol, itmp1,
905 $ itmp2 )
906 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
907 $ THEN
908 tjjs = a( ( icol-1 )*lda+irow )*tscal
909 CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs,
910 $ 1 )
911 ELSE
912 CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
913 $ itmp1, itmp2 )
914 END IF
915 ELSE
916 tjjs = tscal
917 IF( tscal.EQ.one )
918 $ GO TO 110
919 END IF
920*
921* Compute x(j) = x(j) / A(j,j), scaling if necessary.
922*
923 tjj = cabs1( tjjs )
924 IF( tjj.GT.smlnum ) THEN
925*
926* abs(A(j,j)) > SMLNUM:
927*
928 IF( tjj.LT.one ) THEN
929 IF( xj.GT.tjj*bignum ) THEN
930*
931* Scale X by 1/abs(x(j)).
932*
933 rec = one / xj
934 CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
935 xjtmp = xjtmp*rec
936 scale = scale*rec
937 xmax( 1 ) = xmax( 1 )*rec
938 END IF
939 END IF
940* X( J ) = ZLADIV( X( J ), TJJS )
941 CALL dladiv( dble( xjtmp ), dimag( xjtmp ),
942 $ dble( tjjs ), dimag( tjjs ), zr, zi )
943 xjtmp = dcmplx( zr, zi )
944 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
945 $ THEN
946 x( irowx ) = xjtmp
947 END IF
948 ELSE IF( tjj.GT.zero ) THEN
949*
950* 0 < abs(A(j,j)) <= SMLNUM:
951*
952 IF( xj.GT.tjj*bignum ) THEN
953*
954* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
955*
956 rec = ( tjj*bignum ) / xj
957 CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
958 xjtmp = xjtmp*rec
959 scale = scale*rec
960 xmax( 1 ) = xmax( 1 )*rec
961 END IF
962* X( J ) = ZLADIV( X( J ), TJJS )
963 CALL dladiv( dble( xjtmp ), dimag( xjtmp ),
964 $ dble( tjjs ), dimag( tjjs ), zr, zi )
965 xjtmp = dcmplx( zr, zi )
966 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
967 $ THEN
968 x( irowx ) = xjtmp
969 END IF
970 ELSE
971*
972* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
973* scale = 0 and compute a solution to A**T *x = 0.
974*
975 CALL pzlaset( ' ', n, 1, czero, czero, x, ix, jx,
976 $ descx )
977 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
978 $ THEN
979 x( irowx ) = cone
980 END IF
981 xjtmp = cone
982 scale = zero
983 xmax( 1 ) = zero
984 END IF
985 110 CONTINUE
986 ELSE
987*
988* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
989* product has already been divided by 1/A(j,j).
990*
991* X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
992 CALL dladiv( dble( xjtmp ), dimag( xjtmp ),
993 $ dble( tjjs ), dimag( tjjs ), zr, zi )
994 xjtmp = dcmplx( zr, zi )
995 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
996 $ THEN
997 x( irowx ) = xjtmp
998 END IF
999 END IF
1000 xmax( 1 ) = max( xmax( 1 ), cabs1( xjtmp ) )
1001 120 CONTINUE
1002*
1003 ELSE
1004*
1005* Solve A**H * x = b
1006*
1007 DO 140 j = jfirst, jlast, jinc
1008*
1009* Compute x(j) = b(j) - sum A(k,j)*x(k).
1010* k<>j
1011*
1012 CALL infog2l( ix+j-1, jx, descx, nprow, npcol, myrow,
1013 $ mycol, irowx, icolx, itmp1x, itmp2x )
1014 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) ) THEN
1015 xjtmp = x( irowx )
1016 CALL zgebs2d( contxt, 'All', ' ', 1, 1, xjtmp, 1 )
1017 ELSE
1018 CALL zgebr2d( contxt, 'All', ' ', 1, 1, xjtmp, 1,
1019 $ itmp1x, itmp2x )
1020 END IF
1021 xj = cabs1( xjtmp )
1022 uscal = tscal
1023 rec = one / max( xmax( 1 ), one )
1024 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
1025*
1026* If x(j) could overflow, scale x by 1/(2*XMAX).
1027*
1028 rec = rec*half
1029 IF( nounit ) THEN
1030* TJJS = DCONJG( A( J, J ) )*TSCAL
1031 CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
1032 $ myrow, mycol, irow, icol, itmp1,
1033 $ itmp2 )
1034 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
1035 $ THEN
1036 tjjs = dconjg( a( ( icol-1 )*lda+irow ) )*tscal
1037 CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs,
1038 $ 1 )
1039 ELSE
1040 CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
1041 $ itmp1, itmp2 )
1042 END IF
1043 ELSE
1044 tjjs = tscal
1045 END IF
1046 tjj = cabs1( tjjs )
1047 IF( tjj.GT.one ) THEN
1048*
1049* Divide by A(j,j) when scaling x if A(j,j) > 1.
1050*
1051 rec = min( one, rec*tjj )
1052 CALL dladiv( dble( uscal ), dimag( uscal ),
1053 $ dble( tjjs ), dimag( tjjs ), zr, zi )
1054 uscal = dcmplx( zr, zi )
1055 END IF
1056 IF( rec.LT.one ) THEN
1057 CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
1058 xjtmp = xjtmp*rec
1059 scale = scale*rec
1060 xmax( 1 ) = xmax( 1 )*rec
1061 END IF
1062 END IF
1063*
1064 csumj = czero
1065 IF( uscal.EQ.cone ) THEN
1066*
1067* If the scaling needed for A in the dot product is 1,
1068* call PZDOTC to perform the dot product.
1069*
1070 IF( upper ) THEN
1071 CALL pzdotc( j-1, csumj, a, ia, ja+j-1, desca, 1,
1072 $ x, ix, jx, descx, 1 )
1073 ELSE IF( j.LT.n ) THEN
1074 CALL pzdotc( n-j, csumj, a, ia+j, ja+j-1, desca, 1,
1075 $ x, ix+j, jx, descx, 1 )
1076 END IF
1077 IF( mycol.EQ.itmp2x ) THEN
1078 CALL zgebs2d( contxt, 'Row', ' ', 1, 1, csumj, 1 )
1079 ELSE
1080 CALL zgebr2d( contxt, 'Row', ' ', 1, 1, csumj, 1,
1081 $ myrow, itmp2x )
1082 END IF
1083 ELSE
1084*
1085* Otherwise, scale column of A by USCAL before dot
1086* product. Below is not the best way to do it.
1087*
1088 IF( upper ) THEN
1089* DO 180 I = 1, J - 1
1090* CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )*
1091* $ X( I )
1092* 180 CONTINUE
1093 zdum = dconjg( uscal )
1094 CALL pzscal( j-1, zdum, a, ia, ja+j-1, desca, 1 )
1095 CALL pzdotc( j-1, csumj, a, ia, ja+j-1, desca, 1,
1096 $ x, ix, jx, descx, 1 )
1097 CALL dladiv( one, zero,
1098 $ dble( zdum ), dimag( zdum ), zr, zi )
1099 zdum = dcmplx( zr, zi )
1100 CALL pzscal( j-1, zdum, a, ia, ja+j-1, desca, 1 )
1101 ELSE IF( j.LT.n ) THEN
1102* DO 190 I = J + 1, N
1103* CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )*
1104* $ X( I )
1105* 190 CONTINUE
1106 zdum = dconjg( uscal )
1107 CALL pzscal( n-j, zdum, a, ia+j, ja+j-1, desca, 1 )
1108 CALL pzdotc( n-j, csumj, a, ia+j, ja+j-1, desca, 1,
1109 $ x, ix+j, jx, descx, 1 )
1110 CALL dladiv( one, zero,
1111 $ dble( zdum ), dimag( zdum ), zr, zi )
1112 zdum = dcmplx( zr, zi )
1113 CALL pzscal( n-j, zdum, a, ia+j, ja+j-1, desca, 1 )
1114 END IF
1115 IF( mycol.EQ.itmp2x ) THEN
1116 CALL zgebs2d( contxt, 'Row', ' ', 1, 1, csumj, 1 )
1117 ELSE
1118 CALL zgebr2d( contxt, 'Row', ' ', 1, 1, csumj, 1,
1119 $ myrow, itmp2x )
1120 END IF
1121 END IF
1122*
1123 IF( uscal.EQ.dcmplx( tscal ) ) THEN
1124*
1125* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
1126* was not used to scale the dotproduct.
1127*
1128* X( J ) = X( J ) - CSUMJ
1129* XJ = CABS1( X( J ) )
1130 xjtmp = xjtmp - csumj
1131 xj = cabs1( xjtmp )
1132* IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
1133* $ X( IROWX ) = XJTMP
1134 IF( nounit ) THEN
1135* TJJS = DCONJG( A( J, J ) )*TSCAL
1136 CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol,
1137 $ myrow, mycol, irow, icol, itmp1,
1138 $ itmp2 )
1139 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
1140 $ THEN
1141 tjjs = dconjg( a( ( icol-1 )*lda+irow ) )*tscal
1142 CALL zgebs2d( contxt, 'All', ' ', 1, 1, tjjs,
1143 $ 1 )
1144 ELSE
1145 CALL zgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
1146 $ itmp1, itmp2 )
1147 END IF
1148 ELSE
1149 tjjs = tscal
1150 IF( tscal.EQ.one )
1151 $ GO TO 130
1152 END IF
1153*
1154* Compute x(j) = x(j) / A(j,j), scaling if necessary.
1155*
1156 tjj = cabs1( tjjs )
1157 IF( tjj.GT.smlnum ) THEN
1158*
1159* abs(A(j,j)) > SMLNUM:
1160*
1161 IF( tjj.LT.one ) THEN
1162 IF( xj.GT.tjj*bignum ) THEN
1163*
1164* Scale X by 1/abs(x(j)).
1165*
1166 rec = one / xj
1167 CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
1168 xjtmp = xjtmp*rec
1169 scale = scale*rec
1170 xmax( 1 ) = xmax( 1 )*rec
1171 END IF
1172 END IF
1173* X( J ) = ZLADIV( X( J ), TJJS )
1174 CALL dladiv( dble( xjtmp ), dimag( xjtmp ),
1175 $ dble( tjjs ), dimag( tjjs ), zr, zi )
1176 xjtmp = dcmplx( zr, zi )
1177 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1178 $ x( irowx ) = xjtmp
1179 ELSE IF( tjj.GT.zero ) THEN
1180*
1181* 0 < abs(A(j,j)) <= SMLNUM:
1182*
1183 IF( xj.GT.tjj*bignum ) THEN
1184*
1185* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
1186*
1187 rec = ( tjj*bignum ) / xj
1188 CALL pzdscal( n, rec, x, ix, jx, descx, 1 )
1189 xjtmp = xjtmp*rec
1190 scale = scale*rec
1191 xmax( 1 ) = xmax( 1 )*rec
1192 END IF
1193* X( J ) = ZLADIV( X( J ), TJJS )
1194 CALL dladiv( dble( xjtmp ), dimag( xjtmp ),
1195 $ dble( tjjs ), dimag( tjjs ), zr, zi )
1196 xjtmp = dcmplx( zr, zi )
1197 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1198 $ x( irowx ) = xjtmp
1199 ELSE
1200*
1201* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
1202* scale = 0 and compute a solution to A**H *x = 0.
1203*
1204 CALL pzlaset( ' ', n, 1, czero, czero, x, ix, jx,
1205 $ descx )
1206 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1207 $ x( irowx ) = cone
1208 xjtmp = cone
1209 scale = zero
1210 xmax( 1 ) = zero
1211 END IF
1212 130 CONTINUE
1213 ELSE
1214*
1215* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
1216* product has already been divided by 1/A(j,j).
1217*
1218* X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
1219 CALL dladiv( dble( xjtmp ), dimag( xjtmp ),
1220 $ dble( tjjs ), dimag( tjjs ), zr, zi )
1221 xjtmp = dcmplx( zr, zi )
1222 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1223 $ x( irowx ) = xjtmp
1224 END IF
1225 xmax( 1 ) = max( xmax( 1 ), cabs1( xjtmp ) )
1226 140 CONTINUE
1227 END IF
1228 scale = scale / tscal
1229 END IF
1230*
1231* Scale the column norms by 1/TSCAL for return.
1232*
1233 IF( tscal.NE.one ) THEN
1234 CALL dscal( n, one / tscal, cnorm, 1 )
1235 END IF
1236*
1237 RETURN
1238*
1239* End of PZLATTRS
1240*
1241 END
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
pdlabad
subroutine pdlabad(ictxt, small, large)
Definition pdlabad.f:2
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2
pzlaset
subroutine pzlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pzblastst.f:7509
pzlattrs
subroutine pzlattrs(uplo, trans, diag, normin, n, a, ia, ja, desca, x, ix, jx, descx, scale, cnorm, info)
Definition pzlattrs.f:3