Routine: PCLATTRS()  File: SRC\pclattrs.f

 
 
# lines: 1213
  # code: 1213
  # comment: 0
  # blank:0
# Variables:43
# Callers:1
# Callings:1
# Words:426
# Keywords:192
 

 

 
..
     .. Local Scalars ..
     ..
     .. External Functions ..
     ..
     .. External Subroutines ..
     ..
     .. Intrinsic Functions ..
     ..
     .. Statement Functions ..
     ..
     .. Statement Function definitions ..
     ..
     .. Executable Statements ..
     Test the input parameters.
     Quick return if possible
     Determine machine dependent parameters to control overflow.
        Compute the 1-norm of each column, not including the diagonal.
           A is upper triangular.
           A is lower triangular.
     Scale the column norms by TSCAL if the maximum element in CNORM is
     greater than BIGNUM/2.
     Compute a bound on the computed solution vector to see if the
     Level 2 PBLAS routine PCTRSV can be used.
        Compute the growth in A * x = b.
           A is non-unit triangular.
           Compute GROW = 1/G(j) and XBND = 1/M(j).
           Initially, G(0) = max{x(i), i=1,...,n}.
              Exit the loop if the growth factor is too small.
              TJJS = A( J, J )
                 M(j) = G(j-1) / abs(A(j,j))
                 M(j) could overflow, set XBND to 0.
                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
                 G(j) could overflow, set GROW to 0.
           A is unit triangular.
           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
              Exit the loop if the growth factor is too small.
              G(j) = G(j-1)*( 1 + CNORM(j) )
        Compute the growth in A**T * x = b  or  A**H * x = b.
           A is non-unit triangular.
           Compute GROW = 1/G(j) and XBND = 1/M(j).
           Initially, M(0) = max{x(i), i=1,...,n}.
              Exit the loop if the growth factor is too small.
              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
              TJJS = A( J, J )
                 M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
                 M(j) could overflow, set XBND to 0.
           A is unit triangular.
           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
              Exit the loop if the growth factor is too small.
              G(j) = ( 1 + CNORM(j) )*G(j-1)
        Use the Level 2 PBLAS solve if the reciprocal of the bound on
        elements of X is not too small.
        Use a Level 1 PBLAS solve, scaling intermediate results.
           Scale X so that its components are less than or equal to
           BIGNUM in absolute value.
           Solve A * x = b
              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
              XJ = CABS1( X( J ) )
                 TJJS = A( J, J )*TSCAL
                    abs(A(j,j)) > SMLNUM:
                          Scale x by 1/b(j).
                 X( J ) = CLADIV( X( J ), TJJS )
                 XJ = CABS1( X( J ) )
                    0 < abs(A(j,j)) <= SMLNUM:
                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
                       to avoid overflow when dividing by A(j,j).
                          Scale by 1/CNORM(j) to avoid overflow when
                          multiplying x(j) times column j.
                 X( J ) = CLADIV( X( J ), TJJS )
                 XJ = CABS1( X( J ) )
                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
                    scale = 0, and compute a solution to A*x = 0.
              Scale x if necessary to avoid overflow when adding a
              multiple of column j of A.
                    Scale x by 1/(2*abs(x(j))).
                 Scale x by 1/2.
                    Compute the update
                       x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
                    Compute the update
                       x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
           Solve A**T * x = b
              Compute x(j) = b(j) - sum A(k,j)*x(k).
                                    k<>j
              XJ = CABS1( X( J ) )
                 If x(j) could overflow, scale x by 1/(2*XMAX).
                    TJJS = A( J, J )*TSCAL
                       Divide by A(j,j) when scaling x if A(j,j) > 1.
                 If the scaling needed for A in the dot product is 1,
                 call PCDOTU to perform the dot product.
                 Otherwise, scale column of A by USCAL before dot
                 product.  Below is not the best way to do it.
                    DO 130 I = 1, J - 1
                       CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
 130                CONTINUE
                    DO 140 I = J + 1, N
                       CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
  140               CONTINUE
                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
                 was not used to scale the dotproduct.
                 X( J ) = X( J ) - CSUMJ
                 XJ = CABS1( X( J ) )
                  IF( ( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) )
     $               X( IROWX ) = XJTMP
                    TJJS = A( J, J )*TSCAL
                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
                       abs(A(j,j)) > SMLNUM:
                             Scale X by 1/abs(x(j)).
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PCLATTRS( UPLO , TRANS , DIAG , NORMIN , N , A , IA , JA ,
002       $DESCA , X , IX , JX , DESCX , SCALE , CNORM , INFO )
003  
004  *     -- ScaLAPACK routine(version 1.7) --
005  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
006  *     and University of California , Berkeley.
007  *     July 31 , 2001
008  
009  *     .. Scalar Arguments ..
010        CHARACTER DIAG , NORMIN , TRANS , UPLO
011        INTEGER IA , INFO , IX , JA , JX , N
012        REAL SCALE
013  *     ..
014  *     .. Array Arguments ..
015        INTEGER DESCA( * ) , DESCX( * )
016        REAL CNORM( * )
017        COMPLEX A( * ) , X( * )
018  *     ..
019  
020  *     Purpose
021  *     === ====
022  
023  *     PCLATTRS solves one of the triangular systems
024  
025  *     A * x = s*b , A**T * x = s*b , or A**H * x = s*b ,
026  
027  *     with scaling to prevent overflow. Here A is an upper or lower
028  *     triangular matrix , A**T denotes the transpose of A , A**H denotes the
029  *     conjugate transpose of A , x and b are n - element vectors , and s is a
030  *     scaling factor , usually less than or equal to 1 , chosen so that the
031  *     components of x will be less than the overflow threshold. If the
032  *     unscaled problem will not cause overflow , the Level 2 PBLAS routine
033  *     PCTRSV is called. If the matrix A is singular(A(j , j) = 0 for some j)
034  *     then s is set to 0 and a non - trivial solution to A*x = 0 is returned.
035  
036  *     This is very slow relative to PCTRSV. This should only be used
037  *     when scaling is necessary to control overflow , or when it is modified
038  *     to scale better.
039  *     Notes
040  
041  *     === ==
042  
043  *     Each global data object is described by an associated description
044  *     vector. This vector stores the information required to establish
045  *     the mapping between an object element and its corresponding process
046  *     and memory location.
047  
048  *     Let A be a generic term for any 2D block cyclicly distributed array.
049  *     Such a global array has an associated description vector DESCA.
050  *     In the following comments , the character _ should be read as
051  *     "of the global array".
052  
053  *     NOTATION STORED IN EXPLANATION
054  *     --- ------------ -------------- --------------------------------------
055  *     DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case ,
056  *     DTYPE_A = 1.
057  *     CTXT_A(global) DESCA( CTXT_ ) The BLACS context handle , indicating
058  *     the BLACS process grid A is distribu -
059  *     ted over. The context itself is glo -
060  *     bal , but the handle(the integer
061  *     value) may vary.
062  *     M_A(global) DESCA( M_ ) The number of rows in the global
063  *     array A.
064  *     N_A(global) DESCA( N_ ) The number of columns in the global
065  *     array A.
066  *     MB_A(global) DESCA( MB_ ) The blocking factor used to distribute
067  *     the rows of the array.
068  *     NB_A(global) DESCA( NB_ ) The blocking factor used to distribute
069  *     the columns of the array.
070  *     RSRC_A(global) DESCA( RSRC_ ) The process row over which the first
071  *     row of the array A is distributed.
072  *     CSRC_A(global) DESCA( CSRC_ ) The process column over which the
073  *     first column of the array A is
074  *     distributed.
075  *     LLD_A(local) DESCA( LLD_ ) The leading dimension of the local
076  *     array. LLD_A >= MAX(1 , LOCr(M_A)).
077  
078  *     Let K be the number of rows or columns of a distributed matrix ,
079  *     and assume that its process grid has dimension r x c.
080  *     LOCr( K ) denotes the number of elements of K that a process
081  *     would receive if K were distributed over the r processes of its
082  *     process column.
083  *     Similarly , LOCc( K ) denotes the number of elements of K that a
084  *     process would receive if K were distributed over the c processes of
085  *     its process row.
086  *     The values of LOCr() and LOCc() may be determined via a call to the
087  *     ScaLAPACK tool function , NUMROC :
088  *     LOCr( M ) = NUMROC( M , MB_A , MYROW , RSRC_A , NPROW ) ,
089  *     LOCc( N ) = NUMROC( N , NB_A , MYCOL , CSRC_A , NPCOL ).
090  *     An upper bound for these quantities may be computed by :
091  *     LOCr( M ) <= ceil( ceil(M / MB_A) / NPROW )*MB_A
092  *     LOCc( N ) <= ceil( ceil(N / NB_A) / NPCOL )*NB_A
093  
094  *     Arguments
095  *     === ======
096  
097  *     UPLO(global input) CHARACTER*1
098  *     Specifies whether the matrix A is upper or lower triangular.
099  *     = '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        REC = ONE / XJ
267        CALL PCSSCAL( N , REC , X , IX , JX , DESCX , 1 )
268        XJTMP = XJTMP*REC
269        SCALE = SCALE*REC
270        XMAX = XMAX*REC
271        END IF
272        END IF
273  *     X( J ) = CLADIV( X( J ) , TJJS )
274        XJTMP = CLADIV( XJTMP , TJJS )
275        IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
276       $    THEN
277            X( IROWX ) = XJTMP
278        END IF
279        ELSE IF( TJJ.GT.ZERO ) THEN
280  
281  *         0 < abs(A(j , j)) <= SMLNUM :
282  
283            IF( XJ.GT.TJJ*BIGNUM ) THEN
284  
285  *             Scale x by(1 / abs(x(j)))*abs(A(j , j))*BIGNUM.
286  
287                REC =( TJJ*BIGNUM ) / XJ
288                CALL PCSSCAL( N , REC , X , IX , JX , DESCX , 1 )
289                XJTMP = XJTMP*REC
290                SCALE = SCALE*REC
291                XMAX = XMAX*REC
292            END IF
293  *         X( J ) = CLADIV( X( J ) , TJJS )
294            XJTMP = CLADIV( XJTMP , TJJS )
295            IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
296       $        THEN
297                X( IROWX ) = XJTMP
298            END IF
299        ELSE
300  
301  *         A(j , j) = 0 : Set x(1 : n) = 0 , x(j) = 1 , and
302  *         scale = 0 and compute a solution to A**T *x = 0.
303  
304            CALL PCLASET ( ' ' , N , 1 , CZERO , CZERO , X , IX , JX ,
305       $    DESCX )
306            IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
307       $        THEN
308                X( IROWX ) = CONE
309            END IF
310            XJTMP = CONE
311            SCALE = ZERO
312            XMAX = ZERO
313        END IF
314    110 CONTINUE
315        ELSE
316  
317  *         Compute x(j) := x(j) / A(j , j) - CSUMJ if the dot
318  *         product has already been divided by 1 / A(j , j).
319  
320  *         X( J ) = CLADIV( X( J ) , TJJS ) - CSUMJ
321            XJTMP = CLADIV( XJTMP , TJJS ) - CSUMJ
322            IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
323       $        THEN
324                X( IROWX ) = XJTMP
325            END IF
326        END IF
327        XMAX = MAX( XMAX , CABS1( XJTMP ) )
328    120 CONTINUE
329  
330        ELSE
331  
332  *         Solve A**H * x = b
333  
334            DO 140 J = JFIRST , JLAST , JINC
335  
336  *             Compute x(j) = b(j) - sum A(k , j)*x(k).
337  *             k <> j
338  
339                CALL INFOG2L( IX + J - 1 , JX , DESCX , NPROW , NPCOL , MYROW ,
340       $        MYCOL , IROWX , ICOLX , ITMP1X , ITMP2X )
341                IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) ) THEN
342                    XJTMP = X( IROWX )
343                    CALL CGEBS2D( CONTXT , 'All' , ' ' , 1 , 1 , XJTMP , 1 )
344                ELSE
345                    CALL CGEBR2D( CONTXT , 'All' , ' ' , 1 , 1 , XJTMP , 1 ,
346       $            ITMP1X , ITMP2X )
347                END IF
348                XJ = CABS1( XJTMP )
349                USCAL = TSCAL
350                REC = ONE / MAX( XMAX , ONE )
351                IF( CNORM( J ).GT.( BIGNUM - XJ )*REC ) THEN
352  
353  *                 If x(j) could overflow , scale x by 1 / (2*XMAX).
354  
355                    REC = REC*HALF
356                    IF( NOUNIT ) THEN
357  *                     TJJS = CONJG( A( J , J ) )*TSCAL
358                        CALL INFOG2L( IA + J - 1 , JA + J - 1 , DESCA , NPROW , NPCOL ,
359       $                MYROW , MYCOL , IROW , ICOL , ITMP1 ,
360       $                ITMP2 )
361                        IF(( MYROW.EQ.ITMP1 ) .AND.( MYCOL.EQ.ITMP2 ) )
362       $                    THEN
363                            TJJS = CONJG( A(( ICOL - 1 )*LDA + IROW ) )*TSCAL
364                            CALL CGEBS2D( CONTXT , 'All' , ' ' , 1 , 1 , TJJS ,
365       $                    1 )
366                        ELSE
367                            CALL CGEBR2D( CONTXT , 'All' , ' ' , 1 , 1 , TJJS , 1 ,
368       $                    ITMP1 , ITMP2 )
369                        END IF
370                    ELSE
371                        TJJS = TSCAL
372                    END IF
373                    TJJ = CABS1( TJJS )
374                    IF( TJJ.GT.ONE ) THEN
375  
376  *                     Divide by A(j , j) when scaling x if A(j , j) > 1.
377  
378                        REC = MIN( ONE , REC*TJJ )
379                        USCAL = CLADIV( USCAL , TJJS )
380                    END IF
381                    IF( REC.LT.ONE ) THEN
382                        CALL PCSSCAL( N , REC , X , IX , JX , DESCX , 1 )
383                        XJTMP = XJTMP*REC
384                        SCALE = SCALE*REC
385                        XMAX = XMAX*REC
386                    END IF
387                END IF
388  
389                CSUMJ = CZERO
390                IF( USCAL.EQ.CONE ) THEN
391  
392  *                 If the scaling needed for A in the dot product is 1 ,
393  *                 call PCDOTC to perform the dot product.
394  
395                    IF( UPPER ) THEN
396                        CALL PCDOTC( J - 1 , CSUMJ , A , IA , JA + J - 1 , DESCA , 1 ,
397       $                X , IX , JX , DESCX , 1 )
398                    ELSE IF( J.LT.N ) THEN
399                        CALL PCDOTC( N - J , CSUMJ , A , IA + J , JA + J - 1 , DESCA , 1 ,
400       $                X , IX + J , JX , DESCX , 1 )
401                    END IF
402                    IF( MYCOL.EQ.ITMP2X ) THEN
403                        CALL CGEBS2D( CONTXT , 'Row' , ' ' , 1 , 1 , CSUMJ , 1 )
404                    ELSE
405                        CALL CGEBR2D( CONTXT , 'Row' , ' ' , 1 , 1 , CSUMJ , 1 ,
406       $                MYROW , ITMP2X )
407                    END IF
408                ELSE
409  
410  *                 Otherwise , scale column of A by USCAL before dot
411  *                 product. Below is not the best way to do it.
412  
413                    IF( UPPER ) THEN
414  *                     DO 180 I = 1 , J - 1
415  *                     CSUMJ = CSUMJ + ( CONJG( A( I , J ) )*USCAL )*
416  *                     $                          X( I )
417  *                     180                CONTINUE
418                        ZDUM = CONJG( USCAL )
419                        CALL PCSCAL( J - 1 , ZDUM , A , IA , JA + J - 1 , DESCA , 1 )
420                        CALL PCDOTC( J - 1 , CSUMJ , A , IA , JA + J - 1 , DESCA , 1 ,
421       $                X , IX , JX , DESCX , 1 )
422                        ZDUM = CLADIV( CONE , ZDUM )
423                        CALL PCSCAL( J - 1 , ZDUM , A , IA , JA + J - 1 , DESCA , 1 )
424                    ELSE IF( J.LT.N ) THEN
425  *                     DO 190 I = J + 1 , N
426  *                     CSUMJ = CSUMJ + ( CONJG( A( I , J ) )*USCAL )*
427  *                     $                          X( I )
428  *                     190                CONTINUE
429                        ZDUM = CONJG( USCAL )
430                        CALL PCSCAL( N - J , ZDUM , A , IA + J , JA + J - 1 , DESCA , 1 )
431                        CALL PCDOTC( N - J , CSUMJ , A , IA + J , JA + J - 1 , DESCA , 1 ,
432       $                X , IX + J , JX , DESCX , 1 )
433                        ZDUM = CLADIV( CONE , ZDUM )
434                        CALL PCSCAL( N - J , ZDUM , A , IA + J , JA + J - 1 , DESCA , 1 )
435                    END IF
436                    IF( MYCOL.EQ.ITMP2X ) THEN
437                        CALL CGEBS2D( CONTXT , 'Row' , ' ' , 1 , 1 , CSUMJ , 1 )
438                    ELSE
439                        CALL CGEBR2D( CONTXT , 'Row' , ' ' , 1 , 1 , CSUMJ , 1 ,
440       $                MYROW , ITMP2X )
441                    END IF
442                END IF
443  
444                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
445  
446  *                 Compute x(j) :=( x(j) - CSUMJ ) / A(j , j) if 1 / A(j , j)
447  *                 was not used to scale the dotproduct.
448  
449  *                 X( J ) = X( J ) - CSUMJ
450  *                 XJ = CABS1( X( J ) )
451                    XJTMP = XJTMP - CSUMJ
452                    XJ = CABS1( XJTMP )
453  *                 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
454  *                 $ X( IROWX ) = XJTMP
455                    IF( NOUNIT ) THEN
456  *                     TJJS = CONJG( A( J , J ) )*TSCAL
457                        CALL INFOG2L( IA + J - 1 , JA + J - 1 , DESCA , NPROW , NPCOL ,
458       $                MYROW , MYCOL , IROW , ICOL , ITMP1 ,
459       $                ITMP2 )
460                        IF(( MYROW.EQ.ITMP1 ) .AND.( MYCOL.EQ.ITMP2 ) )
461       $                    THEN
462                            TJJS = CONJG( A(( ICOL - 1 )*LDA + IROW ) )*TSCAL
463                            CALL CGEBS2D( CONTXT , 'All' , ' ' , 1 , 1 , TJJS ,
464       $                    1 )
465                        ELSE
466                            CALL CGEBR2D( CONTXT , 'All' , ' ' , 1 , 1 , TJJS , 1 ,
467       $                    ITMP1 , ITMP2 )
468                        END IF
469                    ELSE
470                        TJJS = TSCAL
471                        IF( TSCAL.EQ.ONE )
472       $                    GO TO 130
473                        END IF
474  
475  *                     Compute x(j) = x(j) / A(j , j) , scaling if necessary.
476  
477                        TJJ = CABS1( TJJS )
478                        IF( TJJ.GT.SMLNUM ) THEN
479  
480  *                         abs(A(j , j)) > SMLNUM :
481  
482                            IF( TJJ.LT.ONE ) THEN
483                                IF( XJ.GT.TJJ*BIGNUM ) THEN
484  
485  *                                 Scale X by 1 / abs(x(j)).
486  
487                                    REC = ONE / XJ
488                                    CALL PCSSCAL( N , REC , X , IX , JX , DESCX , 1 )
489                                    XJTMP = XJTMP*REC
490                                    SCALE = SCALE*REC
491                                    XMAX = XMAX*REC
492                                END IF
493                            END IF
494  *                         X( J ) = CLADIV( X( J ) , TJJS )
495                            XJTMP = CLADIV( XJTMP , TJJS )
496                            IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
497       $                        X( IROWX ) = XJTMP
498                            ELSE IF( TJJ.GT.ZERO ) THEN
499  
500  *                             0 < abs(A(j , j)) <= SMLNUM :
501  
502                                IF( XJ.GT.TJJ*BIGNUM ) THEN
503  
504  *                                 Scale x by(1 / abs(x(j)))*abs(A(j , j))*BIGNUM.
505  
506                                    REC =( TJJ*BIGNUM ) / XJ
507                                    CALL PCSSCAL( N , REC , X , IX , JX , DESCX , 1 )
508                                    XJTMP = XJTMP*REC
509                                    SCALE = SCALE*REC
510                                    XMAX = XMAX*REC
511                                END IF
512  *                             X( J ) = CLADIV( X( J ) , TJJS )
513                                XJTMP = CLADIV( XJTMP , TJJS )
514                                IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
515       $                            X( IROWX ) = XJTMP
516                                ELSE
517  
518  *                                 A(j , j) = 0 : Set x(1 : n) = 0 , x(j) = 1 , and
519  *                                 scale = 0 and compute a solution to A**H *x = 0.
520  
521                                    CALL PCLASET ( ' ' , N , 1 , CZERO , CZERO , X , IX , JX ,
522       $                            DESCX )
523                                    IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
524       $                                X( IROWX ) = CONE
525                                        XJTMP = CONE
526                                        SCALE = ZERO
527                                        XMAX = ZERO
528                                    END IF
529    130 CONTINUE
530        ELSE
531  
532  *         Compute x(j) := x(j) / A(j , j) - CSUMJ if the dot
533  *         product has already been divided by 1 / A(j , j).
534  
535  *         X( J ) = CLADIV( X( J ) , TJJS ) - CSUMJ
536            XJTMP = CLADIV( XJTMP , TJJS ) - CSUMJ
537            IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
538       $        X( IROWX ) = XJTMP
539            END IF
540            XMAX = MAX( XMAX , CABS1( XJTMP ) )
541    140 CONTINUE
542        END IF
543        SCALE = SCALE / TSCAL
544        END IF
545  
546  *     Scale the column norms by 1 / TSCAL for return.
547  
548        IF( TSCAL.NE.ONE ) THEN
549            CALL SSCAL( N , ONE / TSCAL , CNORM , 1 )
550        END IF
551  
552        RETURN
553  
554  *     End of PCLATTRS
555  
556        END