OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pclattrs.f
Go to the documentation of this file.
1 SUBROUTINE pclattrs( 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 REAL SCALE
13* ..
14* .. Array Arguments ..
15 INTEGER DESCA( * ), DESCX( * )
16 REAL CNORM( * )
17 COMPLEX A( * ), X( * )
18* ..
19*
20* Purpose
21* =======
22*
23* PCLATTRS 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* PCTRSV 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 PCTRSV. 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 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 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) REAL
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) REAL 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, PCTRSV
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 PCTRSV 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 PCTRSV 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 REAL ZERO, HALF, ONE, TWO
256 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
257 $ two = 2.0e+0 )
258 COMPLEX CZERO, CONE
259 parameter( czero = ( 0.0e+0, 0.0e+0 ),
260 $ cone = ( 1.0e+0, 0.0e+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 REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
274 $ xbnd, xj
275 REAL XMAX( 1 )
276 COMPLEX CSUMJ, TJJS, USCAL, XJTMP, ZDUM
277* ..
278* .. External Functions ..
279 LOGICAL LSAME
280 INTEGER ISAMAX
281 REAL PSLAMCH
282 COMPLEX CLADIV
283 EXTERNAL lsame, isamax, pslamch, cladiv
284* ..
285* .. External Subroutines ..
286 EXTERNAL blacs_gridinfo, sgsum2d, sscal, infog2l,
287 $ pscasum, pslabad, pxerbla, pcamax, pcaxpy,
288 $ pcdotc, pcdotu, pcsscal, pclaset, pcscal,
289 $ pctrsv, cgebr2d, cgebs2d
290* ..
291* .. Intrinsic Functions ..
292 INTRINSIC abs, real, cmplx, conjg, aimag, max, min
293* ..
294* .. Statement Functions ..
295 REAL CABS1, CABS2
296* ..
297* .. Statement Function definitions ..
298 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
299 cabs2( zdum ) = abs( real( zdum ) / 2.e0 ) +
300 $ abs( aimag( zdum ) / 2.e0 )
301* ..
302* .. Executable Statements ..
303*
304 info = 0
305 upper = lsame( uplo, 'U' )
306 notran = lsame( trans, 'N' )
307 nounit = lsame( diag, 'N' )
308*
309 contxt = desca( ctxt_ )
310 rsrc = desca( rsrc_ )
311 csrc = desca( csrc_ )
312 mb = desca( mb_ )
313 nb = desca( nb_ )
314 lda = desca( lld_ )
315 ldx = descx( lld_ )
316*
317* Test the input parameters.
318*
319 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
320 info = -1
321 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
322 $ lsame( trans, 'C' ) ) THEN
323 info = -2
324 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
325 info = -3
326 ELSE IF( .NOT.lsame( normin, 'Y' ) .AND. .NOT.
327 $ lsame( normin, 'N' ) ) THEN
328 info = -4
329 ELSE IF( n.LT.0 ) THEN
330 info = -5
331 END IF
332*
333 CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
334*
335 IF( info.NE.0 ) THEN
336 CALL pxerbla( contxt, 'PCLATTRS', -info )
337 RETURN
338 END IF
339*
340* Quick return if possible
341*
342 IF( n.EQ.0 )
343 $ RETURN
344*
345* Determine machine dependent parameters to control overflow.
346*
347 smlnum = pslamch( contxt, 'Safe minimum' )
348 bignum = one / smlnum
349 CALL pslabad( contxt, smlnum, bignum )
350 smlnum = smlnum / pslamch( contxt, 'Precision' )
351 bignum = one / smlnum
352 scale = one
353*
354*
355 IF( lsame( normin, 'N' ) ) THEN
356*
357* Compute the 1-norm of each column, not including the diagonal.
358*
359 IF( upper ) THEN
360*
361* A is upper triangular.
362*
363 cnorm( 1 ) = zero
364 DO 10 j = 2, n
365 CALL pscasum( j-1, cnorm( j ), a, ia, ja+j-1, desca, 1 )
366 10 CONTINUE
367 ELSE
368*
369* A is lower triangular.
370*
371 DO 20 j = 1, n - 1
372 CALL pscasum( n-j, cnorm( j ), a, ia+j, ja+j-1, desca,
373 $ 1 )
374 20 CONTINUE
375 cnorm( n ) = zero
376 END IF
377 CALL sgsum2d( contxt, 'Row', ' ', n, 1, cnorm, 1, -1, -1 )
378 END IF
379*
380* Scale the column norms by TSCAL if the maximum element in CNORM is
381* greater than BIGNUM/2.
382*
383 imax = isamax( n, cnorm, 1 )
384 tmax = cnorm( imax )
385 IF( tmax.LE.bignum*half ) THEN
386 tscal = one
387 ELSE
388 tscal = half / ( smlnum*tmax )
389 CALL sscal( n, tscal, cnorm, 1 )
390 END IF
391*
392* Compute a bound on the computed solution vector to see if the
393* Level 2 PBLAS routine PCTRSV can be used.
394*
395 xmax( 1 ) = zero
396 CALL pcamax( n, zdum, imax, x, ix, jx, descx, 1 )
397 xmax( 1 ) = cabs2( zdum )
398 CALL sgsum2d( contxt, 'Row', ' ', 1, 1, xmax, 1, -1, -1 )
399 xbnd = xmax( 1 )
400*
401 IF( notran ) THEN
402*
403* Compute the growth in A * x = b.
404*
405 IF( upper ) THEN
406 jfirst = n
407 jlast = 1
408 jinc = -1
409 ELSE
410 jfirst = 1
411 jlast = n
412 jinc = 1
413 END IF
414*
415 IF( tscal.NE.one ) THEN
416 grow = zero
417 GO TO 50
418 END IF
419*
420 IF( nounit ) THEN
421*
422* A is non-unit triangular.
423*
424* Compute GROW = 1/G(j) and XBND = 1/M(j).
425* Initially, G(0) = max{x(i), i=1,...,n}.
426*
427 grow = half / max( xbnd, smlnum )
428 xbnd = grow
429 DO 30 j = jfirst, jlast, jinc
430*
431* Exit the loop if the growth factor is too small.
432*
433 IF( grow.LE.smlnum )
434 $ GO TO 50
435*
436* TJJS = A( J, J )
437 CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol, myrow,
438 $ mycol, irow, icol, itmp1, itmp2 )
439 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
440 tjjs = a( ( icol-1 )*lda+irow )
441 CALL cgebs2d( contxt, 'All', ' ', 1, 1, tjjs, 1 )
442 ELSE
443 CALL cgebr2d( contxt, 'All', ' ', 1, 1, tjjs, 1,
444 $ itmp1, itmp2 )
445 END IF
446 tjj = cabs1( tjjs )
447*
448 IF( tjj.GE.smlnum ) THEN
449*
450* M(j) = G(j-1) / abs(A(j,j))
451*
452 xbnd = min( xbnd, min( one, tjj )*grow )
453 ELSE
454*
455* M(j) could overflow, set XBND to 0.
456*
457 xbnd = zero
458 END IF
459*
460 IF( tjj+cnorm( j ).GE.smlnum ) THEN
461*
462* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
463*
464 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
465 ELSE
466*
467* G(j) could overflow, set GROW to 0.
468*
469 grow = zero
470 END IF
471 30 CONTINUE
472 grow = xbnd
473 ELSE
474*
475* A is unit triangular.
476*
477* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
478*
479 grow = min( one, half / max( xbnd, smlnum ) )
480 DO 40 j = jfirst, jlast, jinc
481*
482* Exit the loop if the growth factor is too small.
483*
484 IF( grow.LE.smlnum )
485 $ GO TO 50
486*
487* G(j) = G(j-1)*( 1 + CNORM(j) )
488*
489 grow = grow*( one / ( one+cnorm( j ) ) )
490 40 CONTINUE
491 END IF
492 50 CONTINUE
493*
494 ELSE
495*
496* Compute the growth in A**T * x = b or A**H * x = b.
497*
498 IF( upper ) THEN
499 jfirst = 1
500 jlast = n
501 jinc = 1
502 ELSE
503 jfirst = n
504 jlast = 1
505 jinc = -1
506 END IF
507*
508 IF( tscal.NE.one ) THEN
509 grow = zero
510 GO TO 80
511 END IF
512*
513 IF( nounit ) THEN
514*
515* A is non-unit triangular.
516*
517* Compute GROW = 1/G(j) and XBND = 1/M(j).
518* Initially, M(0) = max{x(i), i=1,...,n}.
519*
520 grow = half / max( xbnd, smlnum )
521 xbnd = grow
522 DO 60 j = jfirst, jlast, jinc
523*
524* Exit the loop if the growth factor is too small.
525*
526 IF( grow.LE.smlnum )
527 $ GO TO 80
528*
529* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
530*
531 xj = one + cnorm( j )
532 grow = min( grow, xbnd / xj )
533*
534* TJJS = A( J, J )
535 CALL infog2l( ia+j-1, ja+j-1, desca, nprow, npcol, myrow,
536 $ mycol, irow, icol, itmp1, itmp2 )
537 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
538 tjjs = a( ( icol-1 )*lda+irow )
539 CALL cgebs2d( contxt, 'All', ' ', 1, 1, tjjs, 1 )
540 ELSE
541 CALL cgebr2d( contxt, 'all', ' ', 1, 1, tjjs, 1,
542 $ itmp1, itmp2 )
543 END IF
544 tjj = cabs1( tjjs )
545*
546 IF( tjj.GE.smlnum ) THEN
547*
548* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
549*
550 IF( xj.GT.tjj )
551 $ xbnd = xbnd*( tjj / xj )
552 ELSE
553*
554* M(j) could overflow, set XBND to 0.
555*
556 xbnd = zero
557 END IF
558 60 CONTINUE
559 grow = min( grow, xbnd )
560 ELSE
561*
562* A is unit triangular.
563*
564* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
565*
566 grow = min( one, half / max( xbnd, smlnum ) )
567 DO 70 j = jfirst, jlast, jinc
568*
569* Exit the loop if the growth factor is too small.
570*
571 IF( grow.LE.smlnum )
572 $ GO TO 80
573*
574* G(j) = ( 1 + CNORM(j) )*G(j-1)
575*
576 xj = one + cnorm( j )
577 grow = grow / xj
578 70 CONTINUE
579 END IF
580 80 CONTINUE
581 END IF
582*
583 IF( ( grow*tscal ).GT.smlnum ) THEN
584*
585* Use the Level 2 PBLAS solve if the reciprocal of the bound on
586* elements of X is not too small.
587*
588 CALL pctrsv( uplo, trans, diag, n, a, ia, ja, desca, x, ix, jx,
589 $ descx, 1 )
590 ELSE
591*
592* Use a Level 1 PBLAS solve, scaling intermediate results.
593*
594 IF( xmax( 1 ).GT.bignum*half ) THEN
595*
596* Scale X so that its components are less than or equal to
597* BIGNUM in absolute value.
598*
599 scale = ( bignum*half ) / xmax( 1 )
600 CALL pcsscal( n, scale, x, ix, jx, descx, 1 )
601 xmax( 1 ) = bignum
602 ELSE
603 xmax( 1 ) = xmax( 1 )*two
604 END IF
605*
606 IF( notran ) THEN
607*
608* Solve A * x = b
609*
610 DO 100 j = jfirst, jlast, jinc
611*
612* Compute x(j) = b(j) / A(j,j), scaling x if necessary.
613*
614* XJ = CABS1( X( J ) )
615 CALL infog2l( ix+j-1, jx, descx, nprow, npcol, myrow,
616 $ mycol, irowx, icolx, itmp1x, itmp2x )
617 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) ) THEN
618 xjtmp = x( irowx )
619 CALL cgebs2d( contxt, 'All', ' ', 1, 1, xjtmp, 1 )
620 ELSE
621 CALL cgebr2d( contxt, 'All', ' ', 1, 1, XJTMP, 1,
622 $ ITMP1X, ITMP2X )
623 END IF
624 XJ = CABS1( XJTMP )
625 IF( NOUNIT ) THEN
626* TJJS = A( J, J )*TSCAL
627 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
628 $ MYROW, MYCOL, IROW, ICOL, ITMP1, ITMP2 )
629.EQ..AND..EQ. IF( ( MYROWITMP1 ) ( MYCOLITMP2 ) ) THEN
630 TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL
631 CALL CGEBS2D( CONTXT, 'all', ' ', 1, 1, TJJS, 1 )
632 ELSE
633 CALL CGEBR2D( CONTXT, 'all', ' ', 1, 1, TJJS, 1,
634 $ ITMP1, ITMP2 )
635 END IF
636 ELSE
637 TJJS = TSCAL
638.EQ. IF( TSCALONE )
639 $ GO TO 90
640 END IF
641 TJJ = CABS1( TJJS )
642.GT. IF( TJJSMLNUM ) THEN
643*
644* abs(A(j,j)) > SMLNUM:
645*
646.LT. IF( TJJONE ) THEN
647.GT. IF( XJTJJ*BIGNUM ) THEN
648*
649* Scale x by 1/b(j).
650*
651 REC = ONE / XJ
652 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
653 XJTMP = XJTMP*REC
654 SCALE = SCALE*REC
655 XMAX( 1 ) = XMAX( 1 )*REC
656 END IF
657 END IF
658* X( J ) = CLADIV( X( J ), TJJS )
659* XJ = CABS1( X( J ) )
660 XJTMP = CLADIV( XJTMP, TJJS )
661 XJ = CABS1( XJTMP )
662.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
663 $ THEN
664 X( IROWX ) = XJTMP
665 END IF
666.GT. ELSE IF( TJJZERO ) THEN
667*
668* 0 < abs(A(j,j)) <= SMLNUM:
669*
670.GT. IF( XJTJJ*BIGNUM ) THEN
671*
672* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
673* to avoid overflow when dividing by A(j,j).
674*
675 REC = ( TJJ*BIGNUM ) / XJ
676.GT. IF( CNORM( J )ONE ) THEN
677*
678* Scale by 1/CNORM(j) to avoid overflow when
679* multiplying x(j) times column j.
680*
681 REC = REC / CNORM( J )
682 END IF
683 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
684 XJTMP = XJTMP*REC
685 SCALE = SCALE*REC
686 XMAX( 1 ) = XMAX( 1 )*REC
687 END IF
688* X( J ) = CLADIV( X( J ), TJJS )
689* XJ = CABS1( X( J ) )
690 XJTMP = CLADIV( XJTMP, TJJS )
691 XJ = CABS1( XJTMP )
692.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
693 $ THEN
694 X( IROWX ) = XJTMP
695 END IF
696 ELSE
697*
698* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
699* scale = 0, and compute a solution to A*x = 0.
700*
701 CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX,
702 $ DESCX )
703.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
704 $ THEN
705 X( IROWX ) = CONE
706 END IF
707 XJTMP = CONE
708 XJ = ONE
709 SCALE = ZERO
710 XMAX( 1 ) = ZERO
711 END IF
712 90 CONTINUE
713*
714* Scale x if necessary to avoid overflow when adding a
715* multiple of column j of A.
716*
717.GT. IF( XJONE ) THEN
718 REC = ONE / XJ
719.GT. IF( CNORM( J )( BIGNUM-XMAX( 1 ) )*REC ) THEN
720*
721* Scale x by 1/(2*abs(x(j))).
722*
723 REC = REC*HALF
724 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
725 XJTMP = XJTMP*REC
726 SCALE = SCALE*REC
727 END IF
728.GT. ELSE IF( XJ*CNORM( J )( BIGNUM-XMAX( 1 ) ) ) THEN
729*
730* Scale x by 1/2.
731*
732 CALL PCSSCAL( N, HALF, X, IX, JX, DESCX, 1 )
733 XJTMP = XJTMP*HALF
734 SCALE = SCALE*HALF
735 END IF
736*
737 IF( UPPER ) THEN
738.GT. IF( J1 ) THEN
739*
740* Compute the update
741* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
742*
743 ZDUM = -XJTMP*TSCAL
744 CALL PCAXPY( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1, X,
745 $ IX, JX, DESCX, 1 )
746 CALL PCAMAX( J-1, ZDUM, IMAX, X, IX, JX, DESCX, 1 )
747 XMAX( 1 ) = CABS1( ZDUM )
748 CALL SGSUM2D( CONTXT, 'row', ' ', 1, 1, XMAX, 1,
749 $ -1, -1 )
750 END IF
751 ELSE
752.LT. IF( JN ) THEN
753*
754* Compute the update
755* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
756*
757 ZDUM = -XJTMP*TSCAL
758 CALL PCAXPY( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1,
759 $ X, IX+J, JX, DESCX, 1 )
760 CALL PCAMAX( N-J, ZDUM, I, X, IX+J, JX, DESCX, 1 )
761 XMAX( 1 ) = CABS1( ZDUM )
762 CALL SGSUM2D( CONTXT, 'row', ' ', 1, 1, XMAX, 1,
763 $ -1, -1 )
764 END IF
765 END IF
766 100 CONTINUE
767*
768 ELSE IF( LSAME( TRANS, 't' ) ) THEN
769*
770* Solve A**T * x = b
771*
772 DO 120 J = JFIRST, JLAST, JINC
773*
774* Compute x(j) = b(j) - sum A(k,j)*x(k).
775* k<>j
776*
777* XJ = CABS1( X( J ) )
778 CALL INFOG2L( IX+J-1, JX, DESCX, NPROW, NPCOL, MYROW,
779 $ MYCOL, IROWX, ICOLX, ITMP1X, ITMP2X )
780.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) ) THEN
781 XJTMP = X( IROWX )
782 CALL CGEBS2D( CONTXT, 'all', ' ', 1, 1, XJTMP, 1 )
783 ELSE
784 CALL CGEBR2D( CONTXT, 'all', ' ', 1, 1, XJTMP, 1,
785 $ ITMP1X, ITMP2X )
786 END IF
787 XJ = CABS1( XJTMP )
788 USCAL = CMPLX( TSCAL )
789 REC = ONE / MAX( XMAX( 1 ), ONE )
790.GT. IF( CNORM( J )( BIGNUM-XJ )*REC ) THEN
791*
792* If x(j) could overflow, scale x by 1/(2*XMAX).
793*
794 REC = REC*HALF
795 IF( NOUNIT ) THEN
796* TJJS = A( J, J )*TSCAL
797 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
798 $ MYROW, MYCOL, IROW, ICOL, ITMP1,
799 $ ITMP2 )
800.EQ..AND..EQ. IF( ( MYROWITMP1 ) ( MYCOLITMP2 ) )
801 $ THEN
802 TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL
803 CALL CGEBS2D( CONTXT, 'all', ' ', 1, 1, TJJS,
804 $ 1 )
805 ELSE
806 CALL CGEBR2D( CONTXT, 'all', ' ', 1, 1, TJJS, 1,
807 $ ITMP1, ITMP2 )
808 END IF
809 ELSE
810 TJJS = TSCAL
811 END IF
812 TJJ = CABS1( TJJS )
813.GT. IF( TJJONE ) THEN
814*
815* Divide by A(j,j) when scaling x if A(j,j) > 1.
816*
817 REC = MIN( ONE, REC*TJJ )
818 USCAL = CLADIV( USCAL, TJJS )
819 END IF
820.LT. IF( RECONE ) THEN
821 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
822 XJTMP = XJTMP*REC
823 SCALE = SCALE*REC
824 XMAX( 1 ) = XMAX( 1 )*REC
825 END IF
826 END IF
827*
828 CSUMJ = CZERO
829.EQ. IF( USCALCONE ) THEN
830*
831* If the scaling needed for A in the dot product is 1,
832* call PCDOTU to perform the dot product.
833*
834 IF( UPPER ) THEN
835 CALL PCDOTU( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1,
836 $ X, IX, JX, DESCX, 1 )
837.LT. ELSE IF( JN ) THEN
838 CALL PCDOTU( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1,
839 $ X, IX+J, JX, DESCX, 1 )
840 END IF
841.EQ. IF( MYCOLITMP2X ) THEN
842 CALL CGEBS2D( CONTXT, 'row', ' ', 1, 1, CSUMJ, 1 )
843 ELSE
844 CALL CGEBR2D( CONTXT, 'row', ' ', 1, 1, CSUMJ, 1,
845 $ MYROW, ITMP2X )
846 END IF
847 ELSE
848*
849* Otherwise, scale column of A by USCAL before dot
850* product. Below is not the best way to do it.
851*
852 IF( UPPER ) THEN
853* DO 130 I = 1, J - 1
854* CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
855* 130 CONTINUE
856 ZDUM = CONJG( USCAL )
857 CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 )
858 CALL PCDOTU( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1,
859 $ X, IX, JX, DESCX, 1 )
860 ZDUM = CLADIV( ZDUM, USCAL )
861 CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 )
862.LT. ELSE IF( JN ) THEN
863* DO 140 I = J + 1, N
864* CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
865* 140 CONTINUE
866 ZDUM = CONJG( USCAL )
867 CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 )
868 CALL PCDOTU( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1,
869 $ X, IX+J, JX, DESCX, 1 )
870 ZDUM = CLADIV( ZDUM, USCAL )
871 CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 )
872 END IF
873.EQ. IF( MYCOLITMP2X ) THEN
874 CALL CGEBS2D( CONTXT, 'row', ' ', 1, 1, CSUMJ, 1 )
875 ELSE
876 CALL CGEBR2D( CONTXT, 'row', ' ', 1, 1, CSUMJ, 1,
877 $ MYROW, ITMP2X )
878 END IF
879 END IF
880*
881.EQ. IF( USCALCMPLX( TSCAL ) ) THEN
882*
883* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
884* was not used to scale the dotproduct.
885*
886* X( J ) = X( J ) - CSUMJ
887* XJ = CABS1( X( J ) )
888 XJTMP = XJTMP - CSUMJ
889 XJ = CABS1( XJTMP )
890* IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
891* $ X( IROWX ) = XJTMP
892 IF( NOUNIT ) THEN
893* TJJS = A( J, J )*TSCAL
894 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
895 $ MYROW, MYCOL, IROW, ICOL, ITMP1,
896 $ ITMP2 )
897.EQ..AND..EQ. IF( ( MYROWITMP1 ) ( MYCOLITMP2 ) )
898 $ THEN
899 TJJS = A( ( ICOL-1 )*LDA+IROW )*TSCAL
900 CALL CGEBS2D( CONTXT, 'all', ' ', 1, 1, TJJS,
901 $ 1 )
902 ELSE
903 CALL CGEBR2D( CONTXT, 'all', ' ', 1, 1, TJJS, 1,
904 $ ITMP1, ITMP2 )
905 END IF
906 ELSE
907 TJJS = TSCAL
908.EQ. IF( TSCALONE )
909 $ GO TO 110
910 END IF
911*
912* Compute x(j) = x(j) / A(j,j), scaling if necessary.
913*
914 TJJ = CABS1( TJJS )
915.GT. IF( TJJSMLNUM ) THEN
916*
917* abs(A(j,j)) > SMLNUM:
918*
919.LT. IF( TJJONE ) THEN
920.GT. IF( XJTJJ*BIGNUM ) THEN
921*
922* Scale X by 1/abs(x(j)).
923*
924 REC = ONE / XJ
925 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
926 XJTMP = XJTMP*REC
927 SCALE = SCALE*REC
928 XMAX( 1 ) = XMAX( 1 )*REC
929 END IF
930 END IF
931* X( J ) = CLADIV( X( J ), TJJS )
932 XJTMP = CLADIV( XJTMP, TJJS )
933.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
934 $ THEN
935 X( IROWX ) = XJTMP
936 END IF
937.GT. ELSE IF( TJJZERO ) THEN
938*
939* 0 < abs(A(j,j)) <= SMLNUM:
940*
941.GT. IF( XJTJJ*BIGNUM ) THEN
942*
943* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
944*
945 REC = ( TJJ*BIGNUM ) / XJ
946 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
947 XJTMP = XJTMP*REC
948 SCALE = SCALE*REC
949 XMAX( 1 ) = XMAX( 1 )*REC
950 END IF
951* X( J ) = CLADIV( X( J ), TJJS )
952 XJTMP = CLADIV( XJTMP, TJJS )
953.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
954 $ THEN
955 X( IROWX ) = XJTMP
956 END IF
957 ELSE
958*
959* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
960* scale = 0 and compute a solution to A**T *x = 0.
961*
962 CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX,
963 $ DESCX )
964.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
965 $ THEN
966 X( IROWX ) = CONE
967 END IF
968 XJTMP = CONE
969 SCALE = ZERO
970 XMAX( 1 ) = ZERO
971 END IF
972 110 CONTINUE
973 ELSE
974*
975* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
976* product has already been divided by 1/A(j,j).
977*
978* X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
979 XJTMP = CLADIV( XJTMP, TJJS ) - CSUMJ
980.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
981 $ THEN
982 X( IROWX ) = XJTMP
983 END IF
984 END IF
985 XMAX( 1 ) = MAX( XMAX( 1 ), CABS1( XJTMP ) )
986 120 CONTINUE
987*
988 ELSE
989*
990* Solve A**H * x = b
991*
992 DO 140 J = JFIRST, JLAST, JINC
993*
994* Compute x(j) = b(j) - sum A(k,j)*x(k).
995* k<>j
996*
997 CALL INFOG2L( IX+J-1, JX, DESCX, NPROW, NPCOL, MYROW,
998 $ MYCOL, IROWX, ICOLX, ITMP1X, ITMP2X )
999.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) ) THEN
1000 XJTMP = X( IROWX )
1001 CALL CGEBS2D( CONTXT, 'all', ' ', 1, 1, XJTMP, 1 )
1002 ELSE
1003 CALL CGEBR2D( CONTXT, 'all', ' ', 1, 1, XJTMP, 1,
1004 $ ITMP1X, ITMP2X )
1005 END IF
1006 XJ = CABS1( XJTMP )
1007 USCAL = TSCAL
1008 REC = ONE / MAX( XMAX( 1 ), ONE )
1009.GT. IF( CNORM( J )( BIGNUM-XJ )*REC ) THEN
1010*
1011* If x(j) could overflow, scale x by 1/(2*XMAX).
1012*
1013 REC = REC*HALF
1014 IF( NOUNIT ) THEN
1015* TJJS = CONJG( A( J, J ) )*TSCAL
1016 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
1017 $ MYROW, MYCOL, IROW, ICOL, ITMP1,
1018 $ ITMP2 )
1019.EQ..AND..EQ. IF( ( MYROWITMP1 ) ( MYCOLITMP2 ) )
1020 $ THEN
1021 TJJS = CONJG( A( ( ICOL-1 )*LDA+IROW ) )*TSCAL
1022 CALL CGEBS2D( CONTXT, 'all', ' ', 1, 1, TJJS,
1023 $ 1 )
1024 ELSE
1025 CALL CGEBR2D( CONTXT, 'all', ' ', 1, 1, TJJS, 1,
1026 $ ITMP1, ITMP2 )
1027 END IF
1028 ELSE
1029 TJJS = TSCAL
1030 END IF
1031 TJJ = CABS1( TJJS )
1032.GT. IF( TJJONE ) THEN
1033*
1034* Divide by A(j,j) when scaling x if A(j,j) > 1.
1035*
1036 REC = MIN( ONE, REC*TJJ )
1037 USCAL = CLADIV( USCAL, TJJS )
1038 END IF
1039.LT. IF( RECONE ) THEN
1040 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
1041 XJTMP = XJTMP*REC
1042 SCALE = SCALE*REC
1043 XMAX( 1 ) = XMAX( 1 )*REC
1044 END IF
1045 END IF
1046*
1047 CSUMJ = CZERO
1048.EQ. IF( USCALCONE ) THEN
1049*
1050* If the scaling needed for A in the dot product is 1,
1051* call PCDOTC to perform the dot product.
1052*
1053 IF( UPPER ) THEN
1054 CALL PCDOTC( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1,
1055 $ X, IX, JX, DESCX, 1 )
1056.LT. ELSE IF( JN ) THEN
1057 CALL PCDOTC( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1,
1058 $ X, IX+J, JX, DESCX, 1 )
1059 END IF
1060.EQ. IF( MYCOLITMP2X ) THEN
1061 CALL CGEBS2D( CONTXT, 'row', ' ', 1, 1, CSUMJ, 1 )
1062 ELSE
1063 CALL CGEBR2D( CONTXT, 'row', ' ', 1, 1, CSUMJ, 1,
1064 $ MYROW, ITMP2X )
1065 END IF
1066 ELSE
1067*
1068* Otherwise, scale column of A by USCAL before dot
1069* product. Below is not the best way to do it.
1070*
1071 IF( UPPER ) THEN
1072* DO 180 I = 1, J - 1
1073* CSUMJ = CSUMJ + ( CONJG( A( I, J ) )*USCAL )*
1074* $ X( I )
1075* 180 CONTINUE
1076 ZDUM = CONJG( USCAL )
1077 CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 )
1078 CALL PCDOTC( J-1, CSUMJ, A, IA, JA+J-1, DESCA, 1,
1079 $ X, IX, JX, DESCX, 1 )
1080 ZDUM = CLADIV( CONE, ZDUM )
1081 CALL PCSCAL( J-1, ZDUM, A, IA, JA+J-1, DESCA, 1 )
1082.LT. ELSE IF( JN ) THEN
1083* DO 190 I = J + 1, N
1084* CSUMJ = CSUMJ + ( CONJG( A( I, J ) )*USCAL )*
1085* $ X( I )
1086* 190 CONTINUE
1087 ZDUM = CONJG( USCAL )
1088 CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 )
1089 CALL PCDOTC( N-J, CSUMJ, A, IA+J, JA+J-1, DESCA, 1,
1090 $ X, IX+J, JX, DESCX, 1 )
1091 ZDUM = CLADIV( CONE, ZDUM )
1092 CALL PCSCAL( N-J, ZDUM, A, IA+J, JA+J-1, DESCA, 1 )
1093 END IF
1094.EQ. IF( MYCOLITMP2X ) THEN
1095 CALL CGEBS2D( CONTXT, 'row', ' ', 1, 1, CSUMJ, 1 )
1096 ELSE
1097 CALL CGEBR2D( CONTXT, 'row', ' ', 1, 1, CSUMJ, 1,
1098 $ MYROW, ITMP2X )
1099 END IF
1100 END IF
1101*
1102.EQ. IF( USCALCMPLX( TSCAL ) ) THEN
1103*
1104* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
1105* was not used to scale the dotproduct.
1106*
1107* X( J ) = X( J ) - CSUMJ
1108* XJ = CABS1( X( J ) )
1109 XJTMP = XJTMP - CSUMJ
1110 XJ = CABS1( XJTMP )
1111* IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
1112* $ X( IROWX ) = XJTMP
1113 IF( NOUNIT ) THEN
1114* TJJS = CONJG( A( J, J ) )*TSCAL
1115 CALL INFOG2L( IA+J-1, JA+J-1, DESCA, NPROW, NPCOL,
1116 $ MYROW, MYCOL, IROW, ICOL, ITMP1,
1117 $ ITMP2 )
1118.EQ..AND..EQ. IF( ( MYROWITMP1 ) ( MYCOLITMP2 ) )
1119 $ THEN
1120 TJJS = CONJG( A( ( ICOL-1 )*LDA+IROW ) )*TSCAL
1121 CALL CGEBS2D( CONTXT, 'all', ' ', 1, 1, TJJS,
1122 $ 1 )
1123 ELSE
1124 CALL CGEBR2D( CONTXT, 'all', ' ', 1, 1, TJJS, 1,
1125 $ ITMP1, ITMP2 )
1126 END IF
1127 ELSE
1128 TJJS = TSCAL
1129.EQ. IF( TSCALONE )
1130 $ GO TO 130
1131 END IF
1132*
1133* Compute x(j) = x(j) / A(j,j), scaling if necessary.
1134*
1135 TJJ = CABS1( TJJS )
1136.GT. IF( TJJSMLNUM ) THEN
1137*
1138* abs(A(j,j)) > SMLNUM:
1139*
1140.LT. IF( TJJONE ) THEN
1141.GT. IF( XJTJJ*BIGNUM ) THEN
1142*
1143* Scale X by 1/abs(x(j)).
1144*
1145 REC = ONE / XJ
1146 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
1147 XJTMP = XJTMP*REC
1148 SCALE = SCALE*REC
1149 XMAX( 1 ) = XMAX( 1 )*REC
1150 END IF
1151 END IF
1152* X( J ) = CLADIV( X( J ), TJJS )
1153 XJTMP = CLADIV( XJTMP, TJJS )
1154.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
1155 $ X( IROWX ) = XJTMP
1156.GT. ELSE IF( TJJZERO ) THEN
1157*
1158* 0 < abs(A(j,j)) <= SMLNUM:
1159*
1160.GT. IF( XJTJJ*BIGNUM ) THEN
1161*
1162* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
1163*
1164 REC = ( TJJ*BIGNUM ) / XJ
1165 CALL PCSSCAL( N, REC, X, IX, JX, DESCX, 1 )
1166 XJTMP = XJTMP*REC
1167 SCALE = SCALE*REC
1168 XMAX( 1 ) = XMAX( 1 )*REC
1169 END IF
1170* X( J ) = CLADIV( X( J ), TJJS )
1171 XJTMP = CLADIV( XJTMP, TJJS )
1172.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
1173 $ X( IROWX ) = XJTMP
1174 ELSE
1175*
1176* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
1177* scale = 0 and compute a solution to A**H *x = 0.
1178*
1179 CALL PCLASET( ' ', N, 1, CZERO, CZERO, X, IX, JX,
1180 $ DESCX )
1181.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
1182 $ X( IROWX ) = CONE
1183 XJTMP = CONE
1184 SCALE = ZERO
1185 XMAX( 1 ) = ZERO
1186 END IF
1187 130 CONTINUE
1188 ELSE
1189*
1190* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
1191* product has already been divided by 1/A(j,j).
1192*
1193* X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
1194 XJTMP = CLADIV( XJTMP, TJJS ) - CSUMJ
1195.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) )
1196 $ X( IROWX ) = XJTMP
1197 END IF
1198 XMAX( 1 ) = MAX( XMAX( 1 ), CABS1( XJTMP ) )
1199 140 CONTINUE
1200 END IF
1201 SCALE = SCALE / TSCAL
1202 END IF
1203*
1204* Scale the column norms by 1/TSCAL for return.
1205*
1206.NE. IF( TSCALONE ) THEN
1207 CALL SSCAL( N, ONE / TSCAL, CNORM, 1 )
1208 END IF
1209*
1210 RETURN
1211*
1212* End of PCLATTRS
1213*
1214 END
cmplx
float cmplx[2]
Definition pblas.h:136
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
cgebs2d
subroutine cgebs2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1062
pcscal
subroutine pcscal(n, alpha, x, ix, jx, descx, incx)
Definition mpi.f:955
pxerbla
subroutine pxerbla(contxt, srname, info)
Definition mpi.f:1600
cgebr2d
subroutine cgebr2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1103
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition mpi.f:937
pcaxpy
subroutine pcaxpy(n, a, x, ix, jx, descx, incx, y, iy, jy, descy, incy)
Definition mpi.f:1425
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
pclaset
subroutine pclaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pcblastst.f:7508
pclattrs
subroutine pclattrs(uplo, trans, diag, normin, n, a, ia, ja, desca, x, ix, jx, descx, scale, cnorm, info)
Definition pclattrs.f:3
pslabad
subroutine pslabad(ictxt, small, large)
Definition pslabad.f:2