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.NOT..AND..NOT. IF( UPPER LSAME( UPLO, 'l' ) ) THEN
320 INFO = -1
321.NOT..AND..NOT. ELSE IF( NOTRAN LSAME( TRANS, 't.AND..NOT.' )
322 $ LSAME( TRANS, 'c' ) ) THEN
323 INFO = -2
324.NOT..AND..NOT. ELSE IF( NOUNIT LSAME( DIAG, 'u' ) ) THEN
325 INFO = -3
326.NOT. ELSE IF( LSAME( NORMIN, 'y.AND..NOT.' )
327 $ LSAME( NORMIN, 'n' ) ) THEN
328 INFO = -4
329.LT. ELSE IF( N0 ) THEN
330 INFO = -5
331 END IF
332*
333 CALL BLACS_GRIDINFO( CONTXT, NPROW, NPCOL, MYROW, MYCOL )
334*
335.NE. IF( INFO0 ) THEN
336 CALL PXERBLA( CONTXT, 'pclattrs', -INFO )
337 RETURN
338 END IF
339*
340* Quick return if possible
341*
342.EQ. IF( N0 )
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.LE. IF( TMAXBIGNUM*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.NE. IF( TSCALONE ) 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.LE. IF( GROWSMLNUM )
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.EQ..AND..EQ. IF( ( MYROWITMP1 ) ( MYCOLITMP2 ) ) 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.GE. IF( TJJSMLNUM ) 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.GE. IF( TJJ+CNORM( J )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.LE. IF( GROWSMLNUM )
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.NE. IF( TSCALONE ) 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.LE. IF( GROWSMLNUM )
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.EQ..AND..EQ. IF( ( MYROWITMP1 ) ( MYCOLITMP2 ) ) 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.GE. IF( TJJSMLNUM ) THEN
547*
548* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
549*
550.GT. IF( XJTJJ )
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.LE. IF( GROWSMLNUM )
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.GT. IF( ( GROW*TSCAL )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.GT. IF( XMAX( 1 )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.EQ..AND..EQ. IF( ( MYROWITMP1X ) ( MYCOLITMP2X ) ) 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 IF( j.LT.n ) 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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) ) 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 IF( cnorm( j ).GT.( 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 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
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 IF( tjj.GT.one ) 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 IF( rec.LT.one ) 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 IF( uscal.EQ.cone ) 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 ELSE IF( j.LT.n ) 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 IF( mycol.EQ.itmp2x ) 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 ELSE IF( j.LT.n ) 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 IF( mycol.EQ.itmp2x ) 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 IF( uscal.EQ.cmplx( 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 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
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 IF( tscal.EQ.one )
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 IF( tjj.GT.smlnum ) THEN
916*
917* abs(A(j,j)) > SMLNUM:
918*
919 IF( tjj.LT.one ) THEN
920 IF( xj.GT.tjj*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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
934 $ THEN
935 x( irowx ) = xjtmp
936 END IF
937 ELSE IF( tjj.GT.zero ) THEN
938*
939* 0 < abs(A(j,j)) <= SMLNUM:
940*
941 IF( xj.GT.tjj*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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) ) 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 IF( cnorm( j ).GT.( 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 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
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 IF( tjj.GT.one ) 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 IF( rec.LT.one ) 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 IF( uscal.EQ.cone ) 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 ELSE IF( j.LT.n ) 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 IF( mycol.EQ.itmp2x ) 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 ELSE IF( j.LT.n ) 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 IF( mycol.EQ.itmp2x ) 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 IF( uscal.EQ.cmplx( 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 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) )
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 IF( tscal.EQ.one )
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 IF( tjj.GT.smlnum ) THEN
1137*
1138* abs(A(j,j)) > SMLNUM:
1139*
1140 IF( tjj.LT.one ) THEN
1141 IF( xj.GT.tjj*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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
1155 $ x( irowx ) = xjtmp
1156 ELSE IF( tjj.GT.zero ) THEN
1157*
1158* 0 < abs(A(j,j)) <= SMLNUM:
1159*
1160 IF( xj.GT.tjj*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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
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 IF( ( myrow.EQ.itmp1x ) .AND. ( mycol.EQ.itmp2x ) )
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 IF( tscal.NE.one ) 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