|
|
| |
| # lines: |
1213 | | # code: |
1213 | | # comment: | 0 | |
# blank: | 0 |
| # Variables: | 43 |
| # Callers: | 1 |
| # Callings: | 1 |
| # Words: | 429 |
| # 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 PZTRSV 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 ) = ZLADIV( 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 ) = ZLADIV( 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 PZDOTU 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)).
|
|
|
|
001 SUBROUTINE PZLATTRS( 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 DOUBLE PRECISION SCALE
013 * ..
014 * .. Array Arguments ..
015 INTEGER DESCA( * ) , DESCX( * )
016 DOUBLE PRECISION CNORM( * )
017 COMPLEX*16 A( * ) , X( * )
018 * ..
019
020 * Purpose
021 * === ====
022
023 * PZLATTRS 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 * PZTRSV 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 PZTRSV. 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*16 array , dimension(DESCA(LLD_) ,*)
123 * The triangular matrix A. If UPLO = 'U' , the leading n by n
124 * upper triangular part of the array A contains the upper
125 * triangular matrix , and the strictly lower triangular part of
126 * A is not referenced. If UPLO = 'L' , the leading n by n lower
127 * triangular part of the array A contains the lower triangular
128 * matrix , and the strictly upper triangular part of A is not
129 * referenced. If DIAG = 'U' , the diagonal elements of A are
130 * also not referenced and are assumed to be 1.
131
132 * IA(global input) pointer to INTEGER
133 * The global row index of the submatrix of the distributed
134 * matrix A to operate on.
135
136 * JA(global input) pointer to INTEGER
137 * The global column index of the submatrix of the distributed
138 * matrix A to operate on.
139
140 * DESCA(global and local input) INTEGER array of dimension DLEN_.
141 * The array descriptor for the distributed matrix A.
142
143 * X(local input / output) COMPLEX*16 array ,
144 * dimension(DESCX(LLD_) ,*)
145 * On entry , the right hand side b of the triangular system.
146 * On exit , X is overwritten by the solution vector x.
147
148 * IX(global input) pointer to INTEGER
149 * The global row index of the submatrix of the distributed
150 * matrix X to operate on.
151
152 * JX(global input) pointer to INTEGER
153 * The global column index of the submatrix of the distributed
154 * matrix X to operate on.
155
156 * DESCX(global and local input) INTEGER array of dimension DLEN_.
157 * The array descriptor for the distributed matrix X.
158
159 * SCALE(global output) DOUBLE PRECISION
160 * The scaling factor s for the triangular system
161 * A * x = s*b , A**T * x = s*b , or A**H * x = s*b.
162 * If SCALE = 0 , the matrix A is singular or badly scaled , and
163 * the vector x is an exact or approximate solution to A*x = 0.
164
165 * CNORM(global input or global output) DOUBLE PRECISION array ,
166 * dimension(N)
167 * If NORMIN = 'Y' , CNORM is an input argument and CNORM(j)
168 * contains the norm of the off - diagonal part of the j - th column
169 * of A. If TRANS = 'N' , CNORM(j) must be greater than or equal
170 * to the infinity - norm , and if TRANS = 'T' or 'C' , CNORM(j)
171 * must be greater than or equal to the 1 - norm.
172
173 * If NORMIN = 'N' , CNORM is an output argument and CNORM(j)
174 * returns the 1 - norm of the offdiagonal part of the j - th column
175 * of A.
176
177 * INFO(global output) INTEGER
178 * = 0 : successful exit
179 * < 0 : if INFO = - k , the k - th argument had an illegal value
180
181 * Further Details
182 * === ==== =======
183
184 * A rough bound on x is computed ; if that is less than overflow , PZTRSV
185 * is called , otherwise , specific code is used which checks for possible
186 * overflow or divide - by - zero at every operation.
187
188 * A columnwise scheme is used for solving A*x = b. The basic algorithm
189 * if A is lower triangular is
190
191 * x[1 : n] := b[1 : n]
192 * for j = 1 , ... , n
193 * x(j) := x(j) / A(j , j)
194 * x[j + 1 : n] := x[j + 1 : n] - x(j) * A[j + 1 : n , j]
195 * end
196
197 * Define bounds on the components of x after j iterations of the loop :
198 * M(j) = bound on x[1 : j]
199 * G(j) = bound on x[j + 1 : n]
200 * Initially , let M(0) = 0 and G(0) = max{x(i) , i = 1 , ... , n}.
201
202 * Then for iteration j + 1 we have
203 * M(j + 1) <= G(j) / | A(j + 1 , j + 1) |
204 * G(j + 1) <= G(j) + M(j + 1) * | A[j + 2 : n , j + 1] |
205 * <= G(j)( 1 + CNORM(j + 1) / | A(j + 1 , j + 1) | )
206
207 * where CNORM(j + 1) is greater than or equal to the infinity - norm of
208 * column j + 1 of A , not counting the diagonal. Hence
209
210 * G(j) <= G(0) product( 1 + CNORM(i) / | A(i , i) | )
211 * 1 <= i <= j
212 * and
213
214 * |x(j)| <=( G(0) / |A(j , j)| ) product( 1 + CNORM(i) / |A(i , i)| )
215 * 1 <= i < j
216
217 * Since |x(j)| <= M(j) , we use the Level 2 PBLAS routine PZTRSV if the
218 * reciprocal of the largest M(j) , j = 1 , .. , n , is larger than
219 * max(underflow , 1 / overflow).
220
221 * The bound on x(j) is also used to determine when a step in the
222 * columnwise method can be performed without fear of overflow. If
223 * the computed bound is greater than a large constant , x is scaled to
224 * prevent overflow , but if the bound overflows , x is set to 0 , x(j) to
225 * 1, and scale to 0 , and a non - trivial solution to A*x = 0 is found.
226
227 * Similarly , a row - wise scheme is used to solve A**T *x = b or
228 * A**H *x = b. The basic algorithm for A upper triangular is
229
230 * for j = 1 , ... , n
231 * x(j) :=( b(j) - A[1 : j - 1 , j]' * x[1 : j - 1] ) / A(j , j)
232 * end
233
234 * We simultaneously compute two bounds
235 * G(j) = bound on( b(i) - A[1 : i - 1 , i]' * x[1 : i - 1] ) , 1 <= i <= j
236 * M(j) = bound on x(i) , 1 <= i <= j
237
238 * The initial values are G(0) = 0 , M(0) = max{b(i) , i = 1 , .. , n} , and we
239 * add the constraint G(j) >= G(j - 1) and M(j) >= M(j - 1) for j >= 1.
240 * Then the bound on x(j) is
241
242 * M(j) <= M(j - 1) * ( 1 + CNORM(j) ) / | A(j , j) |
243
244 * <= M(0) * product(( 1 + CNORM(i) ) / |A(i , i)| )
245 * 1 <= i <= j
246
247 * and we can safely call PZTRSV if 1 / M(n) and 1 / G(n) are both greater
248 * than max(underflow , 1 / overflow).
249
250 * Last modified by : Mark R. Fahey , August 2000
251
252 * === ==================================================================
253
254 * .. Parameters ..
255 DOUBLE PRECISION ZERO , HALF , ONE , TWO
256 PARAMETER( ZERO = 0.0D + 0 , HALF = 0.5D + 0 , ONE = 1.0D + 0 ,
257 $TWO = 2.0D + 0 )
258 COMPLEX*16 CZERO , CONE
259 PARAMETER( CZERO =( 0.0D + 0 , 0.0D + 0 ) ,
260 $CONE =( 1.0D + 0 , 0.0D + 0 ) )
261 INTEGER BLOCK_CYCLIC_2D , DLEN_ , DTYPE_ , CTXT_ , M_ , N_ ,
262 $MB_ , NB_ , RSRC_ , CSRC_ , LLD_
263 PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
264 $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
265 $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
266 REC = ONE / XJ
267 CALL PZDSCAL( 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 ) = ZLADIV( X( J ) , TJJS )
274 XJTMP = ZLADIV( XJTMP , TJJS )
275 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
275  
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
282
283 IF( XJ.GT.TJJ*BIGNUM ) THEN
284
285 * Scale x by(1 / abs(x(j)))*abs(A(j , j))*BIGNUM.
286
286
287 REC =( TJJ*BIGNUM ) / XJ
288 CALL PZDSCAL( 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 ) = ZLADIV( X( J ) , TJJS )
294 XJTMP = ZLADIV( XJTMP , TJJS )
295 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
295
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
303
304 CALL PZLASET ( ' ' , N , 1 , CZERO , CZERO , X , IX , JX ,
305 $ DESCX )
306 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
306
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 ) = ZLADIV( X( J ) , TJJS ) - CSUMJ
320
321 XJTMP = ZLADIV( XJTMP , TJJS ) - CSUMJ
322 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
322
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
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
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
341
342 XJTMP = X( IROWX )
343 CALL ZGEBS2D( CONTXT , 'All' , ' ' , 1 , 1 , XJTMP , 1 )
344 ELSE
344
345 CALL ZGEBR2D( 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
354
355 REC = REC*HALF
356 IF( NOUNIT ) THEN
357 * TJJS = DCONJG( A( J , J ) )*TSCAL
357
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 ) )
361
362 $ THEN
363 TJJS = DCONJG( A(( ICOL - 1 )*LDA + IROW ) )*TSCAL
364 CALL ZGEBS2D( CONTXT , 'All' , ' ' , 1 , 1 , TJJS ,
365 $ 1 )
366 ELSE
366
367 CALL ZGEBR2D( CONTXT , 'All' , ' ' , 1 , 1 , TJJS , 1 ,
368 $ ITMP1 , ITMP2 )
369 END IF
370 ELSE
370
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
377
378 REC = MIN( ONE , REC*TJJ )
379 USCAL = ZLADIV( USCAL , TJJS )
380 END IF
381 IF( REC.LT.ONE ) THEN
381
382 CALL PZDSCAL( 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 PZDOTC to perform the dot product.
394
394
395 IF( UPPER ) THEN
395
396 CALL PZDOTC( J - 1 , CSUMJ , A , IA , JA + J - 1 , DESCA , 1 ,
397 $ X , IX , JX , DESCX , 1 )
398 ELSE IF( J.LT.N ) THEN
398
399 CALL PZDOTC( 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
402
403 CALL ZGEBS2D( CONTXT , 'Row' , ' ' , 1 , 1 , CSUMJ , 1 )
404 ELSE
404
405 CALL ZGEBR2D( 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
412
413 IF( UPPER ) THEN
414 * DO 180 I = 1 , J - 1
415 * CSUMJ = CSUMJ + ( DCONJG( A( I , J ) )*USCAL )*
416 * $ X( I )
417 * 180 CONTINUE
417
418 ZDUM = DCONJG( USCAL )
419 CALL PZSCAL( J - 1 , ZDUM , A , IA , JA + J - 1 , DESCA , 1 )
420 CALL PZDOTC( J - 1 , CSUMJ , A , IA , JA + J - 1 , DESCA , 1 ,
421 $ X , IX , JX , DESCX , 1 )
422 ZDUM = ZLADIV( CONE , ZDUM )
423 CALL PZSCAL( 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 + ( DCONJG( A( I , J ) )*USCAL )*
427 * $ X( I )
428 * 190 CONTINUE
428
429 ZDUM = DCONJG( USCAL )
430 CALL PZSCAL( N - J , ZDUM , A , IA + J , JA + J - 1 , DESCA , 1 )
431 CALL PZDOTC( N - J , CSUMJ , A , IA + J , JA + J - 1 , DESCA , 1 ,
432 $ X , IX + J , JX , DESCX , 1 )
433 ZDUM = ZLADIV( CONE , ZDUM )
434 CALL PZSCAL( N - J , ZDUM , A , IA + J , JA + J - 1 , DESCA , 1 )
435 END IF
436 IF( MYCOL.EQ.ITMP2X ) THEN
436
437 CALL ZGEBS2D( CONTXT , 'Row' , ' ' , 1 , 1 , CSUMJ , 1 )
438 ELSE
438
439 CALL ZGEBR2D( CONTXT , 'Row' , ' ' , 1 , 1 , CSUMJ , 1 ,
440 $ MYROW , ITMP2X )
441 END IF
442 END IF
443
444 IF( USCAL.EQ.DCMPLX( 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 ) )
450
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 = DCONJG( A( J , J ) )*TSCAL
456
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 ) )
460
461 $ THEN
462 TJJS = DCONJG( A(( ICOL - 1 )*LDA + IROW ) )*TSCAL
463 CALL ZGEBS2D( CONTXT , 'All' , ' ' , 1 , 1 , TJJS ,
464 $ 1 )
465 ELSE
465
466 CALL ZGEBR2D( CONTXT , 'All' , ' ' , 1 , 1 , TJJS , 1 ,
467 $ ITMP1 , ITMP2 )
468 END IF
469 ELSE
469
470 TJJS = TSCAL
471 IF( TSCAL.EQ.ONE )
471
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
481
482 IF( TJJ.LT.ONE ) THEN
482
483 IF( XJ.GT.TJJ*BIGNUM ) THEN
484
485 * Scale X by 1 / abs(x(j)).
486
486
487 REC = ONE / XJ
488 CALL PZDSCAL( 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 ) = ZLADIV( X( J ) , TJJS )
495 XJTMP = ZLADIV( XJTMP , TJJS )
496 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
496
497 $ X( IROWX ) = XJTMP
498 ELSE IF( TJJ.GT.ZERO ) THEN
499
500 * 0 < abs(A(j , j)) <= SMLNUM :
501
501
502 IF( XJ.GT.TJJ*BIGNUM ) THEN
503
504 * Scale x by(1 / abs(x(j)))*abs(A(j , j))*BIGNUM.
505
505
506 REC =( TJJ*BIGNUM ) / XJ
507 CALL PZDSCAL( 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 ) = ZLADIV( X( J ) , TJJS )
513 XJTMP = ZLADIV( XJTMP , TJJS )
514 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
514
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
520
521 CALL PZLASET ( ' ' , N , 1 , CZERO , CZERO , X , IX , JX ,
522 $ DESCX )
523 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
523
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 ) = ZLADIV( X( J ) , TJJS ) - CSUMJ
535
536 XJTMP = ZLADIV( XJTMP , TJJS ) - CSUMJ
537 IF(( MYROW.EQ.ITMP1X ) .AND.( MYCOL.EQ.ITMP2X ) )
537
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
548
549 CALL DSCAL( N , ONE / TSCAL , CNORM , 1 )
550 END IF
551
552 RETURN
553
554 * End of PZLATTRS
555
556 END95
47
|
|
Variables in Routine PZLATTRS()
| Summary Report |
| Data Type | Quantity | Size(byte) |
| CHARACTER | 4 | 4 |
| COMPLEX*16 | 4 | ? |
| DOUBLE PRECISION | 6 | 24 |
| INTEGER | 20 | 80 |
| REAL | 9 | 36 |
| TOTAL | 43 | 144 |
List of Variables
CHARACTER
COMPLEX*16
DOUBLE PRECISION
| CNORM( * ) | HALF | ONE | SCALE | TWO |
| ZERO | | | | |
INTEGER
| BLOCK_CYCLIC_2D | CSRC_ | CTXT_ | DESCA( * ) | DESCX( * ) |
| DLEN_ | DTYPE_ | IA | INFO | IX |
| J | JA | JX | LLD_ | M_ |
| MB_ | N | N_ | NB_ | RSRC_ |
REAL
| CSUMJ | REC | TJJ | TJJS | USCAL |
| XJ | XJTMP | XMAX | ZDUM | |
Variables Dependence Graph
Put the mouse over a right hand side variable to display the corresponding line of the dependence | | - | | - | - | | CSUMJ | <--- | CZEROCSUMJ = CZERO |
| REC | <--- | HALFREC = REC*HALF, ONEREC = ONE / XJ{2REC = ONE / MAX( XMAX, ONE ), 3REC = MIN( ONE, REC*TJJ ), 4REC = ONE / XJ}, RECREC = REC*HALF{2REC = MIN( ONE, REC*TJJ )}, TJJREC = ( TJJ*BIGNUM ) / XJ{2REC = MIN( ONE, REC*TJJ ), 3REC = ( TJJ*BIGNUM ) / XJ}, XJREC = ONE / XJ{2REC = ( TJJ*BIGNUM ) / XJ, 3REC = ONE / XJ, 4REC = ( TJJ*BIGNUM ) / XJ}, XMAXREC = ONE / MAX( XMAX, ONE ) |
| SCALE | <--- | RECSCALE = SCALE*REC{2SCALE = SCALE*REC, 3SCALE = SCALE*REC, 4SCALE = SCALE*REC, 5SCALE = SCALE*REC}, SCALESCALE = SCALE*REC{2SCALE = SCALE*REC, 3SCALE = SCALE*REC, 4SCALE = SCALE*REC, 5SCALE = SCALE*REC, 6SCALE = SCALE / TSCAL}, ZEROSCALE = ZERO{2SCALE = ZERO} |
| TJJ | <--- | TJJSTJJ = CABS1( TJJS ){2TJJ = CABS1( TJJS )} |
| TJJS | <--- | ATJJS = DCONJG( A( ( ICOL-1 )*LDA+IROW ) )*TSCAL{2TJJS = DCONJG( A( ( ICOL-1 )*LDA+IROW ) )*TSCAL} |
| USCAL | <--- | TJJSUSCAL = ZLADIV( USCAL, TJJS ), USCALUSCAL = ZLADIV( USCAL, TJJS ) |
| X | <--- | CONEX( IROWX ) = CONE, XJTMPX( IROWX ) = XJTMP{2X( IROWX ) = XJTMP, 3X( IROWX ) = XJTMP} |
| XJ | <--- | XJTMPXJ = CABS1( XJTMP ){2XJ = CABS1( XJTMP )} |
| XJTMP | <--- | RECXJTMP = XJTMP*REC{2XJTMP = XJTMP*REC, 3XJTMP = XJTMP*REC, 4XJTMP = XJTMP*REC, 5XJTMP = XJTMP*REC}, TJJSXJTMP = ZLADIV( XJTMP, TJJS ){2XJTMP = ZLADIV( XJTMP, TJJS ), 3XJTMP = ZLADIV( XJTMP, TJJS ) - CSUMJ, 4XJTMP = ZLADIV( XJTMP, TJJS ), 5XJTMP = ZLADIV( XJTMP, TJJS ), 6XJTMP = ZLADIV( XJTMP, TJJS ) - CSUMJ}, XXJTMP = X( IROWX ), CONEXJTMP = CONE{2XJTMP = CONE}, XJTMPXJTMP = XJTMP*REC{2XJTMP = ZLADIV( XJTMP, TJJS ), 3XJTMP = XJTMP*REC, 4XJTMP = ZLADIV( XJTMP, TJJS ), 5XJTMP = ZLADIV( XJTMP, TJJS ) - CSUMJ, 6XJTMP = XJTMP*REC, 7XJTMP = XJTMP - CSUMJ, 8XJTMP = XJTMP*REC, 9XJTMP = ZLADIV( XJTMP, TJJS ), 10XJTMP = XJTMP*REC, 11XJTMP = ZLADIV( XJTMP, TJJS ), 12XJTMP = ZLADIV( XJTMP, TJJS ) - CSUMJ}, CSUMJXJTMP = ZLADIV( XJTMP, TJJS ) - CSUMJ{2XJTMP = XJTMP - CSUMJ, 3XJTMP = ZLADIV( XJTMP, TJJS ) - CSUMJ} |
| XMAX | <--- | RECXMAX = XMAX*REC{2XMAX = XMAX*REC, 3XMAX = XMAX*REC, 4XMAX = XMAX*REC, 5XMAX = XMAX*REC}, XJTMPXMAX = MAX( XMAX, CABS1( XJTMP ) ){2XMAX = MAX( XMAX, CABS1( XJTMP ) )}, XMAXXMAX = XMAX*REC{2XMAX = XMAX*REC, 3XMAX = MAX( XMAX, CABS1( XJTMP ) ), 4XMAX = XMAX*REC, 5XMAX = XMAX*REC, 6XMAX = XMAX*REC, 7XMAX = MAX( XMAX, CABS1( XJTMP ) )}, ZEROXMAX = ZERO{2XMAX = ZERO} |
| ZDUM | <--- | USCALZDUM = DCONJG( USCAL ){2ZDUM = DCONJG( USCAL )}, CONEZDUM = ZLADIV( CONE, ZDUM ){2ZDUM = ZLADIV( CONE, ZDUM )}, ZDUMZDUM = ZLADIV( CONE, ZDUM ){2ZDUM = ZLADIV( CONE, ZDUM )} |
|
|
Analysis elements of the routine PZLATTRS()
Put the mouse over each element to display detailed matching information
 Assigned variables |
| | | BLOCK_CYCLIC_2D , CONE , CSRC_ , CSUMJ , CTXT_ , CZERO , DIAG , DLEN_ , DTYPE_ , HALF , INFO , j , LLD_ , M_ , MB_ , N , N_ , NB_ , NORMIN , ONE , REC , RSRC_ , SCALE , TJJ , TJJS , TRANS , TWO , UPLO , USCAL , x , XJ , XJTMP , XMAX , ZDUM , ZERO |
|
| Active variables |
| | | A , BLOCK_CYCLIC_2D , CNORM , CONE , CSRC_ , CSUMJ , CTXT_ , CZERO , DESCA , DESCX , DIAG , DLEN_ , DTYPE_ , HALF , IA , INFO , IX , j , JA , JX , LLD_ , M_ , MB_ , N , N_ , NB_ , NORMIN , one , REC , RSRC_ , SCALE , TJJ , TJJS , TRANS , two , UPLO , USCAL , X , XJ , XJTMP , XMAX , ZDUM , zero |
|
| Accessed arrays [ array name : associated index ] |
| |  A | : ( ICOL-1 )*LDA+IROW , ( ICOL-1 )*LDA+IROW , * , 1:i-1,i , 1:j-1,j , I, J , I, J , i,i , i,i , i,i , J, J , J, J , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j,j , j+1,j+1 , j+1,j+1 , j+1:n,j , j+2:n,j+1 , k,j |
| | CNORM | : * , i , i , i , j , j , j , j , j , J , j+1 , j+1 |
| | DESCA | : * , CSRC_ , CTXT_ , DTYPE_ , LLD_ , LLD_ , M_ , MB_ , N_ , NB_ , RSRC_ |
| | DESCX | : * , LLD_ |
| | X | : * , 1:i-1 , 1:j , 1:j-1 , 1:n , 1:n , 1:n , i , i , I , I , IROWX , IROWX , IROWX , IROWX , IROWX , IROWX , IROWX , IROWX , IROWX , IROWX , j , j , j , j , j , j , j , j , J , j , J , j , j , J , j , j , j , J , J , j , j , J , j , J , j , j , J , j+1:n , j+1:n , k |
|
| Conditional statements [ statement : associated predicate ] |
| | do | : ( 140 J = JFIRST , JLAST , JINC ) , ( it. ) , ( 180 I = 1 , J - 1 ) , ( 190 I = J + 1 , N ) |
| | for | : ( some j) ) , ( any 2D block cyclicly distributed array. ) , ( these quantities may be computed by : ) , ( the distributed matrix A. ) , ( the distributed matrix X. ) , ( the triangular system ) , ( possible ) , ( solving A*x = b. The basic algorithm ) , ( j = 1, ... , n ) , ( iteration j + 1 we have ) , ( A upper triangular is ) , ( j = 1, ... , n ) , ( j >= 1. ) , ( A in the dot product is 1 , ) , ( return. ) |
| | if | : ( the ) , ( the matrix A is singular (A(j , j) = 0 for some j) ) , ( K were distributed over the r processes of its ) , ( K were distributed over the c processes of ) , ( UPLO = 'U' , the leading n by n ) , ( UPLO = 'L' , the leading n by n lower ) , ( DIAG = 'U' , the diagonal elements of A are ) , ( SCALE = 0 , the matrix A is singular or badly scaled , and ) , ( NORMIN = 'Y' , CNORM is an input argument and CNORM(j) ) , ( TRANS = 'N' , CNORM(j) must be greater than or equal ) , ( TRANS = 'T' or 'C' , CNORM(j) ) , ( NORMIN = 'N' , CNORM is an output argument and CNORM(j) ) , ( INFO = - k , the k - th argument had an illegal value ) , ( that is less than overflow , PZTRSV ) , ( A is lower triangular is ) , ( the ) , ( the bound overflows , x is set to 0 , x(j) to ) , ( 1/M(n) and 1 / G(n) are both greater ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( TJJ.GT.ZERO ) , ( XJ.GT.TJJ*BIGNUM ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( the dot ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( (CNORM( J ).GT.( BIGNUM - XJ )*REC ) ) , ( x(j) could overflow , scale x by 1 / (2*XMAX). ) , ( NOUNIT ) , ( (( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) ) , ( TJJ.GT.ONE ) , ( A(j,j) > 1. ) , ( REC.LT.ONE ) , ( USCAL.EQ.CONE ) , ( the scaling needed for A in the dot product is 1 , ) , ( UPPER ) , ( J.LT.N ) , ( MYCOL.EQ.ITMP2X ) , ( UPPER ) , ( J.LT.N ) , ( MYCOL.EQ.ITMP2X ) , ( (USCAL.EQ.DCMPLX( TSCAL ) ) ) , ( 1/A(j,j) ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( NOUNIT ) , ( (( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) ) , ( TSCAL.EQ.ONE ) , ( necessary. ) , ( TJJ.GT.SMLNUM ) , ( TJJ.LT.ONE ) , ( XJ.GT.TJJ*BIGNUM ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( TJJ.GT.ZERO ) , ( XJ.GT.TJJ*BIGNUM ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( the dot ) , ( (( MYROW.EQ.ITMP1X ) .AND. ( MYCOL.EQ.ITMP2X ) ) ) , ( TSCAL.NE.ONE ) |
|
| List of variables |
A( * )
BLOCK_CYCLIC_2D
CNORM( * )
CONE
CSRC_
CSUMJ
CTXT_
| CZERO
DESCA( * )
DESCX( * )
DIAG
DLEN_
DTYPE_
HALF
IA
| INFO
IX
J
JA
JX
LLD_
M_
MB_
| N
N_
NB_
NORMIN
ONE
REC
RSRC_
SCALE
| TJJ
TJJS
TRANS
TWO
UPLO
USCAL
X( * )
XJ
| XJTMP
XMAX
ZDUM
ZERO | | close
| |
A( * )
BLOCK_CYCLIC_2D
CNORM( * )
CONE
CSRC_
CSUMJ
CTXT_
CZERO
DESCA( * )
DESCX( * )
DIAG
DLEN_
DTYPE_
HALF
IA
INFO
IX
J
JA
JX
LLD_
M_
MB_
N
N_
NB_
NORMIN
ONE
REC
RSRC_
SCALE
TJJ
TJJS
TRANS
TWO
UPLO
USCAL
X( * )
XJ
XJTMP
XMAX
ZDUM
ZERO
542
| |
